diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp index 84ec83e474941185aa06420b4e536f70c2c4f539..38ef1f5e86bceff24452643f35e39e1d2f5743e6 100644 --- a/src/scheme/ScalarHybridScheme.cpp +++ b/src/scheme/ScalarHybridScheme.cpp @@ -807,33 +807,42 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH return compute_w; }(); - const NodeValue<const TinyVector<2>> normal_edge_index = [&] { - NodeValue<TinyVector<2>> compute_normal_edge_index{mesh->connectivity()}; + const NodeValuePerCell<const TinyVector<2>> normal_edge_index = [&] { + NodeValuePerCell<TinyVector<2>> compute_normal_edge_index{mesh->connectivity()}; compute_normal_edge_index.fill(zero); parallel_for( - mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - bool l1_OK = false; - - if (node_is_edge[node_id]) { - const auto& node_to_cell = node_to_cell_matrix[node_id]; - CellId cell_id = node_to_cell[0]; - const auto& cell_to_node = cell_to_node_matrix[cell_id]; - for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { - const NodeId node_id_in_cell = cell_to_node[i_node]; - if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) { - cell_id = node_to_cell[1]; - } - } - - const auto& cell_to_face = cell_to_face_matrix[cell_id]; - for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) { - const FaceId face_id = cell_to_face[i_face]; - if (is_boundary_face[face_id]) { - if (not l1_OK) { - compute_normal_edge_index[node_id][0] = face_id; - l1_OK = true; - } else { - compute_normal_edge_index[node_id][1] = face_id; + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { + const auto& cell_to_node = cell_to_node_matrix[cell_id]; + for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { + const NodeId node_id = cell_to_node[i_node]; + bool l1_OK = false; + + if (node_is_edge[node_id]) { + // const auto& node_to_cell = node_to_cell_matrix[node_id]; + // CellId cell_id = node_to_cell[0]; + // const auto& cell_to_node = cell_to_node_matrix[cell_id]; + // for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { + // const NodeId node_id_in_cell = cell_to_node[i_node]; + // if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) { + // cell_id = node_to_cell[1]; + // } + // } + + const auto& cell_to_face = cell_to_face_matrix[cell_id]; + for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) { + const FaceId face_id = cell_to_face[i_face]; + const auto& face_to_node = face_to_node_matrix[face_id]; + for (size_t j_node = 0; j_node < face_to_node.size(); ++j_node) { + if (face_to_node[j_node] == node_id) { + if (is_boundary_face[face_id]) { + if (not l1_OK) { + compute_normal_edge_index(cell_id, i_node)[0] = face_id; + l1_OK = true; + } else { + compute_normal_edge_index(cell_id, i_node)[1] = face_id; + } + } + } } } } @@ -842,113 +851,89 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH return compute_normal_edge_index; }(); - const NodeValue<const TinyVector<Dimension>> normal1_edge = [&] { - NodeValue<TinyVector<Dimension>> compute_normal1_edge{mesh->connectivity()}; + const NodeValuePerCell<const TinyVector<Dimension>> normal1_edge = [&] { + NodeValuePerCell<TinyVector<Dimension>> compute_normal1_edge{mesh->connectivity()}; compute_normal1_edge.fill(zero); parallel_for( - mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - if (node_is_edge[node_id]) { - const FaceId l1 = normal_edge_index[node_id][0]; - - const auto& node_to_cell = node_to_cell_matrix[node_id]; - CellId cell_id = node_to_cell[0]; - const auto& cell_to_node = cell_to_node_matrix[cell_id]; - for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { - const NodeId node_id_in_cell = cell_to_node[i_node]; - if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) { - cell_id = node_to_cell[1]; + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { + const auto& cell_to_node = cell_to_node_matrix[cell_id]; + for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { + const NodeId node_id = cell_to_node[i_node]; + if (node_is_edge[node_id]) { + const FaceId l1 = normal_edge_index(cell_id, i_node)[0]; + + // const auto& node_to_cell = node_to_cell_matrix[node_id]; + // CellId cell_id = node_to_cell[0]; + // const auto& cell_to_node = cell_to_node_matrix[cell_id]; + // for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { + // const NodeId node_id_in_cell = cell_to_node[i_node]; + // if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) { + // cell_id = node_to_cell[1]; + // } + // } + + if (dot(nl[l1], xl[l1] - xj[cell_id]) <= 0) { + compute_normal1_edge(cell_id, i_node) = -nl[l1]; + } else { + compute_normal1_edge(cell_id, i_node) = nl[l1]; } } - - if (dot(nl[l1], xl[l1] - xj[cell_id]) <= 0) { - compute_normal1_edge[node_id] = -nl[l1]; - } else { - compute_normal1_edge[node_id] = nl[l1]; - } } }); return compute_normal1_edge; }(); - const NodeValue<const TinyVector<Dimension>> normal2_edge = [&] { - NodeValue<TinyVector<Dimension>> compute_normal2_edge{mesh->connectivity()}; + const NodeValuePerCell<const TinyVector<Dimension>> normal2_edge = [&] { + NodeValuePerCell<TinyVector<Dimension>> compute_normal2_edge{mesh->connectivity()}; compute_normal2_edge.fill(zero); parallel_for( - mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - if (node_is_edge[node_id]) { - const FaceId l2 = normal_edge_index[node_id][1]; - - const auto& node_to_cell = node_to_cell_matrix[node_id]; - CellId cell_id = node_to_cell[0]; - const auto& cell_to_node = cell_to_node_matrix[cell_id]; - for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { - const NodeId node_id_in_cell = cell_to_node[i_node]; - if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) { - cell_id = node_to_cell[1]; + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { + const auto& cell_to_node = cell_to_node_matrix[cell_id]; + for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { + const NodeId node_id = cell_to_node[i_node]; + if (node_is_edge[node_id]) { + const FaceId l2 = normal_edge_index(cell_id, i_node)[1]; + + // const auto& node_to_cell = node_to_cell_matrix[node_id]; + // CellId cell_id = node_to_cell[0]; + // const auto& cell_to_node = cell_to_node_matrix[cell_id]; + // for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { + // const NodeId node_id_in_cell = cell_to_node[i_node]; + // if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) { + // cell_id = node_to_cell[1]; + // } + // } + if (dot(nl[l2], xl[l2] - xj[cell_id]) <= 0) { + compute_normal2_edge(cell_id, i_node) = -nl[l2]; + } else { + compute_normal2_edge(cell_id, i_node) = nl[l2]; } } - if (dot(nl[l2], xl[l2] - xj[cell_id]) <= 0) { - compute_normal2_edge[node_id] = -nl[l2]; - } else { - compute_normal2_edge[node_id] = nl[l2]; - } } }); return compute_normal2_edge; }(); - const NodeValue<const TinyVector<Dimension>> normal_orth_edge = [&] { - NodeValue<TinyVector<Dimension>> compute_normal_orth_edge{mesh->connectivity()}; + const NodeValuePerCell<const TinyVector<Dimension>> normal_orth_edge = [&] { + NodeValuePerCell<TinyVector<Dimension>> compute_normal_orth_edge{mesh->connectivity()}; compute_normal_orth_edge.fill(zero); parallel_for( - mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - // bool l1_OK = false; - // FaceId l1 = 0; - // FaceId l2 = 0; - // TinyVector<Dimension> n1 = zero; - // TinyVector<Dimension> n2 = zero; - - if (node_is_edge[node_id]) { - // const auto& node_to_cell = node_to_cell_matrix[node_id]; - // CellId cell_id = node_to_cell[0]; - // const auto& cell_to_node = cell_to_node_matrix[cell_id]; - // for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { - // const NodeId node_id_in_cell = cell_to_node[i_node]; - // if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) { - // cell_id = node_to_cell[1]; - // } - // } - - // const auto& cell_to_face = cell_to_face_matrix[cell_id]; - // for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) { - // const FaceId face_id = cell_to_face[i_face]; - // if (is_boundary_face[face_id]) { - // if (not l1_OK) { - // l1 = face_id; - // if (dot(nl[l1], xl[l1] - xj[cell_id]) <= 0) { - // n1 = -nl[l1]; - // } else { - // n1 = nl[l1]; - // } - // l1_OK = true; - - // } else { - // l2 = face_id; - // if (dot(nl[l2], xl[l2] - xj[cell_id]) <= 0) { - // n2 = -nl[l2]; - // } else { - // n2 = nl[l2]; - // } - // } - // } - // } - - compute_normal_orth_edge[node_id][0] = - normal1_edge[node_id][1] * normal2_edge[node_id][2] - normal1_edge[node_id][2] * normal2_edge[node_id][1]; - compute_normal_orth_edge[node_id][1] = - normal1_edge[node_id][2] * normal2_edge[node_id][0] - normal1_edge[node_id][0] * normal2_edge[node_id][2]; - compute_normal_orth_edge[node_id][2] = - normal1_edge[node_id][0] * normal2_edge[node_id][1] - normal1_edge[node_id][1] * normal2_edge[node_id][0]; + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { + const auto& cell_to_node = cell_to_node_matrix[cell_id]; + for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { + const NodeId node_id = cell_to_node[i_node]; + + if (node_is_edge[node_id]) { + compute_normal_orth_edge(cell_id, i_node)[0] = + normal1_edge(cell_id, i_node)[1] * normal2_edge(cell_id, i_node)[2] - + normal1_edge(cell_id, i_node)[2] * normal2_edge(cell_id, i_node)[1]; + compute_normal_orth_edge(cell_id, i_node)[1] = + normal1_edge(cell_id, i_node)[2] * normal2_edge(cell_id, i_node)[0] - + normal1_edge(cell_id, i_node)[0] * normal2_edge(cell_id, i_node)[2]; + compute_normal_orth_edge(cell_id, i_node)[2] = + normal1_edge(cell_id, i_node)[0] * normal2_edge(cell_id, i_node)[1] - + normal1_edge(cell_id, i_node)[1] * normal2_edge(cell_id, i_node)[0]; + } } }); return compute_normal_orth_edge; @@ -1432,7 +1417,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH const auto& face_to_node = face_to_node_matrix[face_id]; const size_t numberOfNodesPerFace = face_to_node.size(); if (not is_boundary_face[face_id]) { - TinyVector<Dimension> t = normal_orth_edge[node_id]; + TinyVector<Dimension> t = normal_orth_edge(cell_id_j, i_node_j); for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) { if (face_to_node[i_face_to_node] == node_id) { @@ -1440,18 +1425,12 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) { CellId cell_id = face_to_cell[i_cell]; if (cell_id_k == cell_id_j) { - S1(cell_id_j, cell_id_k) += (1. - lambda) * dt * (1. - cn_coeff / 2.) * - dot(1. / - (dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t) * - pow(l2Norm(node_kappar[node_id] * t), 2)) * - node_kappar[node_id] * t, - node_kappar[node_id] * Cjr(cell_id_j, i_node_j)); - S2(cell_id_j, cell_id_k) -= (1. - lambda) * dt * (cn_coeff / 2.) * - dot(1. / - (dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t) * - pow(l2Norm(node_kappar[node_id] * t), 2)) * - node_kappar[node_id] * t, - node_kappar[node_id] * Cjr(cell_id_j, i_node_j)); + S1(cell_id_j, cell_id_k) += (1. - lambda) * dt * (1. - cn_coeff / 2.) * mes_l[face_id] / + numberOfNodesPerFace * pow(l2Norm(node_kappar[node_id] * t), 2) / + dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t); + S2(cell_id_j, cell_id_k) -= (1. - lambda) * dt * (cn_coeff / 2.) * mes_l[face_id] / + numberOfNodesPerFace * pow(l2Norm(node_kappar[node_id] * t), 2) / + dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t); S1(cell_id_j, cell_id_k) += lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] / numberOfNodesPerFace * @@ -1468,18 +1447,12 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH node_kappar[node_id] * t, base_xj1_xj2_orth(cell_id_j, i_face)); } else { - S1(cell_id_j, cell_id_k) -= (1. - lambda) * dt * (1. - cn_coeff / 2.) * - dot(1. / - (dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t) * - pow(l2Norm(node_kappar[node_id] * t), 2)) * - node_kappar[node_id] * t, - node_kappar[node_id] * Cjr(cell_id_j, i_node_j)); - S2(cell_id_j, cell_id_k) += (1. - lambda) * dt * (cn_coeff / 2.) * - dot(1. / - (dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t) * - pow(l2Norm(node_kappar[node_id] * t), 2)) * - node_kappar[node_id] * t, - node_kappar[node_id] * Cjr(cell_id_j, i_node_j)); + S1(cell_id_j, cell_id_k) -= (1. - lambda) * dt * (1. - cn_coeff / 2.) * mes_l[face_id] / + numberOfNodesPerFace * pow(l2Norm(node_kappar[node_id] * t), 2) / + dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t); + S2(cell_id_j, cell_id_k) += (1. - lambda) * dt * (cn_coeff / 2.) * mes_l[face_id] / + numberOfNodesPerFace * pow(l2Norm(node_kappar[node_id] * t), 2) / + dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t); S1(cell_id_j, cell_id_k) -= lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] / numberOfNodesPerFace * @@ -1725,50 +1698,41 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH const auto& cell_to_face = cell_to_face_matrix[cell_id]; for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) { - const FaceId face_id = cell_to_face[i_face]; - const auto& face_to_node = face_to_node_matrix[face_id]; - const size_t numberOfNodesPerFace = face_to_node.size(); - if (not is_boundary_face[face_id]) { - FaceId l1 = normal_edge_index[node_id][0]; - FaceId l2 = normal_edge_index[node_id][1]; - TinyVector<Dimension> n1 = normal1_edge[node_id]; - TinyVector<Dimension> n2 = normal2_edge[node_id]; - TinyVector<Dimension> t = normal_orth_edge[node_id]; + const FaceId face_id = cell_to_face[i_face]; + const auto& face_to_node = face_to_node_matrix[face_id]; - // std::cout << "node_id : " << node_id << "; cell_id : " << cell_id - // << "; xj1-xj2 = " << xj1_xj2(cell_id, i_face) << "; t = " << t << "\n"; + if (not is_boundary_face[face_id]) { + const FaceId l1 = normal_edge_index(cell_id, i_node)[0]; + const FaceId l2 = normal_edge_index(cell_id, i_node)[1]; + const TinyVector<Dimension> n1 = normal1_edge(cell_id, i_node); + const TinyVector<Dimension> n2 = normal2_edge(cell_id, i_node); + const TinyVector<Dimension> t = normal_orth_edge(cell_id, i_node); + const size_t numberOfNodesl1 = face_to_node_matrix[l1].size(); + const size_t numberOfNodesl2 = face_to_node_matrix[l2].size(); + const size_t numberOfNodesl3 = face_to_node_matrix[face_id].size(); for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) { if (face_to_node[i_face_to_node] == node_id) { - b[cell_id] += (1. - lambda) * - dot(face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * - node_kappar[node_id] * n1 + - face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) * - node_kappar[node_id] * n2 - - 1. / dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) * - dot(xj1_xj2(cell_id, i_face), - face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * - node_kappar[node_id] * n1 + - face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) * - node_kappar[node_id] * n2) / - pow(l2Norm(node_kappar[node_id] * t), 2) * node_kappar[node_id] * t, - node_kappar[node_id] * Cjr(cell_id, i_node)); - b_prev[cell_id] += - (1. - lambda) * - dot(face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * - node_kappar[node_id] * n1 + - face_boundary_prev_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) * - node_kappar[node_id] * n2 - - 1. / dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) * - dot(xj1_xj2(cell_id, i_face), - face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * - node_kappar[node_id] * n1 + - face_boundary_prev_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) * - node_kappar[node_id] * n2) / - pow(l2Norm(node_kappar[node_id] * t), 2) * node_kappar[node_id] * t, - node_kappar[node_id] * Cjr(cell_id, i_node)); - - b[cell_id] += lambda * mes_l[face_id] / numberOfNodesPerFace * + b[cell_id] += (1. - lambda) * mes_l[l1] * face_boundary_values[l1] / numberOfNodesl1 + + mes_l[l2] * face_boundary_values[l2] / numberOfNodesl2 - + mes_l[face_id] / numberOfNodesl3 * pow(l2Norm(node_kappar[node_id] * t), 2) / + dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) * + dot(xj1_xj2(cell_id, i_face), + face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * + node_kappar[node_id] * n1 + + face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) * + node_kappar[node_id] * n2); + b_prev[cell_id] += (1. - lambda) * mes_l[l1] * face_boundary_prev_values[l1] / numberOfNodesl1 + + mes_l[l2] * face_boundary_prev_values[l2] / numberOfNodesl2 - + mes_l[face_id] / numberOfNodesl3 * pow(l2Norm(node_kappar[node_id] * t), 2) / + dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) * + dot(xj1_xj2(cell_id, i_face), + face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * + node_kappar[node_id] * n1 + + face_boundary_prev_values[l2] / + pow(l2Norm(node_kappar[node_id] * n2), 2) * node_kappar[node_id] * n2); + + b[cell_id] += lambda * mes_l[face_id] / numberOfNodesl3 * dot(face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * node_kappar[node_id] * n1 + face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) * @@ -1782,7 +1746,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH pow(l2Norm(node_kappar[node_id] * t), 2) * node_kappar[node_id] * t, base_xj1_xj2_orth(cell_id, i_face)); b_prev[cell_id] += - lambda * mes_l[face_id] / numberOfNodesPerFace * + lambda * mes_l[face_id] / numberOfNodesl3 * dot(face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * node_kappar[node_id] * n1 + face_boundary_prev_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *