diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index 4b65e08eb4b167e0c1351c68b47c559c4de851f2..8246e2121a55d8ddefde348b6dc2527a9a1b2668 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -518,7 +518,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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]; - if (not is_boundary_node[node_id]) { + if (is_boundary_node[node_id] && not node_is_corner[node_id]) { compute_theta(cell_id, i_node) = dot(inverse(node_betar[node_id]) * Cjr(cell_id, i_node), v[node_id]); } } @@ -530,13 +530,15 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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]; - const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id); - compute_sum_theta[node_id] = 0; - for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { - const CellId cell_id = node_to_cell[i_cell]; - const size_t i_node = node_local_number_in_its_cell[i_cell]; - compute_sum_theta[node_id] += theta(cell_id, i_node); + if ((is_boundary_node[node_id]) && (not node_is_corner[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); + compute_sum_theta[node_id] = 0; + for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { + const CellId cell_id = node_to_cell[i_cell]; + const size_t i_node = node_local_number_in_its_cell[i_cell]; + compute_sum_theta[node_id] += theta(cell_id, i_node); + } } }); return compute_sum_theta; @@ -627,6 +629,20 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand } b[cell_id] += 0.5 * node_boundary_values[node_id] * sum_mes_l; } + } else if (node_is_dirichlet[node_id]) { + if (is_boundary_node[node_id] && (not node_is_corner[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]; + + b[cell_id] += 0; + } + } else if (node_is_corner[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]; + + b[cell_id] += 0; + } + } }; Vector<double> T{mesh->numberOfCells()};