diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index 8246e2121a55d8ddefde348b6dc2527a9a1b2668..c349f4b5e3e3516c10629f6f755c4e38a3e702b9 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -420,10 +420,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand kappa[node_id] = zero; const auto& node_to_cell = node_to_cell_matrix[node_id]; double weight = 0; - for (size_t j = 0; j < node_to_cell.size(); ++j) { - const CellId J = node_to_cell[j]; - double local_weight = 1. / l2Norm(xr[node_id] - xj[J]); - kappa[node_id] += local_weight * cell_kappaj[J]; + for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { + const CellId cell_id = node_to_cell[i_cell]; + double local_weight = 1. / l2Norm(xr[node_id] - xj[cell_id]); + kappa[node_id] += local_weight * cell_kappaj[cell_id]; weight += local_weight; } kappa[node_id] = 1. / weight * kappa[node_id]; @@ -438,10 +438,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand kappa[node_id] = zero; const auto& node_to_cell = node_to_cell_matrix[node_id]; double weight = 0; - for (size_t j = 0; j < node_to_cell.size(); ++j) { - const CellId J = node_to_cell[j]; - double local_weight = 1. / l2Norm(xr[node_id] - xj[J]); - kappa[node_id] += local_weight * cell_kappajb[J]; + for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { + const CellId cell_id = node_to_cell[i_cell]; + double local_weight = 1. / l2Norm(xr[node_id] - xj[cell_id]); + kappa[node_id] += local_weight * cell_kappajb[cell_id]; weight += local_weight; } kappa[node_id] = 1. / weight * kappa[node_id]; @@ -487,6 +487,71 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand return kappa_invBeta; }(); + const NodeValue<const TinyMatrix<Dimension>> corner_betar = [&] { + NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()}; + parallel_for( + mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { + if (node_is_corner[node_id]) { + const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id); + const auto& node_to_cell = node_to_cell_matrix[node_id]; + + size_t i_cell = node_to_cell[0]; + const CellId cell_id = node_to_cell[i_cell]; + const unsigned int i_node = node_local_number_in_its_cell[i_cell]; + + const auto& cell_nodes = cell_to_node_matrix[cell_id]; + const NodeId node_id_p1 = cell_to_node_matrix[cell_id][(i_node + 1) % cell_nodes.size()]; + const NodeId node_id_p2 = cell_to_node_matrix[cell_id][(i_node + 2) % cell_nodes.size()]; + const NodeId node_id_m1 = cell_to_node_matrix[cell_id][(i_node - 1) % cell_nodes.size()]; + + const TinyVector<Dimension> xj1 = 1. / 3. * (xr[node_id] + xr[node_id_m1] + xr[node_id_p2]); + const TinyVector<Dimension> xj2 = 1. / 3. * (xr[node_id] + xr[node_id_p2] + xr[node_id_p1]); + + const TinyVector<Dimension> xrm1 = xr[node_id_m1]; + const TinyVector<Dimension> xrp1 = xr[node_id_p1]; + const TinyVector<Dimension> xrp2 = xr[node_id_p2]; + + TinyVector<Dimension> Cjr1; + TinyVector<Dimension> Cjr2; + + if (Dimension == 2) { + Cjr1[0] = -0.5 * (xrm1[1] - xrp2[1]); + Cjr1[1] = 0.5 * (xrm1[0] - xrp2[0]); + Cjr2[0] = -0.5 * (xrp2[1] - xrp1[1]); + Cjr2[1] = 0.5 * (xrp2[0] - xrp1[0]); + } else { + std::cout << "Dimension " << Dimension << " not implemented. \n"; + throw NotImplementedError("The scheme is not implemented in this dimension."); + } + + beta[node_id] = tensorProduct(Cjr1, (xr[node_id] - xj1)) + tensorProduct(Cjr2, (xr[node_id] - xj2)); + } + }); + return beta; + }(); + + const NodeValue<const TinyMatrix<Dimension>> corner_kappar_invBetar = [&] { + NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()}; + parallel_for( + mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { + if (node_is_corner[node_id]) { + kappa_invBeta[node_id] = node_kappar[node_id] * inverse(corner_betar[node_id]); + } + }); + return kappa_invBeta; + }(); + + const NodeValue<const TinyMatrix<Dimension>> corner_kapparb_invBetar = [&] { + NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()}; + parallel_for( + mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { + if (node_is_corner[node_id]) { + kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(corner_betar[node_id]); + } + }); + return kappa_invBeta; + }(); + const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] { NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()}; parallel_for( @@ -564,7 +629,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); } } - } else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id])) { + } else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id]) && (not node_is_dirichlet[node_id])) { const TinyMatrix<Dimension> I = identity; const TinyVector<Dimension> n = exterior_normal[node_id]; const TinyMatrix<Dimension> nxn = tensorProduct(n, n); @@ -584,6 +649,32 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand P * Cjr(cell_id_k, i_node_j)); } } + } else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) { + for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) { + const CellId cell_id_j = node_to_cell[i_cell_j]; + const size_t i_node_j = node_local_number_in_its_cells[i_cell_j]; + + for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) { + const CellId cell_id_k = node_to_cell[i_cell_k]; + const size_t i_node_k = node_local_number_in_its_cells[i_cell_k]; + + S(cell_id_j, cell_id_k) += + dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); + } + } + } else if (node_is_dirichlet[node_id] && node_is_corner[node_id]) { + for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) { + const CellId cell_id_j = node_to_cell[i_cell_j]; + const size_t i_node_j = node_local_number_in_its_cells[i_cell_j]; + + for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) { + const CellId cell_id_k = node_to_cell[i_cell_k]; + const size_t i_node_k = node_local_number_in_its_cells[i_cell_k]; + + S(cell_id_j, cell_id_k) += + dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); + } + } } }