Use a T4 compatible formula to compute fluxing volumes

Also implement fluxing volume calculation for general faces (polygonal)

src/scheme/FluxingAdvectionSolver.cpp

 #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 size_t face_nb_node = face_to_node.size();
 
-        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);
-
-        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;
-
-        const Rd& b1 = x5 - x0;
-        const Rd& b2 = x4 - x0;
-        const Rd& b3 = x2 - x0;
+        double volume = 0;
 
-        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;
       }
     });