From 53c9a20c2461150471d5164a83d4e860ad3ead64 Mon Sep 17 00:00:00 2001
From: MARMAJOU ISABELLE <id.counilh@wanadoo.fr>
Date: Mon, 10 Feb 2025 17:32:35 +0100
Subject: [PATCH] Add quadrature for polynomial Q3 and P3 cells

---
 src/geometry/LineTransformation.hpp     |   7 +
 src/geometry/SquareQ2Transformation.hpp |   2 +-
 src/geometry/SquareQ3Transformation.hpp |   8 +-
 src/language/utils/IntegrateOnCells.hpp | 361 ++++++++++++++++++------
 src/mesh/MeshData.cpp                   | 196 +++++++++++++
 src/mesh/MeshData.hpp                   | 248 +---------------
 src/scheme/VariationalSolver.cpp        |   4 +-
 7 files changed, 491 insertions(+), 335 deletions(-)

diff --git a/src/geometry/LineTransformation.hpp b/src/geometry/LineTransformation.hpp
index c9f7def71..9e6ffa09e 100644
--- a/src/geometry/LineTransformation.hpp
+++ b/src/geometry/LineTransformation.hpp
@@ -66,6 +66,13 @@ class LineTransformation
     return m_velocity;
   }
 
+  PUGS_INLINE
+  const TinyVector<Dimension>&
+  velocity(const TinyVector<1>&) const
+  {
+    return m_velocity;
+  }
+
   double
   velocityNorm() const
   {
diff --git a/src/geometry/SquareQ2Transformation.hpp b/src/geometry/SquareQ2Transformation.hpp
index 49c0810f8..c5b78f2b4 100644
--- a/src/geometry/SquareQ2Transformation.hpp
+++ b/src/geometry/SquareQ2Transformation.hpp
@@ -15,7 +15,7 @@ class SquareQ2Transformation<2>
   constexpr static size_t NumberOfPoints = 9;
 
  private:
-  // Points are defined as follows (image of the unit triangle with a right angle in a0)
+  // Points are defined as follows (image of the unit square)
   //
   //  a02--a12--a22
   //   |         |
diff --git a/src/geometry/SquareQ3Transformation.hpp b/src/geometry/SquareQ3Transformation.hpp
index 108498b8f..f28fe53be 100644
--- a/src/geometry/SquareQ3Transformation.hpp
+++ b/src/geometry/SquareQ3Transformation.hpp
@@ -15,13 +15,13 @@ class SquareQ3Transformation<2>
   constexpr static size_t NumberOfPoints = 16;
 
  private:
-  // Points are defined as follows (image of the unit triangle with a right angle in a0)
+  // Points are defined as follows (image of the unit square)
   //  a03--a13--a23--a33
   //   |    |    |    |
   //  a02--a12--a22--a32
-  //   |         |    |
-  //  a01  a11  a21--a31
-  //   |         |    |
+  //   |    |    |    |
+  //  a01--a11--a21--a31
+  //   |    |    |    |
   //  a00--a10--a20--a30
   TinyMatrix<Dimension, NumberOfPoints> m_points_images;
 
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index fd184f7c7..bb0667a02 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -8,9 +8,11 @@
 #include <geometry/PrismTransformation.hpp>
 #include <geometry/PyramidTransformation.hpp>
 #include <geometry/SquareQ2Transformation.hpp>
+#include <geometry/SquareQ3Transformation.hpp>
 #include <geometry/SquareTransformation.hpp>
 #include <geometry/TetrahedronTransformation.hpp>
 #include <geometry/TriangleP2Transformation.hpp>
+#include <geometry/TriangleP3Transformation.hpp>
 #include <geometry/TriangleTransformation.hpp>
 #include <language/utils/PugsFunctionAdapter.hpp>
 #include <mesh/CellType.hpp>
@@ -101,8 +103,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
                   "Polynomial meshes are only supported in dimension 2");
 
     if constexpr (is_polynomial_mesh_v<MeshType>) {
-      if (mesh.degree() != 2) {
-        throw NotImplementedError("IntegrateOnCells only support polynomial meshes of degree 2");
+      if (mesh.degree() > 3) {
+        throw NotImplementedError("IntegrateOnCells only support polynomial meshes of degree <= 3");
       }
     }
 
@@ -155,7 +157,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
       }
     }();
 
