diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp index 55c38578f5a025327d59b4f380f56dc9583adb38..61285fab57d4cf8d3fb4b8f38c0d6cb117eab015 100644 --- a/src/scheme/ScalarHybridScheme.cpp +++ b/src/scheme/ScalarHybridScheme.cpp @@ -693,8 +693,15 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH face_normal = -nl[face_id]; } + // compute_normal_base(cell_id, i_face)[0] = + // dot(face_kappal[face_id] * face_normal, 1. / l2Norm(xj1_xj2(cell_id, i_face)) * xj1_xj2(cell_id, + // i_face)); + // compute_normal_base(cell_id, i_face)[1] = + // dot(face_kappal[face_id] * face_normal, + // 1. / l2Norm(xj1_xj2_orth(cell_id, i_face)) * xj1_xj2_orth(cell_id, i_face)); compute_normal_base(cell_id, i_face) = inverse(A) * (face_kappal[face_id] * face_normal); - // std::cout << "alpha_l = " << compute_normal_base(cell_id, i_face)[0] + + // std::cout << cell_id << " " << face_id << "; alpha_l = " << compute_normal_base(cell_id, i_face)[0] // << "; delta_l = " << compute_normal_base(cell_id, i_face)[1] << "\n"; } } @@ -1074,11 +1081,19 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH if (face_to_node[0] == node_id or face_to_node[1] == node_id) { S1(cell_id_j, cell_id_k) += lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * - dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j)); + dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face)); S2(cell_id_j, cell_id_k) -= lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * - dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j)); - // std::cout << "alpha_l = " << normal_face_base(cell_id_j, i_face)[0] << "\n"; + dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face)); + // if (normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) + + // normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) != + // nl[face_id]) { + // std::cout << "(xj1,xj2) = (" << xj[cell_id_k] << "," << xj[cell_id_j] << ") \n" + // << normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) + + // normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) + // << " = " << nl[face_id] << "\n" + // << "\n"; + // } // std::cout << S1(cell_id_j, cell_id_k) << "\n"; } } @@ -1110,7 +1125,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH Q * Cjr(cell_id_j, 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]; @@ -1136,15 +1150,23 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH if (face_to_node[0] == node_id or face_to_node[1] == node_id) { S1(cell_id_j, cell_id_k) += lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * - dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j)); + dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face)); S2(cell_id_j, cell_id_k) -= lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * - dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j)); + dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face)); + // if (normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) + + // normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) != + // nl[face_id]) { + // std::cout << "(xj1,xj2) = (" << xj[cell_id_k] << "," << xj[cell_id_j] << ") \n" + // << normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) + + // normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) + // << " = " << nl[face_id] << "\n" + // << "\n"; + // } } } } } - } 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]; @@ -1170,11 +1192,21 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH if (face_to_node[0] == node_id or face_to_node[1] == node_id) { S1(cell_id_j, cell_id_k) += lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * - dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j)); + dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face)); S2(cell_id_j, cell_id_k) -= lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * - dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j)); + dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face)); + + // if (normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) + + // normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) != + // nl[face_id]) { + // std::cout << "(xj1,xj2) = (" << xj[cell_id_k] << "," << xj[cell_id_j] << ") \n" + // << normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) + + // normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) + // << " = " << nl[face_id] << "\n" + // << "\n"; + // } } } } @@ -1310,15 +1342,14 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH const auto& face_to_node = face_to_node_matrix[face_id]; if (is_boundary_face[face_id]) { if (face_to_node[0] == node_id or face_to_node[1] == node_id) { - b[cell_id_j] += - lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * - node_boundary_values[node_id] * - dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j)); + b[cell_id_j] += lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * + node_boundary_values[node_id] * + dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_face)); // std::cout << "xl - xj2 = " << xl[face_id] - b_prev[cell_id_j] += lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * node_boundary_prev_values[node_id] * - dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j)); + dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_face)); } } } @@ -1350,11 +1381,11 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH b[cell_id_j] += lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * node_boundary_values[node_id] * - dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j)); + dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_face)); b_prev[cell_id_j] += lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * node_boundary_prev_values[node_id] * - dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j)); + dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_face)); } } }