From 3368ffffa00937c9e7d38eaebc4b1e71111122f1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 7 Mar 2024 14:43:50 +0100
Subject: [PATCH] Generalize the use of MeshConcept

---
 src/language/utils/IntegrateCellArray.hpp     |  7 +++---
 src/language/utils/IntegrateCellValue.hpp     |  6 ++---
 src/language/utils/IntegrateOnCells.hpp       | 11 +++++----
 .../ItemArrayVariantFunctionInterpoler.cpp    |  4 ++--
 .../ItemArrayVariantFunctionInterpoler.hpp    |  5 ++--
 .../ItemValueVariantFunctionInterpoler.cpp    |  4 ++--
 .../ItemValueVariantFunctionInterpoler.hpp    |  5 ++--
 src/mesh/DiamondDualMeshBuilder.cpp           |  2 +-
 src/mesh/DiamondDualMeshBuilder.hpp           |  3 ++-
 src/mesh/Dual1DMeshBuilder.cpp                |  2 +-
 src/mesh/Dual1DMeshBuilder.hpp                |  3 ++-
 src/mesh/MedianDualMeshBuilder.cpp            |  2 +-
 src/mesh/MedianDualMeshBuilder.hpp            |  3 ++-
 src/mesh/MeshEdgeBoundary.cpp                 |  6 ++---
 src/mesh/MeshEdgeBoundary.hpp                 | 19 +++++++--------
 src/mesh/MeshEdgeInterface.cpp                |  6 ++---
 src/mesh/MeshEdgeInterface.hpp                | 19 +++++++--------
 src/mesh/MeshFaceBoundary.cpp                 |  4 ++--
 src/mesh/MeshFaceBoundary.hpp                 | 17 +++++++------
 src/mesh/MeshFaceInterface.cpp                |  4 ++--
 src/mesh/MeshFaceInterface.hpp                | 17 +++++++------
 src/mesh/MeshFlatEdgeBoundary.cpp             |  2 +-
 src/mesh/MeshFlatEdgeBoundary.hpp             |  6 ++---
 src/mesh/MeshFlatFaceBoundary.cpp             |  2 +-
 src/mesh/MeshFlatFaceBoundary.hpp             |  6 ++---
 src/mesh/MeshFlatNodeBoundary.cpp             |  4 ++--
 src/mesh/MeshFlatNodeBoundary.hpp             | 17 +++++++------
 src/mesh/MeshLineEdgeBoundary.cpp             |  2 +-
 src/mesh/MeshLineEdgeBoundary.hpp             |  6 ++---
 src/mesh/MeshLineFaceBoundary.cpp             |  2 +-
 src/mesh/MeshLineFaceBoundary.hpp             |  6 ++---
 src/mesh/MeshLineNodeBoundary.cpp             |  4 ++--
 src/mesh/MeshLineNodeBoundary.hpp             |  6 ++---
 src/mesh/MeshNodeBoundary.cpp                 |  8 +++----
 src/mesh/MeshNodeBoundary.hpp                 | 21 ++++++++--------
 src/mesh/MeshNodeBoundaryUtils.hpp            |  3 ++-
 src/mesh/MeshNodeInterface.cpp                |  8 +++----
 src/mesh/MeshNodeInterface.hpp                | 21 ++++++++--------
 src/mesh/MeshRandomizer.cpp                   |  8 +++----
 src/mesh/MeshRandomizer.hpp                   |  4 +++-
 src/mesh/MeshRelaxer.cpp                      |  4 ++--
 src/mesh/MeshRelaxer.hpp                      |  4 +++-
 src/mesh/MeshSmoother.cpp                     |  8 +++----
 src/mesh/MeshSmoother.hpp                     |  4 +++-
 src/mesh/MeshTransformer.cpp                  |  9 +++----
 src/mesh/MeshVariant.hpp                      |  4 ++--
 src/output/GnuplotWriter.cpp                  |  4 ++--
 src/output/GnuplotWriter.hpp                  |  5 ++--
 src/output/GnuplotWriter1D.cpp                |  4 ++--
 src/output/GnuplotWriter1D.hpp                |  5 ++--
 src/output/VTKWriter.cpp                      |  2 +-
 src/output/VTKWriter.hpp                      |  3 ++-
 src/scheme/AcousticSolver.cpp                 | 20 ++++++++--------
 src/scheme/AcousticSolver.hpp                 |  8 ++++---
 src/scheme/CellIntegrator.hpp                 | 12 +++++-----
 src/scheme/DiscreteFunctionIntegrator.cpp     |  8 +++----
 src/scheme/DiscreteFunctionIntegrator.hpp     |  9 +++----
 src/scheme/DiscreteFunctionInterpoler.cpp     |  8 +++----
 src/scheme/DiscreteFunctionInterpoler.hpp     |  9 +++----
 src/scheme/DiscreteFunctionP0.hpp             |  4 ++--
 src/scheme/DiscreteFunctionP0Vector.hpp       |  4 ++--
 .../DiscreteFunctionVectorIntegrator.cpp      |  8 +++----
 .../DiscreteFunctionVectorIntegrator.hpp      |  9 +++----
 .../DiscreteFunctionVectorInterpoler.cpp      |  8 +++----
 .../DiscreteFunctionVectorInterpoler.hpp      |  9 +++----
 src/scheme/EdgeIntegrator.hpp                 | 10 ++++----
 src/scheme/FaceIntegrator.hpp                 | 12 +++++-----
 src/scheme/FluxingAdvectionSolver.cpp         | 24 +++++++++----------
 src/scheme/HyperelasticSolver.cpp             | 20 ++++++++--------
 src/scheme/HyperelasticSolver.hpp             |  8 ++++---
 70 files changed, 275 insertions(+), 256 deletions(-)

diff --git a/src/language/utils/IntegrateCellArray.hpp b/src/language/utils/IntegrateCellArray.hpp
index 41450de3b..380971323 100644
--- a/src/language/utils/IntegrateCellArray.hpp
+++ b/src/language/utils/IntegrateCellArray.hpp
@@ -4,6 +4,7 @@
 #include <language/utils/IntegrateCellValue.hpp>
 #include <mesh/CellType.hpp>
 #include <mesh/ItemArray.hpp>
+#include <mesh/MeshTraits.hpp>
 
 template <typename T>
 class IntegrateCellArray;
@@ -13,7 +14,7 @@ class IntegrateCellArray<OutputType(InputType)>
   static constexpr size_t Dimension = OutputType::Dimension;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static CellArray<OutputType>
   integrate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -33,7 +34,7 @@ class IntegrateCellArray<OutputType(InputType)>
     return cell_array;
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static Table<OutputType>
   integrate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -55,7 +56,7 @@ class IntegrateCellArray<OutputType(InputType)>
     return table;
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static Table<OutputType>
   integrate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/language/utils/IntegrateCellValue.hpp b/src/language/utils/IntegrateCellValue.hpp
index cde3bfacc..d76e9af67 100644
--- a/src/language/utils/IntegrateCellValue.hpp
+++ b/src/language/utils/IntegrateCellValue.hpp
@@ -15,7 +15,7 @@ template <typename OutputType, typename InputType>
 class IntegrateCellValue<OutputType(InputType)>
 {
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static CellValue<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -46,7 +46,7 @@ class IntegrateCellValue<OutputType(InputType)>
       mesh_v->variant());
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static Array<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -76,7 +76,7 @@ class IntegrateCellValue<OutputType(InputType)>
       mesh_v->variant());
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static Array<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index 793575447..2f359b265 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -14,6 +14,7 @@
 #include <mesh/CellType.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemId.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <utils/Array.hpp>
 
 class FunctionSymbolId;
@@ -79,7 +80,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     CellList(const Array<const CellId>& cell_list) : m_cell_list{cell_list} {}
   };
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT>
   static PUGS_INLINE void
   _tensorialIntegrateTo(const FunctionSymbolId& function_symbol_id,
                         const IQuadratureDescriptor& quadrature_descriptor,
@@ -340,7 +341,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     // LCOV_EXCL_STOP
   }
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT>
   static PUGS_INLINE void
   _directIntegrateTo(const FunctionSymbolId& function_symbol_id,
                      const IQuadratureDescriptor& quadrature_descriptor,
@@ -586,7 +587,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
   }
 
  public:
-  template <typename MeshType, typename ArrayT>
+  template <MeshConcept MeshType, typename ArrayT>
   static PUGS_INLINE void
   integrateTo(const FunctionSymbolId& function_symbol_id,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -604,7 +605,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     }
   }
 
-  template <typename MeshType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -623,7 +624,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     return value;
   }
 