-    parallel_for(cell_list.size(), [=, &expression, &tokens, &invalid_cell, &qf,
+    parallel_for(cell_list.size(), [=, &mesh, &connectivity, &expression, &tokens, &invalid_cell, &qf,
                                     &cell_list](typename ListTypeT::index_type index) {
       const int32_t t        = tokens.acquire();
       auto& execution_policy = context_list[t];
@@ -199,34 +201,84 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
             }
           } else {
             static_assert(is_polynomial_mesh_v<MeshType>);
+            if (mesh.degree() == 2) {
+              const auto cell_to_face = cell_to_face_matrix[cell_id];
 
-            const auto cell_to_face = cell_to_face_matrix[cell_id];
+              const Rd& a00 = xr[cell_to_node[0]];
+              const Rd& a01 = xl[cell_to_face[1]][0];
+              const Rd& a02 = xr[cell_to_node[2]];
 
-            const Rd& a00 = xr[cell_to_node[0]];
-            const Rd& a01 = xl[cell_to_face[1]][0];
-            const Rd& a02 = xr[cell_to_node[2]];
+              const Rd& a10 = xl[cell_to_face[2]][0];
+              const Rd& a11 = xj[cell_id];
+              const Rd& a20 = xr[cell_to_node[1]];
 
-            const Rd& a10 = xl[cell_to_face[2]][0];
-            const Rd& a11 = xj[cell_id];
-            const Rd& a20 = xr[cell_to_node[1]];
+              const Rd& a21 = xl[cell_to_face[0]][0];
+              const Rd& a12 = xr[cell_to_node[2]];
+              const Rd& a22 = xr[cell_to_node[2]];
 
-            const Rd& a21 = xl[cell_to_face[0]][0];
-            const Rd& a12 = xr[cell_to_node[2]];
-            const Rd& a22 = xr[cell_to_node[2]];
+              const SquareQ2Transformation<2> T(a00, a01, a02,   //
+                                                a10, a11, a12,   //
+                                                a20, a21, a22);
 
-            const SquareQ2Transformation<2> T(a00, a01, a02,   //
-                                              a10, a11, a12,   //
-                                              a20, a21, a22);
+              for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+                const auto xi = qf.point(i_point);
+                Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+                auto result = expression.execute(execution_policy);
+                value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              }
+            } else if (mesh.degree() == 3) {
+              const auto cell_to_face = cell_to_face_matrix[cell_id];
+              FaceValuePerCell<const bool> cell_face_is_reversed;
 
-            for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
-              const auto xi = qf.point(i_point);
-              Adapter::convertArgs(execution_policy.currentContext(), T(xi));
-              auto result = expression.execute(execution_policy);
-              value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              if constexpr (Dimension > 1) {
+                cell_face_is_reversed = connectivity.cellFaceIsReversed();
+              }
+
+              const Rd& a00 = xr[cell_to_node[0]];
+              const Rd& a30 = xr[cell_to_node[1]];
+              const Rd& a33 = xr[cell_to_node[2]];
+              const Rd& a03 = xr[cell_to_node[2]];
+
+              const Rd& a23 = xr[cell_to_node[2]];
+              const Rd& a13 = xr[cell_to_node[2]];
+
+              const size_t face01_node0 = (cell_face_is_reversed[cell_id][2]) ? 1 : 0;
+
+              const Rd& a10 = xl[cell_to_face[2]][face01_node0];
+              const Rd& a20 = xl[cell_to_face[2]][1 - face01_node0];
+
+              const size_t face12_node0 = (cell_face_is_reversed[cell_id][0]) ? 1 : 0;
+
+              const Rd& a31 = xl[cell_to_face[0]][face12_node0];
+              const Rd& a32 = xl[cell_to_face[0]][1 - face12_node0];
+
+              const size_t face30_node0 = (cell_face_is_reversed[cell_id][1]) ? 1 : 0;
+
+              const Rd& a02 = xl[cell_to_face[1]][face30_node0];
+              const Rd& a01 = xl[cell_to_face[1]][1 - face30_node0];
+
+              const Rd& a11 = 1. / 3 * xr[cell_to_node[0]] + 2. / 3 * xj[cell_id];
+              const Rd& a21 = 1. / 3 * xr[cell_to_node[1]] + 2. / 3 * xj[cell_id];
+
+              const Rd& a12 = 1. / 3 * xr[cell_to_node[3]] + 2. / 3 * xj[cell_id];
+              const Rd& a22 = 1. / 3 * xr[cell_to_node[2]] + 2. / 3 * xj[cell_id];
+
+              const SquareQ3Transformation<2> T(a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31,
+                                                a32, a33);
+
+              for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+                const auto xi = qf.point(i_point);
+                Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+                auto result = expression.execute(execution_policy);
+                value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              }
+
+            } else {
+              throw NotImplementedError("IntegratedOnCells only for mesh_degree<4");
             }
           }
           break;
