diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 941682b453c6cf13cd4774c4ce48968f142dda22..eea42b65b93e20d03a4db2857ba1f2e32447d43f 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -9,6 +9,7 @@
#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/CubeTransformation.hpp>
+#include <geometry/LineCubicTransformation.hpp>
#include <geometry/LineParabolicTransformation.hpp>
#include <geometry/LineTransformation.hpp>
#include <geometry/PrismTransformation.hpp>
@@ -443,10 +444,10 @@ _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
}
}
-template <>
+template <typename ConformTransformationT>
void
_computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
- const LineParabolicTransformation<2>& T,
+ const ConformTransformationT& T,
const TinyVector<2>& Xj,
const double Vi,
const size_t degree,
@@ -541,14 +542,39 @@ _computeEjkMeanByBoundary(const MeshType& mesh,
const FaceId face_id = cell_face_list[i_face];
auto face_node_list = face_to_node_matrix[face_id];
- if (is_reversed) {
- const LineParabolicTransformation<2> T{xr[face_node_list[1]], xl[face_id][0], xr[face_node_list[0]]};
- _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
- } else {
- const LineParabolicTransformation<2> T{xr[face_node_list[0]], xl[face_id][0], xr[face_node_list[1]]};
- _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+#warning rework
+ switch (mesh.degree()) {
+ case 2: {
+ if (is_reversed) {
+ const LineParabolicTransformation<2> T{xr[face_node_list[1]], xl[face_id][0], xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ } else {
+ const LineParabolicTransformation<2> T{xr[face_node_list[0]], xl[face_id][0], xr[face_node_list[1]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ }
+ break;
+ }
+ case 3: {
+ if (is_reversed) {
+ const LineCubicTransformation<2> T{xr[face_node_list[1]], xl[face_id][1], xl[face_id][0],
+ xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ } else {
+ const LineCubicTransformation<2> T{xr[face_node_list[0]], xl[face_id][0], xl[face_id][1],
+ xr[face_node_list[1]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ }
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << "reconstruction on meshes of degree " << mesh.degree();
+ throw NotImplementedError(error_msg.str());
+ }
}
}
}
@@ -615,18 +641,45 @@ _computeEjkMeanByBoundaryInSymmetricCell(const MeshType& mesh,
const FaceId face_id = cell_face_list[i_face];
auto face_node_list = face_to_node_matrix[face_id];
- const auto x0 = symmetrize_coordinates(origin, normal, xr[face_node_list[1]]);
- const auto x1 = symmetrize_coordinates(origin, normal, xl[face_id][0]);
- const auto x2 = symmetrize_coordinates(origin, normal, xr[face_node_list[0]]);
+ switch (mesh.degree()) {
+ case 2: {
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[face_node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xl[face_id][0]);
+ const auto x2 = symmetrize_coordinates(origin, normal, xr[face_node_list[0]]);
- if (is_reversed) {
- const LineParabolicTransformation<2> T{x2, x1, x0};
- _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
- } else {
- const LineParabolicTransformation<2> T{x0, x1, x2};
- _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ if (is_reversed) {
+ const LineParabolicTransformation<2> T{x2, x1, x0};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ } else {
+ const LineParabolicTransformation<2> T{x0, x1, x2};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ }
+ break;
+ }
+ case 3: {
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[face_node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xl[face_id][1]);
+ const auto x2 = symmetrize_coordinates(origin, normal, xl[face_id][0]);
+ const auto x3 = symmetrize_coordinates(origin, normal, xr[face_node_list[0]]);
+
+ if (is_reversed) {
+ const LineCubicTransformation<2> T{x3, x2, x1, x0};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ } else {
+ const LineCubicTransformation<2> T{x0, x1, x2, x3};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ }
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << "reconstruction on meshes of degree " << mesh.degree();
+ throw NotImplementedError(error_msg.str());
+ }
}
}
}
diff --git a/src/scheme/VariationalSolver.cpp b/src/scheme/VariationalSolver.cpp
index 34f946900ae14ceafb9fdc4720cbd492c3157d98..d609cb5c2cc0e1b82e932e63b2e438897296e588 100644
--- a/src/scheme/VariationalSolver.cpp
+++ b/src/scheme/VariationalSolver.cpp
@@ -599,6 +599,10 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
LinearSolver solver;
solver.solveLocalSystem(A, U, b);
+
+ auto error = b - A * U;
+ std::cout << "CG REAL L2 ERROR: " << rang::fgB::yellow << std::sqrt(dot(error, error)) << rang::fg::reset << '\n';
+
this->_forcebcBoundaryConditions(bc_list, U);
const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();