diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp index 264ecef042956e6b77230e00f49cd321b50dcc6a..754c6e3ac6448a96e5e05c0a7f60967243b4798f 100644 --- a/src/scheme/FluxingAdvectionSolver.cpp +++ b/src/scheme/FluxingAdvectionSolver.cpp @@ -1,5 +1,8 @@ #include <scheme/FluxingAdvectionSolver.hpp> +#include <analysis/GaussQuadratureDescriptor.hpp> +#include <analysis/QuadratureManager.hpp> +#include <geometry/PrismTransformation.hpp> #include <language/utils/EvaluateAtPoints.hpp> #include <language/utils/InterpolateItemArray.hpp> #include <language/utils/InterpolateItemValue.hpp> @@ -322,6 +325,10 @@ template <> FaceValue<double> FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume() { + // due to the z component of the jacobian determinant, degree 3 + // polynomials must be exactly integrated + const QuadratureFormula<3>& gauss = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(3)); + const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix(); FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity()); auto old_xr = m_old_mesh->xr(); @@ -329,31 +336,10 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume() parallel_for( m_new_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { const auto& face_to_node = face_to_node_matrix[face_id]; - if (face_to_node.size() == 4) { - const Rd& x0 = old_xr[face_to_node[0]]; - const Rd& x1 = old_xr[face_to_node[1]]; - const Rd& x2 = old_xr[face_to_node[2]]; - const Rd& x3 = old_xr[face_to_node[3]]; - - const Rd& x4 = new_xr[face_to_node[0]]; - const Rd& x5 = new_xr[face_to_node[1]]; - const Rd& x6 = new_xr[face_to_node[2]]; - const Rd& x7 = new_xr[face_to_node[3]]; - - const Rd& a1 = x6 - x1; - const Rd& a2 = x6 - x3; - const Rd& a3 = x6 - x4; - - const Rd& b1 = x7 - x0; - const Rd& b2 = x5 - x0; - const Rd& b3 = x2 - x0; - TinyMatrix<3> M1(a1 + b1, a2, a3); - TinyMatrix<3> M2(b1, a2 + b2, a3); - TinyMatrix<3> M3(a1, b2, a3 + b3); + const size_t face_nb_node = face_to_node.size(); - algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12; - } else if (face_to_node.size() == 3) { + if (face_nb_node == 3) { const Rd& x0 = old_xr[face_to_node[0]]; const Rd& x1 = old_xr[face_to_node[1]]; const Rd& x2 = old_xr[face_to_node[2]]; @@ -362,21 +348,40 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume() const Rd& x4 = new_xr[face_to_node[1]]; const Rd& x5 = new_xr[face_to_node[2]]; - const Rd& a1 = x5 - x1; - const Rd& a2 = x5 - x2; - const Rd& a3 = x5 - x3; + double volume = 0; - const Rd& b1 = x5 - x0; - const Rd& b2 = x4 - x0; - const Rd& b3 = x2 - x0; - - TinyMatrix<3> M1(a1 + b1, a2, a3); - TinyMatrix<3> M2(b1, a2 + b2, a3); - TinyMatrix<3> M3(a1, b2, a3 + b3); + PrismTransformation T(x0, x1, x2, x3, x4, x5); + for (size_t i = 0; i < gauss.numberOfPoints(); ++i) { + volume += gauss.weight(i) * T.jacobianDeterminant(gauss.point(i)); + } - algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12; + algebraic_fluxing_volume[face_id] = volume; } else { - throw NotImplementedError("Not implemented for non quad faces"); + Rd xg_old = zero; + for (size_t i_node = 0; i_node < face_nb_node; ++i_node) { + xg_old += old_xr[face_to_node[i_node]]; + } + xg_old *= 1. / face_nb_node; + Rd xg_new = zero; + for (size_t i_node = 0; i_node < face_nb_node; ++i_node) { + xg_new += new_xr[face_to_node[i_node]]; + } + xg_new *= 1. / face_nb_node; + + double volume = 0; + for (size_t i_node = 0; i_node < face_nb_node; ++i_node) { + PrismTransformation T(xg_old, // + old_xr[face_to_node[i_node]], // + old_xr[face_to_node[(i_node + 1) % face_nb_node]], // + xg_new, // + new_xr[face_to_node[i_node]], // + new_xr[face_to_node[(i_node + 1) % face_nb_node]]); + for (size_t i = 0; i < gauss.numberOfPoints(); ++i) { + volume += gauss.weight(i) * T.jacobianDeterminant(gauss.point(i)); + } + } + + algebraic_fluxing_volume[face_id] = volume; } });