diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index c349f4b5e3e3516c10629f6f755c4e38a3e702b9..c834964f7d9db404fff0e5c6e174c716999459d7 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -176,7 +176,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand mesh_data.xl(), mesh_face_boundary .faceList()); - boundary_condition_list.push_back( DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list}); @@ -254,7 +253,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand using T = std::decay_t<decltype(bc)>; if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) { const auto& face_list = bc.faceList(); - for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { const FaceId face_id = face_list[i_face]; const auto& face_nodes = face_to_node_matrix[face_id]; @@ -495,7 +493,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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]; + size_t i_cell = 0; const CellId cell_id = node_to_cell[i_cell]; const unsigned int i_node = node_local_number_in_its_cell[i_cell]; @@ -520,7 +518,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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."); } @@ -721,17 +718,33 @@ 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]; + if (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]; - b[cell_id] += 0; + b[cell_id_j] += dot(node_boundary_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), + Cjr(cell_id_j, i_node_j)); + } } + } 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]; + 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]; - b[cell_id] += 0; + 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]; + + b[cell_id_j] += + dot(node_boundary_values[node_id] * corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), + Cjr(cell_id_j, i_node_j)); + } } } };