-        }
+        }   // triangle
         case CellType::Quadrangle: {
           if constexpr (is_polygonal_mesh_v<MeshType>) {
             const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
@@ -243,38 +295,89 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
 
             const auto cell_to_face = cell_to_face_matrix[cell_id];
 
-            const Rd& a00 = xr[cell_to_node[0]];
-            const Rd& a01 = xl[cell_to_face[3]][0];
-            const Rd& a02 = xr[cell_to_node[3]];
+            FaceValuePerCell<const bool> cell_face_is_reversed;
+            if constexpr (Dimension > 1) {
+              cell_face_is_reversed = connectivity.cellFaceIsReversed();
+            }
 
-            const Rd& a10 = xl[cell_to_face[0]][0];
-            const Rd& a11 = xj[cell_id];
-            const Rd& a20 = xr[cell_to_node[1]];
+            if (mesh.degree() == 2) {
+              const Rd& a00 = xr[cell_to_node[0]];
+              const Rd& a01 = xl[cell_to_face[3]][0];
+              const Rd& a02 = xr[cell_to_node[3]];
 
-            const Rd& a21 = xl[cell_to_face[1]][0];
-            const Rd& a12 = xl[cell_to_face[2]][0];
-            const Rd& a22 = xr[cell_to_node[2]];
+              const Rd& a10 = xl[cell_to_face[0]][0];
+              const Rd& a11 = xj[cell_id];
+              const Rd& a20 = xr[cell_to_node[1]];
 
-            const SquareQ2Transformation<2> T(a00, a01, a02,   //
-                                              a10, a11, a12,   //
-                                              a20, a21, a22);
+              const Rd& a21 = xl[cell_to_face[1]][0];
+              const Rd& a12 = xl[cell_to_face[2]][0];
+              const Rd& a22 = xr[cell_to_node[2]];
 
-            for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
-              const auto xi = qf.point(i_point);
-              Adapter::convertArgs(execution_policy.currentContext(), T(xi));
-              auto result = expression.execute(execution_policy);
-              value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              const SquareQ2Transformation<2> T(a00, a01, a02,   //
+                                                a10, a11, a12,   //
+                                                a20, a21, a22);
+
+              for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+                const auto xi = qf.point(i_point);
+                Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+                auto result = expression.execute(execution_policy);
+                value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              }
+            } else if (mesh.degree() == 3) {
+              const Rd& a00 = xr[cell_to_node[0]];
+              const Rd& a30 = xr[cell_to_node[1]];
+              const Rd& a33 = xr[cell_to_node[2]];
+              const Rd& a03 = xr[cell_to_node[3]];
+
+              const size_t face01_node0 = (cell_face_is_reversed[cell_id][0]) ? 1 : 0;
+
+              const Rd& a10 = xl[cell_to_face[0]][face01_node0];
+              const Rd& a20 = xl[cell_to_face[0]][1 - face01_node0];
+
+              const size_t face12_node0 = (cell_face_is_reversed[cell_id][1]) ? 1 : 0;
+
+              const Rd& a31 = xl[cell_to_face[1]][face12_node0];
+              const Rd& a32 = xl[cell_to_face[1]][1 - face12_node0];
+
+              const size_t face23_node0 = (cell_face_is_reversed[cell_id][2]) ? 1 : 0;
+
+              const Rd& a23 = xl[cell_to_face[2]][face23_node0];
+              const Rd& a13 = xl[cell_to_face[2]][1 - face23_node0];
+
+              const size_t face30_node0 = (cell_face_is_reversed[cell_id][3]) ? 1 : 0;
+
+              const Rd& a02 = xl[cell_to_face[3]][face30_node0];
+              const Rd& a01 = xl[cell_to_face[3]][1 - face30_node0];
+
+              const Rd& a11 = 1. / 3 * xr[cell_to_node[0]] + 2. / 3 * xj[cell_id];
+              const Rd& a21 = 1. / 3 * xr[cell_to_node[1]] + 2. / 3 * xj[cell_id];
+
+              const Rd& a12 = 1. / 3 * xr[cell_to_node[3]] + 2. / 3 * xj[cell_id];
+              const Rd& a22 = 1. / 3 * xr[cell_to_node[2]] + 2. / 3 * xj[cell_id];
+
+              const SquareQ3Transformation<2> T(a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31,
+                                                a32, a33);
+
+              for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+                const auto xi = qf.point(i_point);
+                Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+                auto result = expression.execute(execution_policy);
+                value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              }
+
+            } else {
+              throw NotImplementedError("IntegratedOnCells only for mesh_degree<4");
             }
           }
           break;
-        }
-          // LCOV_EXCL_START
+        }   // quad
+            //  LCOV_EXCL_START
         default: {
           invalid_cell = std::make_pair(true, cell_id);
           break;
         }
           // LCOV_EXCL_STOP
-        }
+        }   // switch
       } else {
         static_assert(Dimension == 3);
 
@@ -410,9 +513,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           break;
         }
           // LCOV_EXCL_STOP
-        }
-      }
-
+        }   // switch
+      }     // else
       tokens.release(t);
     });
 
