diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index 33baf960fde029b6b5d8043852edf537f7cc9aad..22ecc8509048a45c158247fc492a77ebb0de0f1d 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -414,6 +414,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand const NodeValue<const bool> node_is_angle = [&] { NodeValue<bool> compute_node_is_angle{mesh->connectivity()}; compute_node_is_angle.fill(false); + // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { if (is_boundary_node[node_id]) { @@ -431,10 +432,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand n2 = nl[face_id]; } } - if (n1 != n2) { - compute_node_is_angle[node_id] = true; - } } + if (l2Norm(n1 - n2) > (1.E-15) and l2Norm(n1 + n2) > (1.E-15) and not node_is_corner[node_id]) { + compute_node_is_angle[node_id] = true; + } + // std::cout << node_id << " " << n1 << " " << n2 << " " << compute_node_is_angle[node_id] << "\n"; } }); return compute_node_is_angle; @@ -461,7 +463,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }(); // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { - // std::cout << node_id << " " << exterior_normal[node_id] << "\n"; + // std::cout << node_id << " " << exterior_normal[node_id] << " " << node_is_angle[node_id] << "\n"; //}; const FaceValue<const double> mes_l = [&] { @@ -920,8 +922,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand FaceId face_id = node_to_face[i_face]; sum_mes_l += mes_l[face_id]; } - b[cell_id] += 1. / Vj[cell_id] * 0.5 * node_boundary_values[node_id] * sum_mes_l; - b_prev[cell_id] += 1. / Vj[cell_id] * 0.5 * node_boundary_prev_values[node_id] * sum_mes_l; + b[cell_id] += 0.5 / Vj[cell_id] * node_boundary_values[node_id] * sum_mes_l; + b_prev[cell_id] += 0.5 / Vj[cell_id] * node_boundary_prev_values[node_id] * sum_mes_l; } } else if (node_is_dirichlet[node_id]) { if (not node_is_corner[node_id]) { @@ -980,6 +982,34 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand // std::cout << "g = " << node_boundary_values << "\n"; + NodeValue<TinyVector<Dimension>> ur{mesh->connectivity()}; + ur.fill(zero); + for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++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); + if (not is_boundary_node[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 size_t i_node = node_local_number_in_its_cell[i_cell]; + ur[node_id] += Pj[cell_id] * Cjr(cell_id, i_node); + } + ur[node_id] = -inverse(node_betar[node_id]) * ur[node_id]; + // std::cout << node_id << "; ur = " << ur[node_id] << "\n"; + } else if (is_boundary_node[node_id] and node_is_dirichlet[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 size_t i_node = node_local_number_in_its_cell[i_cell]; + ur[node_id] += (node_boundary_values[node_id] - Pj[cell_id]) * Cjr(cell_id, i_node); + } + if (not node_is_corner[node_id]) { + ur[node_id] = inverse(node_betar[node_id]) * ur[node_id]; + } else { + ur[node_id] = inverse(corner_betar[node_id]) * ur[node_id]; + } + } + // std::cout << node_id << "; ur = " << ur[node_id] << "\n"; + } + LinearSolver solver; solver.solveLocalSystem(A1, T, B);