diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 0b5f7e30b8675e4745040a22878ae26268d9880f..58070b1b39d0641a4d3dffbd6dcda5a75dca3c7e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -732,14 +732,8 @@ SchemeModule::SchemeModule()
                                 -> std::shared_ptr<const DiscreteFunctionVariant> {
                                 return std::visit(
                                   [&](auto&& mesh) -> std::shared_ptr<const DiscreteFunctionVariant> {
-                                    using MeshType = mesh_type_t<decltype(mesh)>;
-                                    if constexpr (is_polygonal_mesh_v<MeshType>) {
-                                      return std::make_shared<DiscreteFunctionVariant>(
-                                        DiscreteFunctionP0(mesh_v,
-                                                           MeshDataManager::instance().getMeshData(*mesh).Vj()));
-                                    } else {
-                                      throw NotImplementedError("unexpected mesh type");
-                                    }
+                                    return std::make_shared<DiscreteFunctionVariant>(
+                                      DiscreteFunctionP0(mesh_v, MeshDataManager::instance().getMeshData(*mesh).Vj()));
                                   },
                                   mesh_v->variant());
                               }
diff --git a/src/mesh/MeshCellZone.cpp b/src/mesh/MeshCellZone.cpp
index 08bdf64ba40c72449cabf2670a9c8c1112d0f2e8..4807bccbcdfe3798f8cb797919333ba4ccc54c84 100644
--- a/src/mesh/MeshCellZone.cpp
+++ b/src/mesh/MeshCellZone.cpp
@@ -2,23 +2,24 @@
 
 #include <mesh/Connectivity.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/PolynomialMesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <size_t Dimension>
-MeshCellZone<Dimension>::MeshCellZone(const Mesh<Dimension>&, const RefCellList& ref_cell_list)
+template <MeshConcept MeshType>
+MeshCellZone::MeshCellZone(const MeshType&, const RefCellList& ref_cell_list)
   : m_cell_list(ref_cell_list.list()), m_zone_name(ref_cell_list.refId().tagName())
 {}
 
-template <size_t Dimension>
-MeshCellZone<Dimension>
-getMeshCellZone(const Mesh<Dimension>& mesh, const IZoneDescriptor& zone_descriptor)
+template <MeshConcept MeshType>
+MeshCellZone
+getMeshCellZone(const MeshType& mesh, const IZoneDescriptor& zone_descriptor)
 {
   for (size_t i_ref_cell_list = 0; i_ref_cell_list < mesh.connectivity().template numberOfRefItemList<ItemType::cell>();
        ++i_ref_cell_list) {
     const auto& ref_cell_list = mesh.connectivity().template refItemList<ItemType::cell>(i_ref_cell_list);
     const RefId& ref          = ref_cell_list.refId();
     if (ref == zone_descriptor) {
-      return MeshCellZone<Dimension>{mesh, ref_cell_list};
+      return MeshCellZone{mesh, ref_cell_list};
     }
   }
 
@@ -28,6 +29,8 @@ getMeshCellZone(const Mesh<Dimension>& mesh, const IZoneDescriptor& zone_descrip
   throw NormalError(ost.str());
 }
 
-template MeshCellZone<1> getMeshCellZone(const Mesh<1>&, const IZoneDescriptor&);
-template MeshCellZone<2> getMeshCellZone(const Mesh<2>&, const IZoneDescriptor&);
-template MeshCellZone<3> getMeshCellZone(const Mesh<3>&, const IZoneDescriptor&);
+template MeshCellZone getMeshCellZone(const Mesh<1>&, const IZoneDescriptor&);
+template MeshCellZone getMeshCellZone(const Mesh<2>&, const IZoneDescriptor&);
+template MeshCellZone getMeshCellZone(const Mesh<3>&, const IZoneDescriptor&);
+
+template MeshCellZone getMeshCellZone(const PolynomialMesh<2>&, const IZoneDescriptor&);
diff --git a/src/mesh/MeshCellZone.hpp b/src/mesh/MeshCellZone.hpp
index 1f8827b1b2998678f9f73ddbbdd2c18a18ed753d..bbf1927fc925dca6f45816e33d48553e94045bf2 100644
--- a/src/mesh/MeshCellZone.hpp
+++ b/src/mesh/MeshCellZone.hpp
@@ -2,13 +2,13 @@
 #define MESH_CELL_ZONE_HPP
 
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
 template <size_t Dimension>
 class Mesh;
 
-template <size_t Dimension>
 class [[nodiscard]] MeshCellZone   // clazy:exclude=copyable-polymorphic
 {
  protected:
@@ -16,9 +16,8 @@ class [[nodiscard]] MeshCellZone   // clazy:exclude=copyable-polymorphic
   std::string m_zone_name;
 
  public:
-  template <size_t MeshDimension>
-  friend MeshCellZone<MeshDimension> getMeshCellZone(const Mesh<MeshDimension>& mesh,
-                                                     const IZoneDescriptor& zone_descriptor);
+  template <MeshConcept MeshTypeT>
+  friend MeshCellZone getMeshCellZone(const MeshTypeT& mesh, const IZoneDescriptor& zone_descriptor);
 
   MeshCellZone& operator=(const MeshCellZone&) = default;
   MeshCellZone& operator=(MeshCellZone&&)      = default;
@@ -30,7 +29,8 @@ class [[nodiscard]] MeshCellZone   // clazy:exclude=copyable-polymorphic
   }
 
  protected:
-  MeshCellZone(const Mesh<Dimension>& mesh, const RefCellList& ref_cell_list);
+  template <MeshConcept MeshType>
+  MeshCellZone(const MeshType& mesh, const RefCellList& ref_cell_list);
 
  public:
   MeshCellZone(const MeshCellZone&) = default;
@@ -40,7 +40,7 @@ class [[nodiscard]] MeshCellZone   // clazy:exclude=copyable-polymorphic
   virtual ~MeshCellZone() = default;
 };
 