@@ -440,8 +542,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
                   "Polynomial meshes are only supported in dimension 2");
 
     if constexpr (is_polynomial_mesh_v<MeshType>) {
-      if (mesh.degree() != 2) {
-        throw NotImplementedError("IntegrateOnCells only support polynomial meshes of degree 2");
+      if (mesh.degree() > 3) {
+        throw NotImplementedError("IntegrateOnCells only support polynomial meshes of degree <= 3");
       }
     }
 
@@ -473,7 +575,12 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
 
     const auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
     const auto cell_type           = connectivity.cellType();
-    const auto xr                  = mesh.xr();
+
+    FaceValuePerCell<const bool> cell_face_is_reversed;
+    if constexpr (Dimension > 1) {
+      cell_face_is_reversed = connectivity.cellFaceIsReversed();
+    }
+    const auto xr = mesh.xr();
 
     const auto cell_to_face_matrix = connectivity.cellToFaceMatrix();
     FaceArray<const Rd> xl;
@@ -485,7 +592,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
 
     auto invalid_cell = std::make_pair(false, CellId{});
 
-    parallel_for(cell_list.size(), [=, &expression, &tokens, &invalid_cell, &quadrature_descriptor,
+    parallel_for(cell_list.size(), [=, &mesh, &expression, &tokens, &invalid_cell, &quadrature_descriptor,
                                     &cell_list](const typename ListTypeT::index_type& index) {
       const int32_t t        = tokens.acquire();
       auto& execution_policy = context_list[t];
@@ -532,24 +639,58 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
             static_assert(is_polynomial_mesh_v<MeshType>);
 
             const auto cell_to_face = cell_to_face_matrix[cell_id];
+            const auto qf           = QuadratureManager::instance().getTriangleFormula(quadrature_descriptor);
 
-            const Rd& a0 = xr[cell_to_node[0]];
-            const Rd& a1 = xr[cell_to_node[1]];
-            const Rd& a2 = xr[cell_to_node[2]];
+            if (mesh.degree() == 2) {
+              const Rd& a0 = xr[cell_to_node[0]];
+              const Rd& a1 = xr[cell_to_node[1]];
+              const Rd& a2 = xr[cell_to_node[2]];
 
-            const Rd& a01 = xl[cell_to_face[2]][0];
-            const Rd& a12 = xl[cell_to_face[0]][0];
-            const Rd& a20 = xl[cell_to_face[1]][0];
+              const Rd& a01 = xl[cell_to_face[2]][0];
+              const Rd& a12 = xl[cell_to_face[0]][0];
+              const Rd& a20 = xl[cell_to_face[1]][0];
 
-            const TriangleP2Transformation<2> T(a0, a1, a2,   //
-                                                a12, a20, a01);
-            const auto qf = QuadratureManager::instance().getTriangleFormula(quadrature_descriptor);
+              const TriangleP2Transformation<2> T(a0, a1, a2,   //
+                                                  a12, a20, a01);
 
-            for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
-              const auto xi = qf.point(i_point);
-              Adapter::convertArgs(execution_policy.currentContext(), T(xi));
-              auto result = expression.execute(execution_policy);
-              value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+                const auto xi = qf.point(i_point);
+                Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+                auto result = expression.execute(execution_policy);
+                value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              }
+            } else if (mesh.degree() == 3) {
+              const Rd& a0 = xr[cell_to_node[0]];
+              const Rd& a1 = xr[cell_to_node[1]];
+              const Rd& a2 = xr[cell_to_node[2]];
+
+              const size_t face01_node0 = (cell_face_is_reversed[cell_id][2]) ? 1 : 0;
+
+              const Rd& a001 = xl[cell_to_face[2]][face01_node0];
+              const Rd& a011 = xl[cell_to_face[2]][1 - face01_node0];
+
+              const size_t face12_node0 = (cell_face_is_reversed[cell_id][0]) ? 1 : 0;
+
+              const Rd& a112 = xl[cell_to_face[0]][face12_node0];
+              const Rd& a122 = xl[cell_to_face[0]][1 - face12_node0];
+
+              const size_t face20_node0 = (cell_face_is_reversed[cell_id][1]) ? 1 : 0;
+
+              const Rd& a220 = xl[cell_to_face[1]][face20_node0];
+              const Rd& a200 = xl[cell_to_face[1]][1 - face20_node0];
+
+              const Rd& a012 = 1. / 3 * (a0 + a1 + a2);
+
+              const TriangleP3Transformation<2> T(a0, a1, a2, a112, a122, a220, a200, a001, a011, a012);
+
+              for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+                const auto xi = qf.point(i_point);
+                Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+                auto result = expression.execute(execution_policy);
+                value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              }
+            } else {
+              throw NotImplementedError("IntegratedOnCells only for mesh_degree<4");
             }
           }
           break;
@@ -570,40 +711,87 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
             static_assert(is_polynomial_mesh_v<MeshType>);
 
             const auto cell_to_face = cell_to_face_matrix[cell_id];
+            if (mesh.degree() == 2) {
+              const Rd& a00 = xr[cell_to_node[0]];
+              const Rd& a01 = xl[cell_to_face[3]][0];
+              const Rd& a02 = xr[cell_to_node[3]];
 
-            const Rd& a00 = xr[cell_to_node[0]];
-            const Rd& a01 = xl[cell_to_face[3]][0];
-            const Rd& a02 = xr[cell_to_node[3]];
+              const Rd& a10 = xl[cell_to_face[0]][0];
+              const Rd& a11 = xj[cell_id];
+              const Rd& a12 = xl[cell_to_face[2]][0];
 
-            const Rd& a10 = xl[cell_to_face[0]][0];
-            const Rd& a11 = xj[cell_id];
-            const Rd& a20 = xr[cell_to_node[1]];
+              const Rd& a20 = xr[cell_to_node[1]];
+              const Rd& a21 = xl[cell_to_face[1]][0];
+              const Rd& a22 = xr[cell_to_node[2]];
 
-            const Rd& a21 = xl[cell_to_face[1]][0];
-            const Rd& a12 = xl[cell_to_face[2]][0];
-            const Rd& a22 = xr[cell_to_node[2]];
+              const SquareQ2Transformation<2> T(a00, a01, a02,   //
+                                                a10, a11, a12,   //
+                                                a20, a21, a22);
+              const auto qf = QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
 
-            const SquareQ2Transformation<2> T(a00, a01, a02,   //
-                                              a10, a11, a12,   //
-                                              a20, a21, a22);
-            const auto qf = QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
+              for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+                const auto xi = qf.point(i_point);
+                Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+                auto result = expression.execute(execution_policy);
+                value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              }
 
-            for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
-              const auto xi = qf.point(i_point);
-              Adapter::convertArgs(execution_policy.currentContext(), T(xi));
-              auto result = expression.execute(execution_policy);
-              value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+            } else if (mesh.degree() == 3) {
+              const Rd& a00 = xr[cell_to_node[0]];
+              const Rd& a30 = xr[cell_to_node[1]];
+              const Rd& a33 = xr[cell_to_node[2]];
+              const Rd& a03 = xr[cell_to_node[3]];
+
+              const size_t face01_node0 = (cell_face_is_reversed[cell_id][0]) ? 1 : 0;
+
+              const Rd& a10 = xl[cell_to_face[0]][face01_node0];
+              const Rd& a20 = xl[cell_to_face[0]][1 - face01_node0];
+
+              const size_t face12_node0 = (cell_face_is_reversed[cell_id][1]) ? 1 : 0;
+
+              const Rd& a31 = xl[cell_to_face[1]][face12_node0];
+              const Rd& a32 = xl[cell_to_face[1]][1 - face12_node0];
+
+              const size_t face23_node0 = (cell_face_is_reversed[cell_id][2]) ? 1 : 0;
+
+              const Rd& a23 = xl[cell_to_face[2]][face23_node0];
+              const Rd& a13 = xl[cell_to_face[2]][1 - face23_node0];
+
+              const size_t face30_node0 = (cell_face_is_reversed[cell_id][3]) ? 1 : 0;
+
+              const Rd& a02 = xl[cell_to_face[3]][face30_node0];
+              const Rd& a01 = xl[cell_to_face[3]][1 - face30_node0];
+
+              const Rd& a11 = 1. / 3 * xr[cell_to_node[0]] + 2. / 3 * xj[cell_id];
+              const Rd& a21 = 1. / 3 * xr[cell_to_node[1]] + 2. / 3 * xj[cell_id];
+
+              const Rd& a12 = 1. / 3 * xr[cell_to_node[3]] + 2. / 3 * xj[cell_id];
+              const Rd& a22 = 1. / 3 * xr[cell_to_node[2]] + 2. / 3 * xj[cell_id];
+
+              const SquareQ3Transformation<2> T(a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31,
+                                                a32, a33);
+              const auto qf = QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
+
+              for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+                const auto xi = qf.point(i_point);
+                Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+                auto result = expression.execute(execution_policy);
+                value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+              }
+            } else {
+              throw NotImplementedError("IntegratedOnCells only for mesh_degree<4");
             }
           }
           break;