-  template <typename MeshType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
index bac8b5c54..a386e8169 100644
--- a/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
+++ b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
@@ -12,7 +12,7 @@
 
 #include <memory>
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 std::shared_ptr<ItemArrayVariant>
 ItemArrayVariantFunctionInterpoler::_interpolate() const
 {
@@ -52,7 +52,7 @@ ItemArrayVariantFunctionInterpoler::_interpolate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::shared_ptr<ItemArrayVariant>
 ItemArrayVariantFunctionInterpoler::_interpolate() const
 {
diff --git a/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp b/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp
index 8631212dc..095614ffe 100644
--- a/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp
+++ b/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp
@@ -5,6 +5,7 @@
 #include <mesh/IZoneDescriptor.hpp>
 #include <mesh/ItemArrayVariant.hpp>
 #include <mesh/ItemType.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 
@@ -15,10 +16,10 @@ class ItemArrayVariantFunctionInterpoler
   const ItemType m_item_type;
   const std::vector<FunctionSymbolId> m_function_id_list;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   std::shared_ptr<ItemArrayVariant> _interpolate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   std::shared_ptr<ItemArrayVariant> _interpolate() const;
 
  public:
diff --git a/src/language/utils/ItemValueVariantFunctionInterpoler.cpp b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
index be073b907..6759e5763 100644
--- a/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
+++ b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
@@ -12,7 +12,7 @@
 
 #include <memory>
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 std::shared_ptr<ItemValueVariant>
 ItemValueVariantFunctionInterpoler::_interpolate() const
 {
@@ -52,7 +52,7 @@ ItemValueVariantFunctionInterpoler::_interpolate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::shared_ptr<ItemValueVariant>
 ItemValueVariantFunctionInterpoler::_interpolate() const
 {
diff --git a/src/language/utils/ItemValueVariantFunctionInterpoler.hpp b/src/language/utils/ItemValueVariantFunctionInterpoler.hpp
index c1c4daddd..046d9bcf5 100644
--- a/src/language/utils/ItemValueVariantFunctionInterpoler.hpp
+++ b/src/language/utils/ItemValueVariantFunctionInterpoler.hpp
@@ -5,6 +5,7 @@
 #include <mesh/IZoneDescriptor.hpp>
 #include <mesh/ItemType.hpp>
 #include <mesh/ItemValueVariant.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 
@@ -15,10 +16,10 @@ class ItemValueVariantFunctionInterpoler
   const ItemType m_item_type;
   const FunctionSymbolId m_function_id;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   std::shared_ptr<ItemValueVariant> _interpolate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   std::shared_ptr<ItemValueVariant> _interpolate() const;
 
  public:
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 43f5528a2..d542cee7f 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -11,7 +11,7 @@
 #include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
 #include <utils/Stringify.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const MeshType& primal_mesh)
 {
diff --git a/src/mesh/DiamondDualMeshBuilder.hpp b/src/mesh/DiamondDualMeshBuilder.hpp
index d1c2b6dbd..ae22c2667 100644
--- a/src/mesh/DiamondDualMeshBuilder.hpp
+++ b/src/mesh/DiamondDualMeshBuilder.hpp
@@ -2,12 +2,13 @@
 #define DIAMOND_DUAL_MESH_BUILDER_HPP
 
 #include <mesh/MeshBuilderBase.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 class DiamondDualMeshBuilder : public MeshBuilderBase
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _buildDualDiamondMeshFrom(const MeshType&);
 
   friend class DualMeshManager;
diff --git a/src/mesh/Dual1DMeshBuilder.cpp b/src/mesh/Dual1DMeshBuilder.cpp
index b39cd000c..64c48423b 100644
--- a/src/mesh/Dual1DMeshBuilder.cpp
+++ b/src/mesh/Dual1DMeshBuilder.cpp
@@ -11,7 +11,7 @@
 #include <mesh/PrimalToDual1DConnectivityDataMapper.hpp>
 #include <utils/Stringify.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 Dual1DMeshBuilder::_buildDual1DMeshFrom(const MeshType& primal_mesh)
 {
diff --git a/src/mesh/Dual1DMeshBuilder.hpp b/src/mesh/Dual1DMeshBuilder.hpp
index 901a46fa5..42810f5c0 100644
--- a/src/mesh/Dual1DMeshBuilder.hpp
+++ b/src/mesh/Dual1DMeshBuilder.hpp
@@ -2,12 +2,13 @@
 #define DUAL_1D_MESH_BUILDER_HPP
 
 #include <mesh/MeshBuilderBase.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 class Dual1DMeshBuilder : public MeshBuilderBase
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _buildDual1DMeshFrom(const MeshType&);
 
   friend class DualMeshManager;
diff --git a/src/mesh/MedianDualMeshBuilder.cpp b/src/mesh/MedianDualMeshBuilder.cpp
index 411005ec9..016f729dc 100644
--- a/src/mesh/MedianDualMeshBuilder.cpp
+++ b/src/mesh/MedianDualMeshBuilder.cpp
@@ -11,7 +11,7 @@
 #include <mesh/PrimalToMedianDualConnectivityDataMapper.hpp>
 #include <utils/Stringify.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 MedianDualMeshBuilder::_buildMedianDualMeshFrom(const MeshType& primal_mesh)
 {
diff --git a/src/mesh/MedianDualMeshBuilder.hpp b/src/mesh/MedianDualMeshBuilder.hpp
index ad31b6aa1..eb33440c5 100644
--- a/src/mesh/MedianDualMeshBuilder.hpp
+++ b/src/mesh/MedianDualMeshBuilder.hpp
@@ -2,12 +2,13 @@
 #define MEDIAN_DUAL_MESH_BUILDER_HPP
 
 #include <mesh/MeshBuilderBase.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 class MedianDualMeshBuilder : public MeshBuilderBase
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _buildMedianDualMeshFrom(const MeshType&);
 
   friend class DualMeshManager;
diff --git a/src/mesh/MeshEdgeBoundary.cpp b/src/mesh/MeshEdgeBoundary.cpp
index 0ad15dbad..65d95e9f6 100644
--- a/src/mesh/MeshEdgeBoundary.cpp
+++ b/src/mesh/MeshEdgeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeBoundary::MeshEdgeBoundary(const MeshType&, const RefEdgeList& ref_edge_list) : m_ref_edge_list(ref_edge_list)
 {}
 
@@ -13,7 +13,7 @@ template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<1>&, const RefEdgeList&);
 template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<2>&, const RefEdgeList&);
 template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<3>&, const RefEdgeList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeBoundary::MeshEdgeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
   constexpr size_t Dimension = MeshType::Dimension;
@@ -58,7 +58,7 @@ MeshEdgeBoundary::MeshEdgeBoundary(const MeshType& mesh, const RefFaceList& ref_
 template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<2>&, const RefFaceList&);
 template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<3>&, const RefFaceList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeBoundary
 getMeshEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshEdgeBoundary.hpp b/src/mesh/MeshEdgeBoundary.hpp
index c65106b2d..411de1885 100644
--- a/src/mesh/MeshEdgeBoundary.hpp
+++ b/src/mesh/MeshEdgeBoundary.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,41 +13,39 @@ class [[nodiscard]] MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
   RefEdgeList m_ref_edge_list;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   friend MeshEdgeBoundary getMeshEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
   MeshEdgeBoundary& operator=(const MeshEdgeBoundary&) = default;
-  MeshEdgeBoundary& operator=(MeshEdgeBoundary&&)      = default;
+  MeshEdgeBoundary& operator=(MeshEdgeBoundary&&) = default;
 
   PUGS_INLINE
-  const RefEdgeList&
-  refEdgeList() const
+  const RefEdgeList& refEdgeList() const
   {
     return m_ref_edge_list;
   }
 
   PUGS_INLINE
-  const Array<const EdgeId>&
-  edgeList() const
+  const Array<const EdgeId>& edgeList() const
   {
     return m_ref_edge_list.list();
   }
 
  protected:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshEdgeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshEdgeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list);
 
  public:
   MeshEdgeBoundary(const MeshEdgeBoundary&) = default;   // LCOV_EXCL_LINE
-  MeshEdgeBoundary(MeshEdgeBoundary&&)      = default;   // LCOV_EXCL_LINE
+  MeshEdgeBoundary(MeshEdgeBoundary &&)     = default;   // LCOV_EXCL_LINE
 
   MeshEdgeBoundary()          = default;
   virtual ~MeshEdgeBoundary() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeBoundary getMeshEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_EDGE_BOUNDARY_HPP
diff --git a/src/mesh/MeshEdgeInterface.cpp b/src/mesh/MeshEdgeInterface.cpp
index 55235d8c7..b74964eee 100644
--- a/src/mesh/MeshEdgeInterface.cpp
+++ b/src/mesh/MeshEdgeInterface.cpp
@@ -5,7 +5,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeInterface::MeshEdgeInterface(const MeshType&, const RefEdgeList& ref_edge_list) : m_ref_edge_list(ref_edge_list)
 {}
 
@@ -13,7 +13,7 @@ template MeshEdgeInterface::MeshEdgeInterface(const Mesh<1>&, const RefEdgeList&
 template MeshEdgeInterface::MeshEdgeInterface(const Mesh<2>&, const RefEdgeList&);
 template MeshEdgeInterface::MeshEdgeInterface(const Mesh<3>&, const RefEdgeList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeInterface::MeshEdgeInterface(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
   constexpr size_t Dimension = MeshType::Dimension;
@@ -58,7 +58,7 @@ MeshEdgeInterface::MeshEdgeInterface(const MeshType& mesh, const RefFaceList& re
 template MeshEdgeInterface::MeshEdgeInterface(const Mesh<2>&, const RefFaceList&);
 template MeshEdgeInterface::MeshEdgeInterface(const Mesh<3>&, const RefFaceList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeInterface
 getMeshEdgeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor)
 {
diff --git a/src/mesh/MeshEdgeInterface.hpp b/src/mesh/MeshEdgeInterface.hpp
index 45bcdbe84..b4a662d04 100644
--- a/src/mesh/MeshEdgeInterface.hpp
+++ b/src/mesh/MeshEdgeInterface.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IInterfaceDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,41 +13,39 @@ class [[nodiscard]] MeshEdgeInterface   // clazy:exclude=copyable-polymorphic
   RefEdgeList m_ref_edge_list;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   friend MeshEdgeInterface getMeshEdgeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
   MeshEdgeInterface& operator=(const MeshEdgeInterface&) = default;
-  MeshEdgeInterface& operator=(MeshEdgeInterface&&)      = default;
+  MeshEdgeInterface& operator=(MeshEdgeInterface&&) = default;
 
   PUGS_INLINE
-  const RefEdgeList&
-  refEdgeList() const
+  const RefEdgeList& refEdgeList() const
   {
     return m_ref_edge_list;
   }
 
   PUGS_INLINE
-  const Array<const EdgeId>&
-  edgeList() const
+  const Array<const EdgeId>& edgeList() const
   {
     return m_ref_edge_list.list();
   }
 
  protected:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshEdgeInterface(const MeshType& mesh, const RefEdgeList& ref_edge_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshEdgeInterface(const MeshType& mesh, const RefFaceList& ref_face_list);
 
  public:
   MeshEdgeInterface(const MeshEdgeInterface&) = default;   // LCOV_EXCL_LINE
-  MeshEdgeInterface(MeshEdgeInterface&&)      = default;   // LCOV_EXCL_LINE
+  MeshEdgeInterface(MeshEdgeInterface &&)     = default;   // LCOV_EXCL_LINE
 
   MeshEdgeInterface()          = default;
   virtual ~MeshEdgeInterface() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeInterface getMeshEdgeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
 #endif   // MESH_EDGE_INTERFACE_HPP
diff --git a/src/mesh/MeshFaceBoundary.cpp b/src/mesh/MeshFaceBoundary.cpp
index d883ba63a..9d0598612 100644
--- a/src/mesh/MeshFaceBoundary.cpp
+++ b/src/mesh/MeshFaceBoundary.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceBoundary::MeshFaceBoundary(const MeshType&, const RefFaceList& ref_face_list) : m_ref_face_list(ref_face_list)
 {}
 
@@ -12,7 +12,7 @@ template MeshFaceBoundary::MeshFaceBoundary(const Mesh<1>&, const RefFaceList&);
 template MeshFaceBoundary::MeshFaceBoundary(const Mesh<2>&, const RefFaceList&);
 template MeshFaceBoundary::MeshFaceBoundary(const Mesh<3>&, const RefFaceList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceBoundary
 getMeshFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFaceBoundary.hpp b/src/mesh/MeshFaceBoundary.hpp
index c44aafa66..443f868c0 100644
--- a/src/mesh/MeshFaceBoundary.hpp
+++ b/src/mesh/MeshFaceBoundary.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,39 +13,37 @@ class [[nodiscard]] MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
   RefFaceList m_ref_face_list;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshFaceBoundary getMeshFaceBoundary(const MeshTypeT& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
   MeshFaceBoundary& operator=(const MeshFaceBoundary&) = default;
-  MeshFaceBoundary& operator=(MeshFaceBoundary&&)      = default;
+  MeshFaceBoundary& operator=(MeshFaceBoundary&&) = default;
 
   PUGS_INLINE
-  const RefFaceList&
-  refFaceList() const
+  const RefFaceList& refFaceList() const
   {
     return m_ref_face_list;
   }
 
   PUGS_INLINE
-  const Array<const FaceId>&
-  faceList() const
+  const Array<const FaceId>& faceList() const
   {
     return m_ref_face_list.list();
   }
 
  protected:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshFaceBoundary(const MeshType& mesh, const RefFaceList& ref_face_list);
 
  public:
   MeshFaceBoundary(const MeshFaceBoundary&) = default;   // LCOV_EXCL_LINE
-  MeshFaceBoundary(MeshFaceBoundary&&)      = default;   // LCOV_EXCL_LINE
+  MeshFaceBoundary(MeshFaceBoundary &&)     = default;   // LCOV_EXCL_LINE
 
   MeshFaceBoundary()          = default;
   virtual ~MeshFaceBoundary() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceBoundary getMeshFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_FACE_BOUNDARY_HPP
diff --git a/src/mesh/MeshFaceInterface.cpp b/src/mesh/MeshFaceInterface.cpp
index 663ac8d44..8dfe8cb49 100644
--- a/src/mesh/MeshFaceInterface.cpp
+++ b/src/mesh/MeshFaceInterface.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceInterface::MeshFaceInterface(const MeshType&, const RefFaceList& ref_face_list) : m_ref_face_list(ref_face_list)
 {}
 
@@ -12,7 +12,7 @@ template MeshFaceInterface::MeshFaceInterface(const Mesh<1>&, const RefFaceList&
 template MeshFaceInterface::MeshFaceInterface(const Mesh<2>&, const RefFaceList&);
 template MeshFaceInterface::MeshFaceInterface(const Mesh<3>&, const RefFaceList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceInterface
 getMeshFaceInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor)
 {
diff --git a/src/mesh/MeshFaceInterface.hpp b/src/mesh/MeshFaceInterface.hpp
index 0d8277e0b..64c75e14a 100644
--- a/src/mesh/MeshFaceInterface.hpp
+++ b/src/mesh/MeshFaceInterface.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IInterfaceDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,39 +13,37 @@ class [[nodiscard]] MeshFaceInterface   // clazy:exclude=copyable-polymorphic
   RefFaceList m_ref_face_list;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   friend MeshFaceInterface getMeshFaceInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
   MeshFaceInterface& operator=(const MeshFaceInterface&) = default;
-  MeshFaceInterface& operator=(MeshFaceInterface&&)      = default;
+  MeshFaceInterface& operator=(MeshFaceInterface&&) = default;
 
   PUGS_INLINE
-  const RefFaceList&
-  refFaceList() const
+  const RefFaceList& refFaceList() const
   {
     return m_ref_face_list;
   }
 
   PUGS_INLINE
-  const Array<const FaceId>&
-  faceList() const
+  const Array<const FaceId>& faceList() const
   {
     return m_ref_face_list.list();
   }
 
  protected:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshFaceInterface(const MeshType& mesh, const RefFaceList& ref_face_list);
 
  public:
   MeshFaceInterface(const MeshFaceInterface&) = default;   // LCOV_EXCL_LINE
-  MeshFaceInterface(MeshFaceInterface&&)      = default;   // LCOV_EXCL_LINE
+  MeshFaceInterface(MeshFaceInterface &&)     = default;   // LCOV_EXCL_LINE
 
   MeshFaceInterface()          = default;
   virtual ~MeshFaceInterface() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceInterface getMeshFaceInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
 #endif   // MESH_FACE_INTERFACE_HPP
diff --git a/src/mesh/MeshFlatEdgeBoundary.cpp b/src/mesh/MeshFlatEdgeBoundary.cpp
index 7a2d4918d..c20b9c1a5 100644
--- a/src/mesh/MeshFlatEdgeBoundary.cpp
+++ b/src/mesh/MeshFlatEdgeBoundary.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatEdgeBoundary<MeshType>
 getMeshFlatEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFlatEdgeBoundary.hpp b/src/mesh/MeshFlatEdgeBoundary.hpp
index f395c4ebf..10ab8edb8 100644
--- a/src/mesh/MeshFlatEdgeBoundary.hpp
+++ b/src/mesh/MeshFlatEdgeBoundary.hpp
@@ -3,7 +3,7 @@
 
 #include <mesh/MeshEdgeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -24,7 +24,7 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
   MeshFlatEdgeBoundary& operator=(const MeshFlatEdgeBoundary&) = default;
   MeshFlatEdgeBoundary& operator=(MeshFlatEdgeBoundary&&)      = default;
 
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshFlatEdgeBoundary<MeshTypeT> getMeshFlatEdgeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -40,7 +40,7 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
   ~MeshFlatEdgeBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatEdgeBoundary<MeshType> getMeshFlatEdgeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshFlatFaceBoundary.cpp b/src/mesh/MeshFlatFaceBoundary.cpp
index 2ecdf6768..632c6a033 100644
--- a/src/mesh/MeshFlatFaceBoundary.cpp
+++ b/src/mesh/MeshFlatFaceBoundary.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatFaceBoundary<MeshType>
 getMeshFlatFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFlatFaceBoundary.hpp b/src/mesh/MeshFlatFaceBoundary.hpp
index 1aa2e854f..0c53a31e0 100644
--- a/src/mesh/MeshFlatFaceBoundary.hpp
+++ b/src/mesh/MeshFlatFaceBoundary.hpp
@@ -3,7 +3,7 @@
 
 #include <mesh/MeshFaceBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -24,7 +24,7 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
   MeshFlatFaceBoundary& operator=(const MeshFlatFaceBoundary&) = default;
   MeshFlatFaceBoundary& operator=(MeshFlatFaceBoundary&&)      = default;
 
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshFlatFaceBoundary<MeshTypeT> getMeshFlatFaceBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -40,7 +40,7 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
   ~MeshFlatFaceBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatFaceBoundary<MeshType> getMeshFlatFaceBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index 983f96738..bdc7c20a2 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshNodeBoundaryUtils.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 MeshFlatNodeBoundary<MeshType>::_checkBoundaryIsFlat(const TinyVector<MeshType::Dimension, double>& normal,
                                                      const TinyVector<MeshType::Dimension, double>& origin,
@@ -323,7 +323,7 @@ MeshFlatNodeBoundary<Mesh<3>>::_getOutgoingNormal(const Mesh<3>& mesh)
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatNodeBoundary<MeshType>
 getMeshFlatNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFlatNodeBoundary.hpp b/src/mesh/MeshFlatNodeBoundary.hpp
index 453804724..2eb28d148 100644
--- a/src/mesh/MeshFlatNodeBoundary.hpp
+++ b/src/mesh/MeshFlatNodeBoundary.hpp
@@ -2,8 +2,9 @@
 #define MESH_FLAT_NODE_BOUNDARY_HPP
 
 #include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -19,23 +20,21 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
   Rd _getNormal(const MeshType& mesh);
 
   void _checkBoundaryIsFlat(const TinyVector<MeshType::Dimension, double>& normal,
-                            const TinyVector<MeshType::Dimension, double>& origin,
-                            const double length,
+                            const TinyVector<MeshType::Dimension, double>& origin, const double length,
                             const MeshType& mesh) const;
 
   Rd _getOutgoingNormal(const MeshType& mesh);
 
  public:
-  const Rd&
-  outgoingNormal() const
+  const Rd& outgoingNormal() const
   {
     return m_outgoing_normal;
   }
 
   MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&)      = default;
+  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
 
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshFlatNodeBoundary<MeshTypeT> getMeshFlatNodeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -51,11 +50,11 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
  public:
   MeshFlatNodeBoundary()                            = default;
   MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary(MeshFlatNodeBoundary&&)      = default;
+  MeshFlatNodeBoundary(MeshFlatNodeBoundary &&)     = default;
   ~MeshFlatNodeBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatNodeBoundary<MeshType> getMeshFlatNodeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshLineEdgeBoundary.cpp b/src/mesh/MeshLineEdgeBoundary.cpp
index f1e0333a9..40de23d1e 100644
--- a/src/mesh/MeshLineEdgeBoundary.cpp
+++ b/src/mesh/MeshLineEdgeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineEdgeBoundary<MeshType>
 getMeshLineEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshLineEdgeBoundary.hpp b/src/mesh/MeshLineEdgeBoundary.hpp
index 39c97a194..0fc5db0ba 100644
--- a/src/mesh/MeshLineEdgeBoundary.hpp
+++ b/src/mesh/MeshLineEdgeBoundary.hpp
@@ -4,7 +4,7 @@
 #include <algebra/TinyMatrix.hpp>
 #include <mesh/MeshEdgeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -18,7 +18,7 @@ class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // cl
   const Rd m_direction;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshLineEdgeBoundary<MeshTypeT> getMeshLineEdgeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -44,7 +44,7 @@ class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // cl
   ~MeshLineEdgeBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineEdgeBoundary<MeshType> getMeshLineEdgeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshLineFaceBoundary.cpp b/src/mesh/MeshLineFaceBoundary.cpp
index 81ea991ed..3608e7aef 100644
--- a/src/mesh/MeshLineFaceBoundary.cpp
+++ b/src/mesh/MeshLineFaceBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineFaceBoundary<MeshType>
 getMeshLineFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshLineFaceBoundary.hpp b/src/mesh/MeshLineFaceBoundary.hpp
index 203faac8e..206961a08 100644
--- a/src/mesh/MeshLineFaceBoundary.hpp
+++ b/src/mesh/MeshLineFaceBoundary.hpp
@@ -4,7 +4,7 @@
 #include <algebra/TinyMatrix.hpp>
 #include <mesh/MeshFaceBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -17,7 +17,7 @@ class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // cl
   const Rd m_direction;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshLineFaceBoundary<MeshTypeT> getMeshLineFaceBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -43,7 +43,7 @@ class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // cl
   ~MeshLineFaceBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineFaceBoundary<MeshType> getMeshLineFaceBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshLineNodeBoundary.cpp b/src/mesh/MeshLineNodeBoundary.cpp
index 1228fa5aa..1cd3bb64d 100644
--- a/src/mesh/MeshLineNodeBoundary.cpp
+++ b/src/mesh/MeshLineNodeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshNodeBoundaryUtils.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 MeshLineNodeBoundary<MeshType>::_checkBoundaryIsLine(const TinyVector<MeshType::Dimension, double>& direction,
                                                      const TinyVector<MeshType::Dimension, double>& origin,
@@ -105,7 +105,7 @@ MeshLineNodeBoundary<Mesh<3>>::_getDirection(const Mesh<3>& mesh)
   return direction;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineNodeBoundary<MeshType>
 getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshLineNodeBoundary.hpp b/src/mesh/MeshLineNodeBoundary.hpp
index bb6fc4b04..fb75ef9a4 100644
--- a/src/mesh/MeshLineNodeBoundary.hpp
+++ b/src/mesh/MeshLineNodeBoundary.hpp
@@ -4,7 +4,7 @@
 #include <algebra/TinyMatrix.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -24,7 +24,7 @@ class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // cl
                             const MeshType& mesh) const;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshLineNodeBoundary<MeshTypeT> getMeshLineNodeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -58,7 +58,7 @@ class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // cl
   ~MeshLineNodeBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineNodeBoundary<MeshType> getMeshLineNodeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshNodeBoundary.cpp b/src/mesh/MeshNodeBoundary.cpp
index e750580ce..8cdcad2fb 100644
--- a/src/mesh/MeshNodeBoundary.cpp
+++ b/src/mesh/MeshNodeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
   constexpr size_t Dimension = MeshType::Dimension;
@@ -52,7 +52,7 @@ MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
 {
   constexpr size_t Dimension = MeshType::Dimension;
@@ -98,7 +98,7 @@ MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeBoundary::MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list) : m_ref_node_list(ref_node_list)
 {
   if (ref_node_list.type() != RefItemListBase::Type::boundary) {
@@ -121,7 +121,7 @@ template MeshNodeBoundary::MeshNodeBoundary(const Mesh<1>&, const RefNodeList&);
 template MeshNodeBoundary::MeshNodeBoundary(const Mesh<2>&, const RefNodeList&);
 template MeshNodeBoundary::MeshNodeBoundary(const Mesh<3>&, const RefNodeList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeBoundary
 getMeshNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 0e869e2ee..36d79a308 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,43 +13,41 @@ class [[nodiscard]] MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   RefNodeList m_ref_node_list;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshNodeBoundary getMeshNodeBoundary(const MeshTypeT& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
   MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
-  MeshNodeBoundary& operator=(MeshNodeBoundary&&)      = default;
+  MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
 
   PUGS_INLINE
-  const RefNodeList&
-  refNodeList() const
+  const RefNodeList& refNodeList() const
   {
     return m_ref_node_list;
   }
 
   PUGS_INLINE
-  const Array<const NodeId>&
-  nodeList() const
+  const Array<const NodeId>& nodeList() const
   {
     return m_ref_node_list.list();
   }
 
  protected:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list);
 
  public:
   MeshNodeBoundary(const MeshNodeBoundary&) = default;
-  MeshNodeBoundary(MeshNodeBoundary&&)      = default;
+  MeshNodeBoundary(MeshNodeBoundary &&)     = default;
 
   MeshNodeBoundary()          = default;
   virtual ~MeshNodeBoundary() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeBoundary getMeshNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshNodeBoundaryUtils.hpp b/src/mesh/MeshNodeBoundaryUtils.hpp
index a62f9d8f3..f187a210d 100644
--- a/src/mesh/MeshNodeBoundaryUtils.hpp
+++ b/src/mesh/MeshNodeBoundaryUtils.hpp
@@ -2,11 +2,12 @@
 #define MESH_NODE_BOUNDARY_UTILS_HPP
 
 #include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 
 #include <array>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::array<TinyVector<MeshType::Dimension>, MeshType::Dimension*(MeshType::Dimension - 1)> getBounds(
   const MeshType& mesh,
   const RefNodeList& ref_node_list);
diff --git a/src/mesh/MeshNodeInterface.cpp b/src/mesh/MeshNodeInterface.cpp
index 54da43007..57eed14cb 100644
--- a/src/mesh/MeshNodeInterface.cpp
+++ b/src/mesh/MeshNodeInterface.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshTraits.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface::MeshNodeInterface(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
   static_assert(is_polygonal_mesh_v<MeshType>);
@@ -53,7 +53,7 @@ MeshNodeInterface::MeshNodeInterface(const MeshType& mesh, const RefFaceList& re
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface::MeshNodeInterface(const MeshType& mesh, const RefEdgeList& ref_edge_list)
 {
   static_assert(is_polygonal_mesh_v<MeshType>);
@@ -100,7 +100,7 @@ MeshNodeInterface::MeshNodeInterface(const MeshType& mesh, const RefEdgeList& re
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface::MeshNodeInterface(const MeshType&, const RefNodeList& ref_node_list) : m_ref_node_list(ref_node_list)
 {
   if (ref_node_list.type() != RefItemListBase::Type::interface) {
@@ -123,7 +123,7 @@ template MeshNodeInterface::MeshNodeInterface(const Mesh<1>&, const RefNodeList&
 template MeshNodeInterface::MeshNodeInterface(const Mesh<2>&, const RefNodeList&);
 template MeshNodeInterface::MeshNodeInterface(const Mesh<3>&, const RefNodeList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface
 getMeshNodeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor)
 {
diff --git a/src/mesh/MeshNodeInterface.hpp b/src/mesh/MeshNodeInterface.hpp
index 4eaafd06b..bb2105649 100644
--- a/src/mesh/MeshNodeInterface.hpp
+++ b/src/mesh/MeshNodeInterface.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IInterfaceDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,43 +13,41 @@ class [[nodiscard]] MeshNodeInterface   // clazy:exclude=copyable-polymorphic
   RefNodeList m_ref_node_list;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   friend MeshNodeInterface getMeshNodeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
   MeshNodeInterface& operator=(const MeshNodeInterface&) = default;
-  MeshNodeInterface& operator=(MeshNodeInterface&&)      = default;
+  MeshNodeInterface& operator=(MeshNodeInterface&&) = default;
 
   PUGS_INLINE
-  const RefNodeList&
-  refNodeList() const
+  const RefNodeList& refNodeList() const
   {
     return m_ref_node_list;
   }
 
   PUGS_INLINE
-  const Array<const NodeId>&
-  nodeList() const
+  const Array<const NodeId>& nodeList() const
   {
     return m_ref_node_list.list();
   }
 
  protected:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeInterface(const MeshType& mesh, const RefFaceList& ref_face_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeInterface(const MeshType& mesh, const RefEdgeList& ref_edge_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeInterface(const MeshType&, const RefNodeList& ref_node_list);
 
  public:
   MeshNodeInterface(const MeshNodeInterface&) = default;
-  MeshNodeInterface(MeshNodeInterface&&)      = default;
+  MeshNodeInterface(MeshNodeInterface &&)     = default;
 
   MeshNodeInterface()          = default;
   virtual ~MeshNodeInterface() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface getMeshNodeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
 #endif   // MESH_NODE_INTERFACE_HPP
diff --git a/src/mesh/MeshRandomizer.cpp b/src/mesh/MeshRandomizer.cpp
index d6b8c2f1b..8c94807db 100644
--- a/src/mesh/MeshRandomizer.cpp
+++ b/src/mesh/MeshRandomizer.cpp
@@ -16,7 +16,7 @@
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <utils/RandomEngine.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshRandomizerHandler::MeshRandomizer
 {
  private:
@@ -274,7 +274,7 @@ class MeshRandomizerHandler::MeshRandomizer
   ~MeshRandomizer() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshRandomizerHandler::MeshRandomizer<MeshType>::AxisBoundaryCondition
 {
  public:
@@ -313,7 +313,7 @@ class MeshRandomizerHandler::MeshRandomizer<Mesh<1>>::AxisBoundaryCondition
   ~AxisBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshRandomizerHandler::MeshRandomizer<MeshType>::FixedBoundaryCondition
 {
  private:
@@ -331,7 +331,7 @@ class MeshRandomizerHandler::MeshRandomizer<MeshType>::FixedBoundaryCondition
   ~FixedBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshRandomizerHandler::MeshRandomizer<MeshType>::SymmetryBoundaryCondition
 {
  public:
diff --git a/src/mesh/MeshRandomizer.hpp b/src/mesh/MeshRandomizer.hpp
index 048e04bfd..6db719cb8 100644
--- a/src/mesh/MeshRandomizer.hpp
+++ b/src/mesh/MeshRandomizer.hpp
@@ -1,6 +1,8 @@
 #ifndef MESH_RANDOMIZER_HPP
 #define MESH_RANDOMIZER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 #include <memory>
 #include <vector>
 
@@ -11,7 +13,7 @@ class FunctionSymbolId;
 class MeshRandomizerHandler
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   class MeshRandomizer;
 
  public:
diff --git a/src/mesh/MeshRelaxer.cpp b/src/mesh/MeshRelaxer.cpp
index 7dcc86bc0..65302649d 100644
--- a/src/mesh/MeshRelaxer.cpp
+++ b/src/mesh/MeshRelaxer.cpp
@@ -5,12 +5,12 @@
 #include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::shared_ptr<const MeshVariant>
 MeshRelaxer::_relax(const MeshType& source_mesh, const MeshType& destination_mesh, const double& theta) const
 {
   static_assert(is_polygonal_mesh_v<MeshType>);
-  if (source_mesh.shared_connectivity() == destination_mesh.shared_connectivity()) {
+  if (source_mesh.connectivity().id() == destination_mesh.connectivity().id()) {
     using ConnectivityType               = typename MeshType::Connectivity;
     constexpr size_t Dimension           = MeshType::Dimension;
     const ConnectivityType& connectivity = source_mesh.connectivity();
diff --git a/src/mesh/MeshRelaxer.hpp b/src/mesh/MeshRelaxer.hpp
index 1e2b08206..c771842cf 100644
--- a/src/mesh/MeshRelaxer.hpp
+++ b/src/mesh/MeshRelaxer.hpp
@@ -1,6 +1,8 @@
 #ifndef MESH_RELAXER_HPP
 #define MESH_RELAXER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 class MeshVariant;
 
 #include <memory>
@@ -8,7 +10,7 @@ class MeshVariant;
 class MeshRelaxer
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   std::shared_ptr<const MeshVariant> _relax(const MeshType& source_mesh,
                                             const MeshType& destination_mesh,
                                             const double& theta) const;
diff --git a/src/mesh/MeshSmoother.cpp b/src/mesh/MeshSmoother.cpp
index eafc9b94b..7950c732a 100644
--- a/src/mesh/MeshSmoother.cpp
+++ b/src/mesh/MeshSmoother.cpp
@@ -20,7 +20,7 @@
 
 #include <variant>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshSmootherHandler::MeshSmoother
 {
  private:
@@ -321,7 +321,7 @@ class MeshSmootherHandler::MeshSmoother
   ~MeshSmoother() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshSmootherHandler::MeshSmoother<MeshType>::AxisBoundaryCondition
 {
  public:
@@ -360,7 +360,7 @@ class MeshSmootherHandler::MeshSmoother<Mesh<1>>::AxisBoundaryCondition
   ~AxisBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshSmootherHandler::MeshSmoother<MeshType>::FixedBoundaryCondition
 {
  private:
@@ -378,7 +378,7 @@ class MeshSmootherHandler::MeshSmoother<MeshType>::FixedBoundaryCondition
   ~FixedBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshSmootherHandler::MeshSmoother<MeshType>::SymmetryBoundaryCondition
 {
  public:
diff --git a/src/mesh/MeshSmoother.hpp b/src/mesh/MeshSmoother.hpp
index c5875561a..d43e380a9 100644
--- a/src/mesh/MeshSmoother.hpp
+++ b/src/mesh/MeshSmoother.hpp
@@ -1,6 +1,8 @@
 #ifndef MESH_SMOOTHER_HPP
 #define MESH_SMOOTHER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 #include <memory>
 #include <vector>
 
@@ -13,7 +15,7 @@ class DiscreteFunctionVariant;
 class MeshSmootherHandler
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   class MeshSmoother;
 
  public:
diff --git a/src/mesh/MeshTransformer.cpp b/src/mesh/MeshTransformer.cpp
index d44d3251c..81f723a6c 100644
--- a/src/mesh/MeshTransformer.cpp
+++ b/src/mesh/MeshTransformer.cpp
@@ -1,6 +1,7 @@
 #include <mesh/MeshTransformer.hpp>
 
 #include <mesh/Connectivity.hpp>
+#include <mesh/ItemArray.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
@@ -13,15 +14,15 @@ class MeshTransformer::MeshTransformation<OutputType(InputType)>
   static constexpr size_t Dimension = OutputType::Dimension;
 
  public:
-  template <typename MeshType>
-  static std::shared_ptr<const MeshType>
-  transform(const FunctionSymbolId& function_symbol_id, const MeshType& mesh)
+  template <size_t Dimension>
+  static std::shared_ptr<const Mesh<Dimension>>
+  transform(const FunctionSymbolId& function_symbol_id, const Mesh<Dimension>& mesh)
   {
     NodeValue<OutputType> xr(mesh.connectivity());
     NodeValue<const InputType> given_xr = mesh.xr();
     EvaluateAtPoints<OutputType(InputType)>::evaluateTo(function_symbol_id, given_xr, xr);
 
-    return std::make_shared<const MeshType>(mesh.shared_connectivity(), xr);
+    return std::make_shared<const Mesh<Dimension>>(mesh.shared_connectivity(), xr);
   }
 };
 
diff --git a/src/mesh/MeshVariant.hpp b/src/mesh/MeshVariant.hpp
index 7278878a9..c37d2dd01 100644
--- a/src/mesh/MeshVariant.hpp
+++ b/src/mesh/MeshVariant.hpp
@@ -40,7 +40,7 @@ class MeshVariant
   const IConnectivity& connectivity() const;
   std::shared_ptr<const IConnectivity> shared_connectivity() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE std::shared_ptr<const MeshType>
   get() const
   {
@@ -66,7 +66,7 @@ class MeshVariant
 
   MeshVariant() = delete;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshVariant(const std::shared_ptr<const MeshType>& p_mesh) : m_p_mesh_variant{p_mesh}
   {}
 
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 5415d59f8..991ec7453 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -148,7 +148,7 @@ GnuplotWriter::_writeData(const ItemDataT& item_data,
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 GnuplotWriter::_writeDataAtNodes(const MeshType& mesh,
                                  const OutputNamedItemDataSet& output_named_item_data_set,
@@ -240,7 +240,7 @@ GnuplotWriter::_writeNodeData(const NodeArray<DataType>& node_array, NodeId node
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 GnuplotWriter::_write(const MeshType& mesh,
                       const OutputNamedItemDataSet& output_named_item_data_set,
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
index 9cf22a7ef..3069b79da 100644
--- a/src/output/GnuplotWriter.hpp
+++ b/src/output/GnuplotWriter.hpp
@@ -5,6 +5,7 @@
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <output/OutputNamedItemValueSet.hpp>
 
 class MeshVariant;
@@ -40,12 +41,12 @@ class GnuplotWriter final : public WriterBase
                   [[maybe_unused]] NodeId node_id,
                   std::ostream& fout) const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _writeDataAtNodes(const MeshType& mesh,
                          const OutputNamedItemDataSet& output_named_item_data_set,
                          std::ostream& fout) const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_data_set,
               std::optional<double> time) const;
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index a31d5432f..84e8aa080 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -142,7 +142,7 @@ GnuplotWriter1D::_itemDataNbRow(const ItemArray<DataType, item_type>& item_array
   return item_array.sizeOfArrays();
 }
 
-template <typename MeshType, ItemType item_type>
+template <MeshConcept MeshType, ItemType item_type>
 void
 GnuplotWriter1D::_writeItemDatas(const MeshType& mesh,
                                  const OutputNamedItemDataSet& output_named_item_data_set,
@@ -275,7 +275,7 @@ GnuplotWriter1D::_writeItemDatas(const MeshType& mesh,
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 GnuplotWriter1D::_write(const MeshType& mesh,
                         const OutputNamedItemDataSet& output_named_item_data_set,
diff --git a/src/output/GnuplotWriter1D.hpp b/src/output/GnuplotWriter1D.hpp
index 61b697dd1..666619182 100644
--- a/src/output/GnuplotWriter1D.hpp
+++ b/src/output/GnuplotWriter1D.hpp
@@ -5,6 +5,7 @@
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <output/OutputNamedItemValueSet.hpp>
 
 class MeshVariant;
@@ -37,14 +38,14 @@ class GnuplotWriter1D final : public WriterBase
   template <typename DataType, ItemType item_type>
   size_t _itemDataNbRow(const ItemArray<DataType, item_type>&) const;
 
-  template <typename MeshType, ItemType item_type>
+  template <MeshConcept MeshType, ItemType item_type>
   void _writeItemDatas(const MeshType& mesh,
                        const OutputNamedItemDataSet& output_named_item_data_set,
                        std::ostream& fout) const;
 
   void _writePreamble(const OutputNamedItemDataSet& output_named_item_value_set, std::ostream& fout) const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_value_set,
               std::optional<double> time) const;
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 889b6b554..cd7834d71 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -355,7 +355,7 @@ VTKWriter::_write_node_data(std::ofstream& os,
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 VTKWriter::_write(const MeshType& mesh,
                   const OutputNamedItemDataSet& given_output_named_item_data_set,
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 1141daec8..6c40f1005 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -5,6 +5,7 @@
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <output/OutputNamedItemValueSet.hpp>
 
 #include <optional>
@@ -93,7 +94,7 @@ class VTKWriter final : public WriterBase
                         const ItemDataT& item_data,
                         SerializedDataList& serialized_data_list) const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_data_set,
               std::optional<double> time) const;
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 43476df57..90998d3ab 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -48,7 +48,7 @@ acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
     c.meshVariant()->variant());
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler::IAcousticSolver
 {
  private:
@@ -514,7 +514,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   ~AcousticSolver()                = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 AcousticSolverHandler::AcousticSolver<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
                                                                   const MeshType& mesh,
@@ -579,7 +579,7 @@ AcousticSolverHandler::AcousticSolver<MeshType>::_applyPressureBC(const Boundary
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 AcousticSolverHandler::AcousticSolver<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
                                                                   NodeValue<Rdxd>& Ar,
@@ -610,7 +610,7 @@ AcousticSolverHandler::AcousticSolver<MeshType>::_applySymmetryBC(const Boundary
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 AcousticSolverHandler::AcousticSolver<MeshType>::_applyVelocityBC(const BoundaryConditionList& bc_list,
                                                                   NodeValue<Rdxd>& Ar,
@@ -647,7 +647,7 @@ AcousticSolverHandler::AcousticSolver<MeshType>::_applyVelocityBC(const Boundary
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 AcousticSolverHandler::AcousticSolver<MeshType>::_applyExternalVelocityBC(const BoundaryConditionList& bc_list,
                                                                           const DiscreteScalarFunction& p,
@@ -676,7 +676,7 @@ AcousticSolverHandler::AcousticSolver<MeshType>::_applyExternalVelocityBC(const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::FixedBoundaryCondition
 {
  private:
@@ -694,7 +694,7 @@ class AcousticSolverHandler::AcousticSolver<MeshType>::FixedBoundaryCondition
   ~FixedBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::PressureBoundaryCondition
 {
  private:
@@ -748,7 +748,7 @@ class AcousticSolverHandler::AcousticSolver<Mesh<1>>::PressureBoundaryCondition
   ~PressureBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::ExternalFSIVelocityBoundaryCondition
 {
  private:
@@ -797,7 +797,7 @@ class AcousticSolverHandler::AcousticSolver<MeshType>::ExternalFSIVelocityBounda
   ~ExternalFSIVelocityBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::VelocityBoundaryCondition
 {
  private:
@@ -826,7 +826,7 @@ class AcousticSolverHandler::AcousticSolver<MeshType>::VelocityBoundaryCondition
   ~VelocityBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::SymmetryBoundaryCondition
 {
  public:
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 75e08981d..73f9522ba 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -1,6 +1,8 @@
 #ifndef ACOUSTIC_SOLVER_HPP
 #define ACOUSTIC_SOLVER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 #include <memory>
 #include <tuple>
 #include <vector>
@@ -59,14 +61,14 @@ class AcousticSolverHandler
           const std::shared_ptr<const DiscreteFunctionVariant>& p,
           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
-    IAcousticSolver()                             = default;
-    IAcousticSolver(IAcousticSolver&&)            = default;
+    IAcousticSolver()                  = default;
+    IAcousticSolver(IAcousticSolver&&) = default;
     IAcousticSolver& operator=(IAcousticSolver&&) = default;
 
     virtual ~IAcousticSolver() = default;
   };
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   class AcousticSolver;
 
   std::unique_ptr<IAcousticSolver> m_acoustic_solver;
diff --git a/src/scheme/CellIntegrator.hpp b/src/scheme/CellIntegrator.hpp
index 26544b7c4..33241f9b8 100644
--- a/src/scheme/CellIntegrator.hpp
+++ b/src/scheme/CellIntegrator.hpp
@@ -75,7 +75,7 @@ class CellIntegrator
     CellList(const ArrayT<CellId>& cell_list) : m_cell_list{cell_list} {}
   };
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _tensorialIntegrateTo(std::function<OutputType(InputType)> function,
                         const IQuadratureDescriptor& quadrature_descriptor,
@@ -303,7 +303,7 @@ class CellIntegrator
     // LCOV_EXCL_STOP
   }
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _directIntegrateTo(std::function<OutputType(InputType)> function,
                      const IQuadratureDescriptor& quadrature_descriptor,
@@ -514,7 +514,7 @@ class CellIntegrator
   }
 
  public:
-  template <typename MeshType, typename OutputArrayT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   integrateTo(const std::function<OutputType(InputType)>& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -532,7 +532,7 @@ class CellIntegrator
     }
   }
 
-  template <typename MeshType, typename OutputArrayT, typename FunctionType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename FunctionType>
   static PUGS_INLINE void
   integrateTo(const FunctionType& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -542,7 +542,7 @@ class CellIntegrator
     integrateTo(std::function(function), quadrature_descriptor, mesh, value);
   }
 
-  template <typename MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const std::function<OutputType(InputType)>& function,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -561,7 +561,7 @@ class CellIntegrator
     return value;
   }
 
-  template <typename MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE auto
   integrate(const FunctionType& function,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/scheme/DiscreteFunctionIntegrator.cpp b/src/scheme/DiscreteFunctionIntegrator.cpp
index 92dd526d2..5caf8e501 100644
--- a/src/scheme/DiscreteFunctionIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionIntegrator.cpp
@@ -8,7 +8,7 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionIntegrator::_integrateOnZoneList() const
 {
@@ -67,7 +67,7 @@ DiscreteFunctionIntegrator::_integrateOnZoneList() const
   return DiscreteFunctionP0<ValueType>(m_mesh, cell_value);
 }
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionIntegrator::_integrateGlobally() const
 {
@@ -82,7 +82,7 @@ DiscreteFunctionIntegrator::_integrateGlobally() const
                                          template integrate<MeshType>(m_function_id, *m_quadrature_descriptor, *mesh));
 }
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionIntegrator::_integrate() const
 {
@@ -93,7 +93,7 @@ DiscreteFunctionIntegrator::_integrate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 DiscreteFunctionVariant
 DiscreteFunctionIntegrator::_integrate() const
 {
diff --git a/src/scheme/DiscreteFunctionIntegrator.hpp b/src/scheme/DiscreteFunctionIntegrator.hpp
index b832ee6f4..f294d39e1 100644
--- a/src/scheme/DiscreteFunctionIntegrator.hpp
+++ b/src/scheme/DiscreteFunctionIntegrator.hpp
@@ -4,6 +4,7 @@
 #include <analysis/IQuadratureDescriptor.hpp>
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 
 #include <memory>
 
@@ -18,16 +19,16 @@ class DiscreteFunctionIntegrator
   std::shared_ptr<const IQuadratureDescriptor> m_quadrature_descriptor;
   const FunctionSymbolId m_function_id;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _integrateOnZoneList() const;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _integrateGlobally() const;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _integrate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionVariant _integrate() const;
 
  public:
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index f8e5a58ca..8f2e0eef7 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -8,7 +8,7 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionInterpoler::_interpolateOnZoneList() const
 {
@@ -68,7 +68,7 @@ DiscreteFunctionInterpoler::_interpolateOnZoneList() const
   return DiscreteFunctionP0<ValueType>{m_mesh, cell_value};
 }
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionInterpoler::_interpolateGlobally() const
 {
@@ -98,7 +98,7 @@ DiscreteFunctionInterpoler::_interpolateGlobally() const
   }
 }
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionInterpoler::_interpolate() const
 {
@@ -109,7 +109,7 @@ DiscreteFunctionInterpoler::_interpolate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 DiscreteFunctionVariant
 DiscreteFunctionInterpoler::_interpolate() const
 {
diff --git a/src/scheme/DiscreteFunctionInterpoler.hpp b/src/scheme/DiscreteFunctionInterpoler.hpp
index 73c3e55f2..95350e34c 100644
--- a/src/scheme/DiscreteFunctionInterpoler.hpp
+++ b/src/scheme/DiscreteFunctionInterpoler.hpp
@@ -3,6 +3,7 @@
 
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
 class DiscreteFunctionVariant;
@@ -18,16 +19,16 @@ class DiscreteFunctionInterpoler
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const FunctionSymbolId m_function_id;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _interpolateOnZoneList() const;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _interpolateGlobally() const;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _interpolate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionVariant _interpolate() const;
 
  public:
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index 56642754b..df43f5290 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -735,12 +735,12 @@ class DiscreteFunctionP0
     Assert(mesh_v->connectivity().id() == cell_value.connectivity_ptr()->id());
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionP0(const std::shared_ptr<const MeshType>& p_mesh)
     : m_mesh_v{p_mesh->meshVariant()}, m_cell_values{m_mesh_v->connectivity()}
   {}
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionP0(const std::shared_ptr<const MeshType>& p_mesh, const CellValue<DataType>& cell_value)
     : m_mesh_v{p_mesh->meshVariant()}, m_cell_values{cell_value}
   {
diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
index aaa549173..da3904b2a 100644
--- a/src/scheme/DiscreteFunctionP0Vector.hpp
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -240,12 +240,12 @@ class DiscreteFunctionP0Vector
     Assert(m_mesh->connectivity().id() == cell_arrays.connectivity_ptr()->id());
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& p_mesh, size_t size)
     : m_mesh{p_mesh->meshVariant()}, m_cell_arrays{m_mesh->connectivity(), size}
   {}
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& p_mesh, const CellArray<DataType>& cell_arrays)
     : m_mesh{p_mesh->meshVariant()}, m_cell_arrays{cell_arrays}
   {
diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.cpp b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
index 74bb43d66..fccc6f79d 100644
--- a/src/scheme/DiscreteFunctionVectorIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
@@ -7,7 +7,7 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorIntegrator::_integrateOnZoneList() const
 {
@@ -68,7 +68,7 @@ DiscreteFunctionVectorIntegrator::_integrateOnZoneList() const
   return DiscreteFunctionP0Vector<DataType>(m_mesh, cell_array);
 }
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorIntegrator::_integrateGlobally() const
 {
@@ -83,7 +83,7 @@ DiscreteFunctionVectorIntegrator::_integrateGlobally() const
                                                                                       *m_quadrature_descriptor, *mesh));
 }
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorIntegrator::_integrate() const
 {
@@ -94,7 +94,7 @@ DiscreteFunctionVectorIntegrator::_integrate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorIntegrator::_integrate() const
 {
diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.hpp b/src/scheme/DiscreteFunctionVectorIntegrator.hpp
index e89943cfb..b9b78525f 100644
--- a/src/scheme/DiscreteFunctionVectorIntegrator.hpp
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.hpp
@@ -4,6 +4,7 @@
 #include <analysis/IQuadratureDescriptor.hpp>
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
 #include <memory>
@@ -21,16 +22,16 @@ class DiscreteFunctionVectorIntegrator
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const std::vector<FunctionSymbolId> m_function_id_list;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _integrateOnZoneList() const;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _integrateGlobally() const;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _integrate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionVariant _integrate() const;
 
  public:
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.cpp b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
index fffead074..661eb483c 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
@@ -7,7 +7,7 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorInterpoler::_interpolateOnZoneList() const
 {
@@ -65,7 +65,7 @@ DiscreteFunctionVectorInterpoler::_interpolateOnZoneList() const
   return DiscreteFunctionP0Vector<DataType>(p_mesh, cell_array);
 }
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorInterpoler::_interpolateGlobally() const
 {
@@ -81,7 +81,7 @@ DiscreteFunctionVectorInterpoler::_interpolateGlobally() const
                                               ItemType::cell>(m_function_id_list, mesh_data.xj()));
 }
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorInterpoler::_interpolate() const
 {
@@ -92,7 +92,7 @@ DiscreteFunctionVectorInterpoler::_interpolate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorInterpoler::_interpolate() const
 {
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.hpp b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
index 2999f1220..ee8536d2d 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.hpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
@@ -3,6 +3,7 @@
 
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
 class DiscreteFunctionVariant;
@@ -19,16 +20,16 @@ class DiscreteFunctionVectorInterpoler
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const std::vector<FunctionSymbolId> m_function_id_list;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _interpolateOnZoneList() const;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _interpolateGlobally() const;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _interpolate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionVariant _interpolate() const;
 
  public:
diff --git a/src/scheme/EdgeIntegrator.hpp b/src/scheme/EdgeIntegrator.hpp
index 2bf287c3d..7557a7ed8 100644
--- a/src/scheme/EdgeIntegrator.hpp
+++ b/src/scheme/EdgeIntegrator.hpp
@@ -67,7 +67,7 @@ class EdgeIntegrator
     EdgeList(const ArrayT<EdgeId>& edge_list) : m_edge_list{edge_list} {}
   };
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _integrateTo(std::function<OutputType(InputType)> function,
                const IQuadratureDescriptor& quadrature_descriptor,
@@ -116,7 +116,7 @@ class EdgeIntegrator
   }
 
  public:
-  template <typename MeshType, typename OutputArrayT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   integrateTo(const std::function<OutputType(InputType)>& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -129,7 +129,7 @@ class EdgeIntegrator
                                          AllEdges<Dimension>{mesh.connectivity()}, value);
   }
 
-  template <typename MeshType, typename OutputArrayT, typename FunctionType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename FunctionType>
   static PUGS_INLINE void
   integrateTo(const FunctionType& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -139,7 +139,7 @@ class EdgeIntegrator
     integrateTo(std::function(function), quadrature_descriptor, mesh, value);
   }
 
-  template <typename MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const std::function<OutputType(InputType)>& function,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -152,7 +152,7 @@ class EdgeIntegrator
     return value;
   }
 
-  template <typename MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE auto
   integrate(const FunctionType& function,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/scheme/FaceIntegrator.hpp b/src/scheme/FaceIntegrator.hpp
index 657ca7a88..3fc6b25f7 100644
--- a/src/scheme/FaceIntegrator.hpp
+++ b/src/scheme/FaceIntegrator.hpp
@@ -69,7 +69,7 @@ class FaceIntegrator
     FaceList(const ArrayT<FaceId>& face_list) : m_face_list{face_list} {}
   };
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _tensorialIntegrateTo(std::function<OutputType(InputType)> function,
                         const IQuadratureDescriptor& quadrature_descriptor,
@@ -167,7 +167,7 @@ class FaceIntegrator
     }
   }
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _directIntegrateTo(std::function<OutputType(InputType)> function,
                      const IQuadratureDescriptor& quadrature_descriptor,
@@ -265,7 +265,7 @@ class FaceIntegrator
   }
 
  public:
-  template <typename MeshType, typename OutputArrayT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   integrateTo(const std::function<OutputType(InputType)>& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -283,7 +283,7 @@ class FaceIntegrator
     }
   }
 
-  template <typename MeshType, typename OutputArrayT, typename FunctionType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename FunctionType>
   static PUGS_INLINE void
   integrateTo(const FunctionType& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -293,7 +293,7 @@ class FaceIntegrator
     integrateTo(std::function(function), quadrature_descriptor, mesh, value);
   }
 
-  template <typename MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const std::function<OutputType(InputType)>& function,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -312,7 +312,7 @@ class FaceIntegrator
     return value;
   }
 
-  template <typename MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE auto
   integrate(const FunctionType& function,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index 509156b1e..a2ebb2b25 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -27,7 +27,7 @@
 #include <variant>
 #include <vector>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class FluxingAdvectionSolver
 {
  private:
@@ -207,7 +207,7 @@ class FluxingAdvectionSolver
   ~FluxingAdvectionSolver() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 FluxingAdvectionSolver<MeshType>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
 {
@@ -385,7 +385,7 @@ FluxingAdvectionSolver<Mesh<3>>::_computeAlgebraicFluxingVolume()
   return algebraic_fluxing_volume;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 FaceValue<double>
 FluxingAdvectionSolver<MeshType>::_computeFluxingVolume(FaceValue<double> algebraic_fluxing_volumes)
 {
@@ -400,7 +400,7 @@ FluxingAdvectionSolver<MeshType>::_computeFluxingVolume(FaceValue<double> algebr
   return algebraic_fluxing_volumes;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 FluxingAdvectionSolver<MeshType>::_computeCycleNumber(FaceValue<double> fluxing_volumes)
 {
@@ -443,7 +443,7 @@ FluxingAdvectionSolver<MeshType>::_computeCycleNumber(FaceValue<double> fluxing_
   m_cycle_fluxing_volume = fluxing_volumes;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 FluxingAdvectionSolver<MeshType>::_computeGeometricalData()
 {
@@ -453,7 +453,7 @@ FluxingAdvectionSolver<MeshType>::_computeGeometricalData()
   this->_computeCycleNumber(fluxing_volumes);
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 template <typename CellDataType>
 void
 FluxingAdvectionSolver<MeshType>::_remapOne(const CellValue<const double>& step_Vj,
@@ -633,7 +633,7 @@ FluxingAdvectionSolver<MeshType>::_remapOne(const CellValue<const double>& step_
   old_q = new_q;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 FluxingAdvectionSolver<MeshType>::_remapAllQuantities()
 {
@@ -696,7 +696,7 @@ FluxingAdvectionSolver<MeshType>::_remapAllQuantities()
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
 FluxingAdvectionSolver<MeshType>::remap(
   const std::vector<std::shared_ptr<const VariableBCDescriptor>>& quantity_list_with_bc)
@@ -748,7 +748,7 @@ FluxingAdvectionSolver<MeshType>::remap(
   return new_variables;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 template <typename DataType>
 class FluxingAdvectionSolver<MeshType>::InflowValueBoundaryCondition
 {
@@ -779,7 +779,7 @@ class FluxingAdvectionSolver<MeshType>::InflowValueBoundaryCondition
   ~InflowValueBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class FluxingAdvectionSolver<MeshType>::InflowArrayBoundaryCondition
 {
  private:
@@ -809,7 +809,7 @@ class FluxingAdvectionSolver<MeshType>::InflowArrayBoundaryCondition
   ~InflowArrayBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class FluxingAdvectionSolver<MeshType>::OutflowBoundaryCondition
 {
  private:
@@ -830,7 +830,7 @@ class FluxingAdvectionSolver<MeshType>::OutflowBoundaryCondition
   ~OutflowBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class FluxingAdvectionSolver<MeshType>::SymmetryBoundaryCondition
 {
  private:
diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index 5309385ec..d88067f98 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -48,7 +48,7 @@ hyperelastic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
     c.meshVariant()->variant());
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticSolverHandler::IHyperelasticSolver
 {
  private:
@@ -549,7 +549,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
   ~HyperelasticSolver()                    = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
                                                                           const MeshType& mesh,
@@ -614,7 +614,7 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyPressureBC(const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyNormalStressBC(const BoundaryConditionList& bc_list,
                                                                               const MeshType& mesh,
@@ -679,7 +679,7 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyNormalStressBC(co
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
                                                                           NodeValue<Rdxd>& Ar,
@@ -710,7 +710,7 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applySymmetryBC(const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyVelocityBC(const BoundaryConditionList& bc_list,
                                                                           NodeValue<Rdxd>& Ar,
@@ -747,7 +747,7 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyVelocityBC(const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::FixedBoundaryCondition
 {
  private:
@@ -765,7 +765,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::FixedBoundaryCond
   ~FixedBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::PressureBoundaryCondition
 {
  private:
@@ -819,7 +819,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Mesh<1>>::PressureBoundaryCo
   ~PressureBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::NormalStressBoundaryCondition
 {
  private:
@@ -873,7 +873,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Mesh<1>>::NormalStressBounda
   ~NormalStressBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::VelocityBoundaryCondition
 {
  private:
@@ -902,7 +902,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::VelocityBoundaryC
   ~VelocityBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::SymmetryBoundaryCondition
 {
  public:
diff --git a/src/scheme/HyperelasticSolver.hpp b/src/scheme/HyperelasticSolver.hpp
index 75963a519..0fe58e857 100644
--- a/src/scheme/HyperelasticSolver.hpp
+++ b/src/scheme/HyperelasticSolver.hpp
@@ -1,6 +1,8 @@
 #ifndef HYPERELASTIC_SOLVER_HPP
 #define HYPERELASTIC_SOLVER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 #include <memory>
 #include <tuple>
 #include <vector>
@@ -65,14 +67,14 @@ class HyperelasticSolverHandler
           const std::shared_ptr<const DiscreteFunctionVariant>& p,
           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
-    IHyperelasticSolver()                                 = default;
-    IHyperelasticSolver(IHyperelasticSolver&&)            = default;
+    IHyperelasticSolver()                      = default;
+    IHyperelasticSolver(IHyperelasticSolver&&) = default;
     IHyperelasticSolver& operator=(IHyperelasticSolver&&) = default;
 
     virtual ~IHyperelasticSolver() = default;
   };
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   class HyperelasticSolver;
 
   std::unique_ptr<IHyperelasticSolver> m_hyperelastic_solver;
-- 
GitLab