-template <size_t Dimension>
-MeshCellZone<Dimension> getMeshCellZone(const Mesh<Dimension>& mesh, const IZoneDescriptor& zone_descriptor);
+template <MeshConcept MeshType>
+MeshCellZone getMeshCellZone(const MeshType& mesh, const IZoneDescriptor& zone_descriptor);
 
 #endif   // MESH_CELL_ZONE_HPP
diff --git a/src/mesh/MeshData.cpp b/src/mesh/MeshData.cpp
index 40e3878624b1ae84d1923063df6bc6b707926104..286b960c0b3ac52db55ec244259d6f7bc5a074c9 100644
--- a/src/mesh/MeshData.cpp
+++ b/src/mesh/MeshData.cpp
@@ -4,6 +4,7 @@
 #include <mesh/ItemValueVariant.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshVariant.hpp>
+#include <mesh/PolynomialMesh.hpp>
 #include <output/NamedItemValueVariant.hpp>
 #include <output/VTKWriter.hpp>
 #include <utils/Exceptions.hpp>
@@ -25,3 +26,9 @@ MeshData<Mesh<Dimension>>::_storeBadMesh()
 template void MeshData<Mesh<1>>::_storeBadMesh();
 template void MeshData<Mesh<2>>::_storeBadMesh();
 template void MeshData<Mesh<3>>::_storeBadMesh();
+
+void
+MeshData<PolynomialMesh<2>>::_storeBadMesh()
+{
+  throw NotImplementedError("storeBadMesh");
+}
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 4fb87bc5ad4aee87c7dd92e9232ec4508145a14f..1684abb2f63093197181f8a1d6cf18fe6750c0fd 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -2,15 +2,16 @@
 #define MESH_DATA_HPP
 
 #include <algebra/TinyVector.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
 #include <mesh/ItemValue.hpp>
+#include <mesh/Mesh.hpp>
 #include <mesh/MeshTraits.hpp>
+#include <mesh/PolynomialMesh.hpp>
 #include <mesh/SubItemValuePerItem.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PugsUtils.hpp>
 
-template <size_t Dimension>
-class Mesh;
-
 template <MeshConcept MeshType>
 class MeshData;
 
@@ -701,4 +702,222 @@ class MeshData<Mesh<Dimension>>
   ~MeshData() {}
 };
 