-        }
-          // LCOV_EXCL_START
+        }   // case quad
+            //  LCOV_EXCL_START
         default: {
           invalid_cell = std::make_pair(true, cell_id);
           break;
         }
           // LCOV_EXCL_STOP
-        }
+
+        }   // switch
       } else {
         static_assert(Dimension == 3);
 
@@ -730,9 +918,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           break;
         }
           // LCOV_EXCL_STOP
-        }
-      }
-
+        }   // switch
+      }     // fin test dimension
       tokens.release(t);
     });
 
diff --git a/src/mesh/MeshData.cpp b/src/mesh/MeshData.cpp
index 286b960c0..0eadfdf30 100644
--- a/src/mesh/MeshData.cpp
+++ b/src/mesh/MeshData.cpp
@@ -1,5 +1,8 @@
 #include <mesh/MeshData.hpp>
 
+#include <geometry/LineCubicTransformation.hpp>
+#include <geometry/LineParabolicTransformation.hpp>
+#include <geometry/LineTransformation.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemValueVariant.hpp>
 #include <mesh/Mesh.hpp>
@@ -32,3 +35,196 @@ MeshData<PolynomialMesh<2>>::_storeBadMesh()
 {
   throw NotImplementedError("storeBadMesh");
 }
+
+template <size_t MeshDegree>
+auto getTransformation(const NodeValue<const TinyVector<2>>& xr,
+                       const FaceArray<const TinyVector<2>>& xl,
+                       const ItemToItemMatrix<ItemType::face, ItemType::node>& face_to_node_matrix,
+                       const FaceId face_id);
+
+template <>
+auto
+getTransformation<1>(const NodeValue<const TinyVector<2>>& xr,
+                     const FaceArray<const TinyVector<2>>&,
+                     const ItemToItemMatrix<ItemType::face, ItemType::node>& face_to_node_matrix,
+                     const FaceId face_id)
+{
+  return LineTransformation<2>{xr[face_to_node_matrix[face_id][0]],   //
+                               xr[face_to_node_matrix[face_id][1]]};
+}
+
+template <>
+auto
+getTransformation<2>(const NodeValue<const TinyVector<2>>& xr,
+                     const FaceArray<const TinyVector<2>>& xl,
+                     const ItemToItemMatrix<ItemType::face, ItemType::node>& face_to_node_matrix,
+                     const FaceId face_id)
+{
+  return LineParabolicTransformation<2>{xr[face_to_node_matrix[face_id][0]],   //
+                                        xl[face_id][0],                        //
+                                        xr[face_to_node_matrix[face_id][1]]};
+}
+
+template <>
+auto
+getTransformation<3>(const NodeValue<const TinyVector<2>>& xr,
+                     const FaceArray<const TinyVector<2>>& xl,
+                     const ItemToItemMatrix<ItemType::face, ItemType::node>& face_to_node_matrix,
+                     const FaceId face_id)
+{
+  return LineCubicTransformation<2>{xr[face_to_node_matrix[face_id][0]],   //
+                                    xl[face_id][0],                        //
+                                    xl[face_id][1],                        //
+                                    xr[face_to_node_matrix[face_id][1]]};
+}
+
+template <size_t MeshDegree>
+void
+MeshData<PolynomialMesh<2>>::_computeCellCentroidForMeshDegree()
+{
+  FaceValue<Rd> face_integral{m_mesh.connectivity()};
+
+  using R1 = TinyVector<1>;
+
+  const auto xr = m_mesh.xr();
+  const auto xl = m_mesh.xl();
+
+  auto face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+
+  constexpr size_t total_degree = 3 * MeshDegree - 1;
+
+  QuadratureFormula<1> qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(total_degree));
+  parallel_for(
+    m_mesh.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+      auto T = getTransformation<MeshDegree>(xr, xl, face_to_node_matrix, face_id);
+
+      Rd integral = zero;
+      for (size_t i_ksi = 0; i_ksi < qf.numberOfPoints(); ++i_ksi) {
+        const R1& ksi    = qf.point(i_ksi);
+        const Rd x       = T(ksi);
+        const Rd x_prime = T.velocity(ksi);
+        const Rd x_prime_perp{x_prime[1], -x_prime[0]};
+
+        integral += qf.weight(i_ksi) * Rd{x[0] * x[0] * (x_prime_perp[0]),   //
+                                          x[1] * x[1] * (x_prime_perp[1])};
+      }
+
+      face_integral[face_id] = integral;
+    });
+
+  auto cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
+  auto cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
+
+  CellValue<const double> Vj = this->Vj();
+  CellValue<Rd> cell_centroid{m_mesh.connectivity()};
+  parallel_for(
+    m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+      auto cell_faces = cell_to_face_matrix[cell_id];
+      Rd xj           = zero;
+      for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
+        const FaceId face_id     = cell_faces[i_face];
+        const double orientation = (cell_face_is_reversed[cell_id][i_face]) ? -1 : 1;
+
+        xj += orientation * face_integral[face_id];
+      }
+      cell_centroid[cell_id] = 1 / (2 * Vj[cell_id]) * xj;
+    });
+
+  m_cell_centroid = cell_centroid;
+}
+
+void
+MeshData<PolynomialMesh<2>>::_computeCellCentroid()
+{
+  switch (m_mesh.degree()) {
+  case 1: {
+    this->_computeCellCentroidForMeshDegree<1>();
+    break;
+  }
+  case 2: {
+    this->_computeCellCentroidForMeshDegree<2>();
+    break;
+  }
+  case 3: {
+    this->_computeCellCentroidForMeshDegree<3>();
+    break;
+  }
+  default: {
+    throw NotImplementedError("mesh degree > 3");
+  }
+  }
+}
+template <size_t MeshDegree>
+void
+MeshData<PolynomialMesh<2>>::_computeCellVolumeForMeshDegree()
+{
+  FaceValue<double> face_integral{m_mesh.connectivity()};
+
+  using R1 = TinyVector<1>;
+
+  const auto xr = m_mesh.xr();
+  const auto xl = m_mesh.xl();
+
+  auto face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+
+  constexpr size_t total_degree = 2 * MeshDegree - 1;
+
+  QuadratureFormula<1> qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(total_degree));
+  parallel_for(
+    m_mesh.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+      auto T = getTransformation<MeshDegree>(xr, xl, face_to_node_matrix, face_id);
+
+      double integral = 0;
+      for (size_t i_ksi = 0; i_ksi < qf.numberOfPoints(); ++i_ksi) {
+        const R1& ksi    = qf.point(i_ksi);
+        const Rd x       = T(ksi);
+        const Rd x_prime = T.velocity(ksi);
+        const Rd x_prime_perp{x_prime[1], -x_prime[0]};
+
+        integral += qf.weight(i_ksi) * dot(x, x_prime_perp);
+      }
+
+      face_integral[face_id] = inv_Dimension * integral;
+    });
+
+  auto cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
+  auto cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
+
+  CellValue<double> Vj{m_mesh.connectivity()};
+  parallel_for(
+    m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+      auto cell_faces = cell_to_face_matrix[cell_id];
+      double volume   = 0;
+      for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
+        const FaceId face_id     = cell_faces[i_face];
+        const double orientation = (cell_face_is_reversed[cell_id][i_face]) ? -1 : 1;
+
+        volume += orientation * face_integral[face_id];
+      }
+      Vj[cell_id] = volume;
+    });
+
+  m_Vj = Vj;
+}
+
+void
+MeshData<PolynomialMesh<2>>::_computeCellVolume()
+{
+  switch (m_mesh.degree()) {
+  case 1: {
+    this->_computeCellVolumeForMeshDegree<1>();
+    break;
+  }
+  case 2: {
+    this->_computeCellVolumeForMeshDegree<2>();
+    break;
+  }
+  case 3: {
+    this->_computeCellVolumeForMeshDegree<3>();
+    break;
+  }
+  default: {
+    throw NotImplementedError("mesh degree > 3");
+  }
+  }
+}
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index d2ff74768..7f8ae995d 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -721,249 +721,13 @@ class MeshData<PolynomialMesh<2>>
 
   void _storeBadMesh();
 
