diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index d1aa9d9941305428c3ed9a750ab3512928f9d7d2..90bcff7d4c4db5b3a4908595bf5e9075985620dc 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -305,39 +305,39 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { if (is_boundary_node[node_id]) { - const auto& node_to_cell = node_to_cell_matrix[node_id]; - + const auto& node_to_cell = node_to_cell_matrix[node_id]; + const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id); for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { + const unsigned int i_node = node_local_number_in_its_cell[i_cell]; const CellId cell_id = node_to_cell[i_cell]; const auto& cell_nodes = cell_to_node_matrix[cell_id]; - const NodeId prev_node_id = cell_to_node_matrix[cell_id][(node_id - 1) % cell_nodes.size()]; - const NodeId next_node_id = cell_to_node_matrix[cell_id][(node_id + 1) % cell_nodes.size()]; - + const NodeId prev_node_id = cell_to_node_matrix[cell_id][(i_node - 1) % cell_nodes.size()]; + const NodeId next_node_id = cell_to_node_matrix[cell_id][(i_node + 1) % cell_nodes.size()]; if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) { compute_node_is_corner[node_id] = true; } } } }); - return compute_node_is_corner; }(); const NodeValue<const TinyVector<Dimension>> exterior_normal = [&] { - NodeValue<TinyVector<Dimension>> compute_exterior_normal; - parallel_for( - mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - compute_exterior_normal[node_id] = zero; + NodeValue<TinyVector<Dimension>> compute_exterior_normal{mesh->connectivity()}; + compute_exterior_normal.fill(zero); + for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { + if (is_boundary_node[node_id]) { const auto& node_to_cell = node_to_cell_matrix[node_id]; const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id); for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { - const CellId cell_id = node_to_cell[i_cell]; - const unsigned int i_node_loc = node_local_number_in_its_cell[cell_id]; - compute_exterior_normal[node_id] += Cjr(cell_id, i_node_loc); + const CellId cell_id = node_to_cell[i_cell]; + const unsigned int i_node = node_local_number_in_its_cell[i_cell]; + compute_exterior_normal[node_id] += Cjr(cell_id, i_node); } - compute_exterior_normal[node_id] = - 1. / l2Norm(compute_exterior_normal[node_id]) * compute_exterior_normal[node_id]; - }); + const double norm_exterior_normal = l2Norm(compute_exterior_normal[node_id]); + compute_exterior_normal[node_id] *= 1. / norm_exterior_normal; + } + } return compute_exterior_normal; }(); @@ -352,8 +352,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }(); const NodeValue<const double> node_boundary_values = [&] { - NodeValue<double> compute_node_boundary_values; - NodeValue<double> sum_mes_l; + NodeValue<double> compute_node_boundary_values{mesh->connectivity()}; + NodeValue<double> sum_mes_l{mesh->connectivity()}; compute_node_boundary_values.fill(0); sum_mes_l.fill(0); for (const auto& boundary_condition : boundary_condition_list) { @@ -376,11 +376,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand } } } - for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { - if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id])) { - compute_node_boundary_values[node_id] /= sum_mes_l[node_id]; - } - } + } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) { const auto& face_list = bc.faceList(); const auto& value_list = bc.valueList(); @@ -399,15 +395,17 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand } } } - for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { - if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) { - compute_node_boundary_values[node_id] /= sum_mes_l[node_id]; - } - } } }, boundary_condition); } + for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { + if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id])) { + compute_node_boundary_values[node_id] /= sum_mes_l[node_id]; + } else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) { + compute_node_boundary_values[node_id] /= sum_mes_l[node_id]; + } + } return compute_node_boundary_values; }(); @@ -416,7 +414,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues(); const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] { - NodeValue<TinyMatrix<Dimension>> kappa; + NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()}; parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { kappa[node_id] = zero; @@ -434,7 +432,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }(); const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] { - NodeValue<TinyMatrix<Dimension>> kappa; + NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()}; parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { kappa[node_id] = zero; @@ -452,41 +450,45 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }(); const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] { - NodeValue<TinyMatrix<Dimension>> beta; + NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()}; parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId 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]; beta[node_id] = zero; - for (size_t j = 0; j < node_to_cell.size(); ++j) { - const CellId J = node_to_cell[j]; - const unsigned int R = node_local_number_in_its_cell[j]; - beta[node_id] += tensorProduct(Cjr(J, R), xr[node_id] - xj[J]); + for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { + const CellId cell_id = node_to_cell[i_cell]; + const unsigned int i_node = node_local_number_in_its_cell[i_cell]; + beta[node_id] += tensorProduct(Cjr(cell_id, i_node), xr[node_id] - xj[cell_id]); } }); return beta; }(); const NodeValue<const TinyMatrix<Dimension>> kappar_invBetar = [&] { - NodeValue<TinyMatrix<Dimension>> kappa_invBeta; + NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()}; parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - kappa_invBeta[node_id] = node_kappar[node_id] * inverse(node_betar[node_id]); + if (not node_is_corner[node_id]) { + kappa_invBeta[node_id] = node_kappar[node_id] * inverse(node_betar[node_id]); + } }); return kappa_invBeta; }(); const NodeValue<const TinyMatrix<Dimension>> kapparb_invBetar = [&] { - NodeValue<TinyMatrix<Dimension>> kappa_invBeta; + NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()}; parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(node_betar[node_id]); + if (not node_is_corner[node_id]) { + kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(node_betar[node_id]); + } }); return kappa_invBeta; }(); const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] { - NodeValue<TinyVector<Dimension>> compute_sum_Cjr; + NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()}; parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { const auto& node_to_cell = node_to_cell_matrix[node_id]; @@ -502,7 +504,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }(); const NodeValue<const TinyVector<Dimension>> v = [&] { - NodeValue<TinyVector<Dimension>> compute_v; + NodeValue<TinyVector<Dimension>> compute_v{mesh->connectivity()}; parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { compute_v[node_id] = node_kapparb[node_id] * exterior_normal[node_id]; }); @@ -510,20 +512,23 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }(); const NodeValuePerCell<const double> theta = [&] { - NodeValuePerCell<double> compute_theta; - parallel_for( - mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { - const auto& cell_nodes = cell_to_node_matrix[cell_id]; - for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) { - const NodeId node_id = cell_nodes[i_node]; - compute_theta(cell_id, i_node) = dot(inverse(node_betar[node_id]) * Cjr(cell_id, i_node), v[node_id]); - } - }); + NodeValuePerCell<double> compute_theta{mesh->connectivity()}; + // parallel_for( + // mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + const auto& cell_nodes = cell_to_node_matrix[cell_id]; + for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) { + const NodeId node_id = cell_nodes[i_node]; + std::cout << "node_id = " << node_id << "\n"; + compute_theta(cell_id, i_node) = dot(inverse(node_betar[node_id]) * Cjr(cell_id, i_node), v[node_id]); + } + } + // }); return compute_theta; }(); const NodeValue<const double> sum_theta = [&] { - NodeValue<double> compute_sum_theta; + NodeValue<double> compute_sum_theta{mesh->connectivity()}; parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { const auto& node_to_cell = node_to_cell_matrix[node_id];