+template <>
+class MeshData<PolynomialMesh<2>>
+{
+ public:
+  static constexpr size_t Dimension = 2;
+
+  using MeshType = PolynomialMesh<Dimension>;
+
+  using Rd = TinyVector<Dimension>;
+
+  static constexpr double inv_Dimension = 1. / Dimension;
+
+ private:
+  const MeshType& m_mesh;
+  CellValue<const Rd> m_cell_centroid;
+  CellValue<const double> m_Vj;
+
+  void _storeBadMesh();
+
+  PUGS_INLINE
+  void
+  _computeCellCentroid()
+  {
+    FaceValue<Rd> face_integral{m_mesh.connectivity()};
+
+    if (m_mesh.degree() != 2) {
+      throw NotImplementedError("only available for 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;
+  }
+
+  PUGS_INLINE
+  void
+  _computeCellVolume()
+  {
+    FaceValue<double> face_integral{m_mesh.connectivity()};
+
+    if (m_mesh.degree() != 2) {
+      throw NotImplementedError("only available for 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;
+  }
+
+  void
+  _checkCellVolume()
+  {
+    Assert(m_Vj.isBuilt());
+
+    bool is_valid = [&] {
+      for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
+        if (m_Vj[j] <= 0) {
+          return false;
+        }
+      }
+      return true;
+    }();
+
+    if (not parallel::allReduceAnd(is_valid)) {
+      this->_storeBadMesh();
+    }
+  }
+
+ public:
+  PUGS_INLINE
+  const MeshType&
+  mesh() const
+  {
+    return m_mesh;
+  }
+
+  PUGS_INLINE
+  CellValue<const Rd>
+  xj()
+  {
+    if (not m_cell_centroid.isBuilt()) {
+      this->_computeCellCentroid();
+    }
+    return m_cell_centroid;
+  }
+
+  PUGS_INLINE
+  CellValue<const double>
+  Vj()
+  {
+    if (not m_Vj.isBuilt()) {
+      this->_computeCellVolume();
+      this->_checkCellVolume();
+    }
+    return m_Vj;
+  }
+
+ private:
+  // MeshData **must** be constructed through MeshDataManager
+  friend class MeshDataManager;
+  MeshData(const MeshType& mesh) : m_mesh(mesh) {}
+
+ public:
+  MeshData() = delete;
+
+  MeshData(const MeshData&) = delete;
+  MeshData(MeshData&&)      = delete;
+
+  ~MeshData() {}
+};
+
 #endif   // MESH_DATA_HPP
diff --git a/src/mesh/MeshDataManager.cpp b/src/mesh/MeshDataManager.cpp
index e95a717d8f90b9aaba86f9c3702f15a06a0df54a..0e53860748244de6b4e0ecba2b00eb21b76ee86b 100644
--- a/src/mesh/MeshDataManager.cpp
+++ b/src/mesh/MeshDataManager.cpp
@@ -5,6 +5,7 @@
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshDataVariant.hpp>
+#include <mesh/PolynomialMesh.hpp>
 #include <utils/Exceptions.hpp>
 
 #include <sstream>
@@ -60,3 +61,5 @@ MeshDataManager::getMeshData(const MeshType& mesh)
 template MeshData<Mesh<1>>& MeshDataManager::getMeshData(const Mesh<1>&);
 template MeshData<Mesh<2>>& MeshDataManager::getMeshData(const Mesh<2>&);
 template MeshData<Mesh<3>>& MeshDataManager::getMeshData(const Mesh<3>&);
+
+template MeshData<PolynomialMesh<2>>& MeshDataManager::getMeshData(const PolynomialMesh<2>&);
diff --git a/src/mesh/MeshDataVariant.hpp b/src/mesh/MeshDataVariant.hpp
index 1f5e77325df2900b5fe0969f4d71347980470217..3edd8a9d1b13155c96b9f0beeeb3e15e929adf61 100644
--- a/src/mesh/MeshDataVariant.hpp
+++ b/src/mesh/MeshDataVariant.hpp
@@ -16,12 +16,17 @@
 template <size_t Dimension>
 class Mesh;
 
+template <size_t Dimension>
+class PolynomialMesh;
+
 class MeshDataVariant
 {
  private:
   using Variant = std::variant<std::shared_ptr<MeshData<Mesh<1>>>,   //
                                std::shared_ptr<MeshData<Mesh<2>>>,   //
-                               std::shared_ptr<MeshData<Mesh<3>>>>;
+                               std::shared_ptr<MeshData<Mesh<3>>>,
+
+                               std::shared_ptr<MeshData<PolynomialMesh<2>>>>;
 
   Variant m_p_mesh_data_variant;
 
diff --git a/src/mesh/MeshSmoother.cpp b/src/mesh/MeshSmoother.cpp
index 1157f5acda265f1c0b9eb3783d7fe91749a5984d..9c6b9b6bc15fa3e2ee0b766fe16eb862c7d7c630 100644
--- a/src/mesh/MeshSmoother.cpp
+++ b/src/mesh/MeshSmoother.cpp
@@ -245,8 +245,8 @@ class MeshSmootherHandler::MeshSmoother
     is_displaced.fill(false);
 
     for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) {
-      MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]);
-      const auto cell_list              = cell_zone.cellList();
+      MeshCellZone cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]);
+      const auto cell_list   = cell_zone.cellList();
       CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
       is_zone_cell.fill(false);
       parallel_for(
diff --git a/src/mesh/MeshTransformer.cpp b/src/mesh/MeshTransformer.cpp
index 13f6f79e4a9d251964fface2b13503e2cd384f8f..9e8c3ce0f3eb4c100c0a6a73808c29d093992402 100644
--- a/src/mesh/MeshTransformer.cpp
+++ b/src/mesh/MeshTransformer.cpp
@@ -56,11 +56,7 @@ MeshTransformer::transform(const FunctionSymbolId& function_id, std::shared_ptr<
       using MeshType             = mesh_type_t<decltype(mesh)>;
       constexpr size_t Dimension = MeshType::Dimension;
       using TransformT           = TinyVector<Dimension>(TinyVector<Dimension>);
-      // if constexpr (is_polygonal_mesh_v<MeshType>) {
       return std::make_shared<MeshVariant>(MeshTransformation<TransformT>::transform(function_id, *mesh));
-      // } else {
-      //   throw NotImplementedError("unexpected mesh type " + demangle<MeshType>());
-      // }
     },
     mesh_v->variant());
 }
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index efe1fcc062ad5645c4e8900371becc5bfcb1b702..271b5783f9d68ffb4c22ee7243a35fb6ed8bf859 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -193,12 +193,7 @@ DiscreteFunctionInterpoler::interpolate() const
   return std::visit(
     [&](auto&& mesh) -> DiscreteFunctionVariant {
       using MeshType = mesh_type_t<decltype(mesh)>;
-
-      if constexpr (is_polygonal_mesh_v<MeshType>) {
-        return this->_interpolate<MeshType>();
-      } else {
-        throw NormalError("unexpected mesh type " + demangle<MeshType>());
-      }
+      return this->_interpolate<MeshType>();
     },
     m_mesh->variant());
 }
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index dc22547f0cc357ede017ed185939cebf8e3f760c..257e5f161ef241829f84c07649ca1eb088d2b81f 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -716,18 +716,13 @@ class DiscreteFunctionP0
       [&](auto&& p_mesh) -> DataType {
         const auto& mesh = *p_mesh;
 
-        using MeshType = mesh_type_t<decltype(mesh)>;
-        if constexpr (is_polygonal_mesh_v<MeshType>) {
-          CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(mesh).Vj();
-
-          CellValue<std::remove_const_t<DataType>> f_v{mesh.connectivity()};
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { f_v[cell_id] = cell_volume[cell_id] * f[cell_id]; });
-
-          return sum(f_v);
-        } else {
-          throw NotImplementedError("unexpected mesh type " + demangle<MeshType>());
-        }
+        CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+        CellValue<std::remove_const_t<DataType>> f_v{mesh.connectivity()};
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { f_v[cell_id] = cell_volume[cell_id] * f[cell_id]; });
+
+        return sum(f_v);
       },
       f.m_mesh_v->variant());
   }