-  PUGS_INLINE
-  void
-  _computeCellCentroid()
-  {
-    FaceValue<Rd> face_integral{m_mesh.connectivity()};
-
-    if (m_mesh.degree() == 1) {
-      using R1 = TinyVector<1>;
-
-      auto w0 = [](const R1 x) { return 0.5 * (1 - x[0]); };
-      auto w1 = [](const R1 x) { return 0.5 * (1 + x[0]); };
-
-      auto w0_prime = [](const R1) { return -0.5; };
-      auto w1_prime = [](const R1) { return 0.5; };
-
-      const auto xr = m_mesh.xr();
-
-      auto face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
-
-      const size_t total_degree = 2;
-
-      QuadratureFormula<1> qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(total_degree));
-
-      parallel_for(
-        m_mesh.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
-          const Rd a0 = xr[face_to_node_matrix[face_id][0]];
-          const Rd a1 = xr[face_to_node_matrix[face_id][1]];
-
-          Rd integral = zero;
-          for (size_t i_ksi = 0; i_ksi < qf.numberOfPoints(); ++i_ksi) {
-            const R1& ksi    = qf.point(i_ksi);
-            const Rd x       = w0(ksi) * a0 + w1(ksi) * a1;
-            const Rd x_prime = w0_prime(ksi) * a0 + w1_prime(ksi) * a1;
-            const Rd x_prime_perp{x_prime[1], -x_prime[0]};
-
-            integral += qf.weight(i_ksi) * Rd{x[0] * x[0] * (x_prime_perp[0]),   //
-                                              x[1] * x[1] * (x_prime_perp[1])};
-          }
-
-          face_integral[face_id] = integral;
-        });
-
-      auto cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
-      auto cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
-
-      CellValue<const double> Vj = this->Vj();
-      CellValue<Rd> cell_centroid{m_mesh.connectivity()};
-      parallel_for(
-        m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-          auto cell_faces = cell_to_face_matrix[cell_id];
-          Rd xj           = zero;
-          for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
-            const FaceId face_id     = cell_faces[i_face];
-            const double orientation = (cell_face_is_reversed[cell_id][i_face]) ? -1 : 1;
-
-            xj += orientation * face_integral[face_id];
-          }
-          cell_centroid[cell_id] = 1 / (2 * Vj[cell_id]) * xj;
-        });
-
-      m_cell_centroid = cell_centroid;
-    } else if (m_mesh.degree() == 2) {
-      using R1 = TinyVector<1>;
-
-      auto w0 = [](const R1 x) { return 0.5 * (x[0] - 1) * x[0]; };
-      auto w1 = [](const R1 x) { return (1 - x[0]) * (1 + x[0]); };
-      auto w2 = [](const R1 x) { return 0.5 * (x[0] + 1) * x[0]; };
-
-      auto w0_prime = [](const R1 x) { return x[0] - 0.5; };
-      auto w1_prime = [](const R1 x) { return -2 * x[0]; };
-      auto w2_prime = [](const R1 x) { return x[0] + 0.5; };
-
-      const auto xr = m_mesh.xr();
-      const auto xl = m_mesh.xl();
-
-      auto face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
-
-      const size_t total_degree = 3 * m_mesh.degree() - 1;
-
-      QuadratureFormula<1> qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(total_degree));
-
-      parallel_for(
-        m_mesh.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
-          const Rd a0 = xr[face_to_node_matrix[face_id][0]];
-          const Rd a1 = xl[face_id][0];
-          const Rd a2 = xr[face_to_node_matrix[face_id][1]];
-
-          Rd integral = zero;
-          for (size_t i_ksi = 0; i_ksi < qf.numberOfPoints(); ++i_ksi) {
-            const R1& ksi    = qf.point(i_ksi);
-            const Rd x       = w0(ksi) * a0 + w1(ksi) * a1 + w2(ksi) * a2;
-            const Rd x_prime = w0_prime(ksi) * a0 + w1_prime(ksi) * a1 + w2_prime(ksi) * a2;
-            const Rd x_prime_perp{x_prime[1], -x_prime[0]};
-
-            integral += qf.weight(i_ksi) * Rd{x[0] * x[0] * (x_prime_perp[0]),   //
-                                              x[1] * x[1] * (x_prime_perp[1])};
-          }
-
-          face_integral[face_id] = integral;
-        });
-
-      auto cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
-      auto cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
-
-      CellValue<const double> Vj = this->Vj();
-      CellValue<Rd> cell_centroid{m_mesh.connectivity()};
-      parallel_for(
-        m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-          auto cell_faces = cell_to_face_matrix[cell_id];
-          Rd xj           = zero;
-          for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
-            const FaceId face_id     = cell_faces[i_face];
-            const double orientation = (cell_face_is_reversed[cell_id][i_face]) ? -1 : 1;
-
-            xj += orientation * face_integral[face_id];
-          }
-          cell_centroid[cell_id] = 1 / (2 * Vj[cell_id]) * xj;
-        });
-
-      m_cell_centroid = cell_centroid;
-    } else {
-      throw NotImplementedError("only available for degree <= 2");
-    }
-  }
-
-  PUGS_INLINE
-  void
-  _computeCellVolume()
-  {
-    FaceValue<double> face_integral{m_mesh.connectivity()};
-
-    if (m_mesh.degree() == 1) {
-      using R1 = TinyVector<1>;
-
-      auto w0 = [](const R1 x) { return 0.5 * (1 - x[0]); };
-      auto w1 = [](const R1 x) { return 0.5 * (1 + x[0]); };
-
-      auto w0_prime = [](const R1) { return -0.5; };
-      auto w1_prime = [](const R1) { return 0.5; };
-
-      const auto xr = m_mesh.xr();
+  template <size_t MeshDegree>
+  void _computeCellCentroidForMeshDegree();
+  void _computeCellCentroid();
 
-      auto face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
-
-      const size_t total_degree = 2 * m_mesh.degree() - 1;
-      QuadratureFormula<1> qf   = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(total_degree));
-
-      parallel_for(
-        m_mesh.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
-          const Rd a0 = xr[face_to_node_matrix[face_id][0]];
-          const Rd a1 = xr[face_to_node_matrix[face_id][1]];
-
-          double integral = 0;
-          for (size_t i_ksi = 0; i_ksi < qf.numberOfPoints(); ++i_ksi) {
-            const R1& ksi    = qf.point(i_ksi);
-            const Rd x       = w0(ksi) * a0 + w1(ksi) * a1;
-            const Rd x_prime = w0_prime(ksi) * a0 + w1_prime(ksi) * a1;
-            const Rd x_prime_perp{x_prime[1], -x_prime[0]};
-
-            integral += qf.weight(i_ksi) * dot(x, x_prime_perp);
-          }
-
-          face_integral[face_id] = inv_Dimension * integral;
-        });
-
-      auto cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
-      auto cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
-
-      CellValue<double> Vj{m_mesh.connectivity()};
-      parallel_for(
-        m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-          auto cell_faces = cell_to_face_matrix[cell_id];
-          double volume   = 0;
-          for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
-            const FaceId face_id     = cell_faces[i_face];
-            const double orientation = (cell_face_is_reversed[cell_id][i_face]) ? -1 : 1;
-
-            volume += orientation * face_integral[face_id];
-          }
-          Vj[cell_id] = volume;
-        });
-
-      m_Vj = Vj;
-    } else if (m_mesh.degree() == 2) {
-      using R1 = TinyVector<1>;
-
-      auto w0 = [](const R1 x) { return 0.5 * (x[0] - 1) * x[0]; };
-      auto w1 = [](const R1 x) { return (1 - x[0]) * (1 + x[0]); };
-      auto w2 = [](const R1 x) { return 0.5 * (x[0] + 1) * x[0]; };
-
-      auto w0_prime = [](const R1 x) { return x[0] - 0.5; };
-      auto w1_prime = [](const R1 x) { return -2 * x[0]; };
-      auto w2_prime = [](const R1 x) { return x[0] + 0.5; };
-
-      const auto xr = m_mesh.xr();
-      const auto xl = m_mesh.xl();
-
-      auto face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
-
-      const size_t total_degree = 2 * m_mesh.degree() - 1;
-      QuadratureFormula<1> qf   = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(total_degree));
-
-      parallel_for(
-        m_mesh.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
-          const Rd a0 = xr[face_to_node_matrix[face_id][0]];
-          const Rd a1 = xl[face_id][0];
-          const Rd a2 = xr[face_to_node_matrix[face_id][1]];
-
-          double integral = 0;
-          for (size_t i_ksi = 0; i_ksi < qf.numberOfPoints(); ++i_ksi) {
-            const R1& ksi    = qf.point(i_ksi);
-            const Rd x       = w0(ksi) * a0 + w1(ksi) * a1 + w2(ksi) * a2;
-            const Rd x_prime = w0_prime(ksi) * a0 + w1_prime(ksi) * a1 + w2_prime(ksi) * a2;
-            const Rd x_prime_perp{x_prime[1], -x_prime[0]};
-
-            integral += qf.weight(i_ksi) * dot(x, x_prime_perp);
-          }
-
-          face_integral[face_id] = inv_Dimension * integral;
-        });
-
-      auto cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
-      auto cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
-
-      CellValue<double> Vj{m_mesh.connectivity()};
-      parallel_for(
-        m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-          auto cell_faces = cell_to_face_matrix[cell_id];
-          double volume   = 0;
-          for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
-            const FaceId face_id     = cell_faces[i_face];
-            const double orientation = (cell_face_is_reversed[cell_id][i_face]) ? -1 : 1;
-
-            volume += orientation * face_integral[face_id];
-          }
-          Vj[cell_id] = volume;
-        });
-
-      m_Vj = Vj;
-    } else {
-      throw NotImplementedError("only available for degree <= 2");
-    }
-  }
+  template <size_t MeshDegree>
+  void _computeCellVolumeForMeshDegree();
+  void _computeCellVolume();
 
   void
   _checkCellVolume()
diff --git a/src/scheme/VariationalSolver.cpp b/src/scheme/VariationalSolver.cpp
index 5b7f1245d..043cc4bbd 100644
--- a/src/scheme/VariationalSolver.cpp
+++ b/src/scheme/VariationalSolver.cpp
@@ -215,8 +215,10 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
             const Rd x1 = [&]() {
               if (mesh.degree() == 2) {
                 return mesh.xl()[face_id][0];
-              } else {
+              } else if (mesh.degree() == 1) {
                 return 0.5 * (x0 + x2);
+              } else {
+                throw NotImplementedError("degree>2");
               }
             }();
 
-- 
GitLab