From 83670ac0666617809288dc54d66a9046eebb9e4b Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 5 Mar 2024 00:58:22 +0100
Subject: [PATCH] Simplify template code

-add MeshConcept to simplify compilation error
-add mesh_type_t template utility to simplify
-simplify code (dealing with constness of mesh types was not that easy)
-check constness (at compile time) to help code writing in boundary utilities
---
 src/language/modules/SchemeModule.cpp         |  2 +-
 src/language/utils/IntegrateCellValue.hpp     |  6 ++---
 .../ItemArrayVariantFunctionInterpoler.cpp    |  2 +-
 .../ItemValueVariantFunctionInterpoler.cpp    |  2 +-
 src/mesh/DiamondDualMeshBuilder.cpp           |  2 +-
 src/mesh/DualMeshManager.cpp                  |  5 ++--
 src/mesh/MedianDualMeshBuilder.cpp            |  2 +-
 src/mesh/MeshFlatEdgeBoundary.cpp             | 12 ++++-----
 src/mesh/MeshFlatEdgeBoundary.hpp             | 10 ++++---
 src/mesh/MeshFlatFaceBoundary.cpp             | 12 ++++-----
 src/mesh/MeshFlatFaceBoundary.hpp             | 10 ++++---
 src/mesh/MeshFlatNodeBoundary.cpp             | 26 +++++++++----------
 src/mesh/MeshFlatNodeBoundary.hpp             | 10 ++++---
 src/mesh/MeshLineEdgeBoundary.cpp             |  9 +++----
 src/mesh/MeshLineEdgeBoundary.hpp             | 10 ++++---
 src/mesh/MeshLineFaceBoundary.cpp             |  7 +++--
 src/mesh/MeshLineFaceBoundary.hpp             |  9 ++++---
 src/mesh/MeshLineNodeBoundary.cpp             | 16 ++++++------
 src/mesh/MeshLineNodeBoundary.hpp             |  9 ++++---
 src/mesh/MeshRandomizer.cpp                   | 12 ++++-----
 src/mesh/MeshRelaxer.cpp                      |  4 +--
 src/mesh/MeshSmoother.cpp                     |  8 +++---
 src/mesh/MeshTraits.hpp                       | 23 +++++++++++++++-
 src/mesh/MeshTransformer.cpp                  |  3 ++-
 src/output/GnuplotWriter.cpp                  |  7 ++---
 src/output/GnuplotWriter1D.cpp                |  5 ++--
 src/scheme/AcousticSolver.cpp                 |  4 +--
 src/scheme/DiscreteFunctionIntegrator.cpp     |  5 ++--
 src/scheme/DiscreteFunctionInterpoler.cpp     |  4 ++-
 .../DiscreteFunctionVectorIntegrator.cpp      |  4 ++-
 .../DiscreteFunctionVectorInterpoler.cpp      |  4 ++-
 src/scheme/FluxingAdvectionSolver.cpp         | 13 +++++-----
 src/scheme/HyperelasticSolver.cpp             |  6 ++---
 33 files changed, 152 insertions(+), 111 deletions(-)

diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index e8d4bb390..7b404fcf6 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -637,7 +637,7 @@ SchemeModule::SchemeModule()
                                 -> std::shared_ptr<const DiscreteFunctionVariant> {
                                 return std::visit(
                                   [&](auto&& mesh) {
-                                    using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+                                    using MeshType = mesh_type_t<decltype(mesh)>;
                                     if constexpr (is_polygonal_mesh_v<MeshType>) {
                                       return std::make_shared<DiscreteFunctionVariant>(
                                         DiscreteFunctionP0(mesh_v,
diff --git a/src/language/utils/IntegrateCellValue.hpp b/src/language/utils/IntegrateCellValue.hpp
index ecd304ca7..cde3bfacc 100644
--- a/src/language/utils/IntegrateCellValue.hpp
+++ b/src/language/utils/IntegrateCellValue.hpp
@@ -36,7 +36,7 @@ class IntegrateCellValue<OutputType(InputType)>
   {
     return std::visit(
       [&](auto&& p_mesh) -> CellValue<OutputType> {
-        using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+        using MeshType = mesh_type_t<decltype(p_mesh)>;
         if constexpr ((is_polygonal_mesh_v<MeshType>)and(MeshType::Dimension == InputType::Dimension)) {
           return integrate(function_symbol_id, quadrature_descriptor, *p_mesh);
         } else {
@@ -66,7 +66,7 @@ class IntegrateCellValue<OutputType(InputType)>
   {
     return std::visit(
       [&](auto&& p_mesh) -> Array<OutputType> {
-        using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+        using MeshType = mesh_type_t<decltype(p_mesh)>;
         if constexpr ((is_polygonal_mesh_v<MeshType>)and(MeshType::Dimension == InputType::Dimension)) {
           return integrate(function_symbol_id, quadrature_descriptor, *p_mesh, list_of_cells);
         } else {
@@ -95,7 +95,7 @@ class IntegrateCellValue<OutputType(InputType)>
   {
     return std::visit(
       [&](auto&& p_mesh) -> Array<OutputType> {
-        using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+        using MeshType = mesh_type_t<decltype(p_mesh)>;
         if constexpr ((is_polygonal_mesh_v<MeshType>)and(MeshType::Dimension == InputType::Dimension)) {
           return integrate(function_symbol_id, quadrature_descriptor, *p_mesh, list_of_cells);
         } else {
diff --git a/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
index a288765cc..1c4b2b9f0 100644
--- a/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
+++ b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
@@ -147,7 +147,7 @@ ItemArrayVariantFunctionInterpoler::interpolate() const
 {
   return std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(mesh)>;
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         return this->_interpolate<MeshType>();
       } else {
diff --git a/src/language/utils/ItemValueVariantFunctionInterpoler.cpp b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
index 1273a8da4..74d07b577 100644
--- a/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
+++ b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
@@ -133,7 +133,7 @@ ItemValueVariantFunctionInterpoler::interpolate() const
 {
   return std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(mesh)>;
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         return this->_interpolate<MeshType>();
       } else {
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 36cb042d3..f6ecba946 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -45,7 +45,7 @@ DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const MeshV
 {
   std::visit(
     [&](auto&& p_mesh) {
-      using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         if constexpr (MeshType::Dimension > 1) {
           this->_buildDualDiamondMeshFrom(*p_mesh);
diff --git a/src/mesh/DualMeshManager.cpp b/src/mesh/DualMeshManager.cpp
index 846033489..a09413776 100644
--- a/src/mesh/DualMeshManager.cpp
+++ b/src/mesh/DualMeshManager.cpp
@@ -5,6 +5,7 @@
 #include <mesh/Dual1DMeshBuilder.hpp>
 #include <mesh/MedianDualMeshBuilder.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
 #include <utils/Exceptions.hpp>
 #include <utils/PugsAssert.hpp>
@@ -77,7 +78,7 @@ DualMeshManager::getMedianDualMesh(const std::shared_ptr<const MeshVariant>& mes
 {
   return std::visit(
     [&](auto&& p_mesh) {
-      using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
       if constexpr (MeshType::Dimension == 1) {
         return this->getDual1DMesh(mesh_v);
       } else {
@@ -100,7 +101,7 @@ DualMeshManager::getDiamondDualMesh(const std::shared_ptr<const MeshVariant>& me
 {
   return std::visit(
     [&](auto&& p_mesh) {
-      using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
       if constexpr (MeshType::Dimension == 1) {
         return this->getDual1DMesh(mesh_v);
       } else {
diff --git a/src/mesh/MedianDualMeshBuilder.cpp b/src/mesh/MedianDualMeshBuilder.cpp
index a895dd088..5a22cf0b2 100644
--- a/src/mesh/MedianDualMeshBuilder.cpp
+++ b/src/mesh/MedianDualMeshBuilder.cpp
@@ -46,7 +46,7 @@ MedianDualMeshBuilder::MedianDualMeshBuilder(const std::shared_ptr<const MeshVar
 {
   std::visit(
     [&](auto&& p_mesh) {
-      using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         if constexpr (MeshType::Dimension > 1) {
           this->_buildMedianDualMeshFrom(*p_mesh);
diff --git a/src/mesh/MeshFlatEdgeBoundary.cpp b/src/mesh/MeshFlatEdgeBoundary.cpp
index 8b5669a34..7a2d4918d 100644
--- a/src/mesh/MeshFlatEdgeBoundary.cpp
+++ b/src/mesh/MeshFlatEdgeBoundary.cpp
@@ -5,16 +5,16 @@
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
 template <typename MeshType>
-MeshFlatEdgeBoundary<const MeshType>
+MeshFlatEdgeBoundary<MeshType>
 getMeshFlatEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   MeshEdgeBoundary mesh_edge_boundary          = getMeshEdgeBoundary(mesh, boundary_descriptor);
   MeshFlatNodeBoundary mesh_flat_node_boundary = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshFlatEdgeBoundary<const MeshType>{mesh, mesh_edge_boundary.refEdgeList(),
-                                              mesh_flat_node_boundary.outgoingNormal()};
+  return MeshFlatEdgeBoundary<MeshType>{mesh, mesh_edge_boundary.refEdgeList(),
+                                        mesh_flat_node_boundary.outgoingNormal()};
 }
 
-template MeshFlatEdgeBoundary<const Mesh<1>> getMeshFlatEdgeBoundary(const Mesh<1>&, const IBoundaryDescriptor&);
-template MeshFlatEdgeBoundary<const Mesh<2>> getMeshFlatEdgeBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
-template MeshFlatEdgeBoundary<const Mesh<3>> getMeshFlatEdgeBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
+template MeshFlatEdgeBoundary<Mesh<1>> getMeshFlatEdgeBoundary(const Mesh<1>&, const IBoundaryDescriptor&);
+template MeshFlatEdgeBoundary<Mesh<2>> getMeshFlatEdgeBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
+template MeshFlatEdgeBoundary<Mesh<3>> getMeshFlatEdgeBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFlatEdgeBoundary.hpp b/src/mesh/MeshFlatEdgeBoundary.hpp
index 60b54257f..f395c4ebf 100644
--- a/src/mesh/MeshFlatEdgeBoundary.hpp
+++ b/src/mesh/MeshFlatEdgeBoundary.hpp
@@ -7,6 +7,8 @@ template <typename MeshType>
 class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
+  static_assert(not std::is_const_v<MeshType>, "MeshType must be non-const");
+
   using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
@@ -23,8 +25,8 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
   MeshFlatEdgeBoundary& operator=(MeshFlatEdgeBoundary&&)      = default;
 
   template <typename MeshTypeT>
-  friend MeshFlatEdgeBoundary<const MeshTypeT> getMeshFlatEdgeBoundary(const MeshTypeT& mesh,
-                                                                       const IBoundaryDescriptor& boundary_descriptor);
+  friend MeshFlatEdgeBoundary<MeshTypeT> getMeshFlatEdgeBoundary(const MeshTypeT& mesh,
+                                                                 const IBoundaryDescriptor& boundary_descriptor);
 
  private:
   MeshFlatEdgeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list, const Rd& outgoing_normal)
@@ -39,7 +41,7 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
 };
 
 template <typename MeshType>
-MeshFlatEdgeBoundary<const MeshType> getMeshFlatEdgeBoundary(const MeshType& mesh,
-                                                             const IBoundaryDescriptor& boundary_descriptor);
+MeshFlatEdgeBoundary<MeshType> getMeshFlatEdgeBoundary(const MeshType& mesh,
+                                                       const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_FLAT_EDGE_BOUNDARY_HPP
diff --git a/src/mesh/MeshFlatFaceBoundary.cpp b/src/mesh/MeshFlatFaceBoundary.cpp
index 12d7f439f..2ecdf6768 100644
--- a/src/mesh/MeshFlatFaceBoundary.cpp
+++ b/src/mesh/MeshFlatFaceBoundary.cpp
@@ -5,16 +5,16 @@
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
 template <typename MeshType>
-MeshFlatFaceBoundary<const MeshType>
+MeshFlatFaceBoundary<MeshType>
 getMeshFlatFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   MeshFaceBoundary mesh_face_boundary          = getMeshFaceBoundary(mesh, boundary_descriptor);
   MeshFlatNodeBoundary mesh_flat_node_boundary = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshFlatFaceBoundary<const MeshType>{mesh, mesh_face_boundary.refFaceList(),
-                                              mesh_flat_node_boundary.outgoingNormal()};
+  return MeshFlatFaceBoundary<MeshType>{mesh, mesh_face_boundary.refFaceList(),
+                                        mesh_flat_node_boundary.outgoingNormal()};
 }
 
-template MeshFlatFaceBoundary<const Mesh<1>> getMeshFlatFaceBoundary(const Mesh<1>&, const IBoundaryDescriptor&);
-template MeshFlatFaceBoundary<const Mesh<2>> getMeshFlatFaceBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
-template MeshFlatFaceBoundary<const Mesh<3>> getMeshFlatFaceBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
+template MeshFlatFaceBoundary<Mesh<1>> getMeshFlatFaceBoundary(const Mesh<1>&, const IBoundaryDescriptor&);
+template MeshFlatFaceBoundary<Mesh<2>> getMeshFlatFaceBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
+template MeshFlatFaceBoundary<Mesh<3>> getMeshFlatFaceBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFlatFaceBoundary.hpp b/src/mesh/MeshFlatFaceBoundary.hpp
index 3620d0f3d..1aa2e854f 100644
--- a/src/mesh/MeshFlatFaceBoundary.hpp
+++ b/src/mesh/MeshFlatFaceBoundary.hpp
@@ -7,6 +7,8 @@ template <typename MeshType>
 class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
+  static_assert(not std::is_const_v<MeshType>, "MeshType must be non-const");
+
   using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
@@ -23,8 +25,8 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
   MeshFlatFaceBoundary& operator=(MeshFlatFaceBoundary&&)      = default;
 
   template <typename MeshTypeT>
-  friend MeshFlatFaceBoundary<const MeshTypeT> getMeshFlatFaceBoundary(const MeshTypeT& mesh,
-                                                                       const IBoundaryDescriptor& boundary_descriptor);
+  friend MeshFlatFaceBoundary<MeshTypeT> getMeshFlatFaceBoundary(const MeshTypeT& mesh,
+                                                                 const IBoundaryDescriptor& boundary_descriptor);
 
  private:
   MeshFlatFaceBoundary(const MeshType& mesh, const RefFaceList& ref_face_list, const Rd& outgoing_normal)
@@ -39,7 +41,7 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
 };
 
 template <typename MeshType>
-MeshFlatFaceBoundary<const MeshType> getMeshFlatFaceBoundary(const MeshType& mesh,
-                                                             const IBoundaryDescriptor& boundary_descriptor);
+MeshFlatFaceBoundary<MeshType> getMeshFlatFaceBoundary(const MeshType& mesh,
+                                                       const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_FLAT_FACE_BOUNDARY_HPP
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index df7631988..983f96738 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -34,7 +34,7 @@ MeshFlatNodeBoundary<MeshType>::_checkBoundaryIsFlat(const TinyVector<MeshType::
 
 template <>
 TinyVector<1, double>
-MeshFlatNodeBoundary<const Mesh<1>>::_getNormal(const Mesh<1>& mesh)
+MeshFlatNodeBoundary<Mesh<1>>::_getNormal(const Mesh<1>& mesh)
 {
   using R = TinyVector<1, double>;
 
@@ -60,7 +60,7 @@ MeshFlatNodeBoundary<const Mesh<1>>::_getNormal(const Mesh<1>& mesh)
 
 template <>
 TinyVector<2, double>
-MeshFlatNodeBoundary<const Mesh<2>>::_getNormal(const Mesh<2>& mesh)
+MeshFlatNodeBoundary<Mesh<2>>::_getNormal(const Mesh<2>& mesh)
 {
   using R2 = TinyVector<2, double>;
 
@@ -88,7 +88,7 @@ MeshFlatNodeBoundary<const Mesh<2>>::_getNormal(const Mesh<2>& mesh)
 
 template <>
 TinyVector<3, double>
-MeshFlatNodeBoundary<const Mesh<3>>::_getFarestNode(const Mesh<3>& mesh, const Rd& x0, const Rd& x1)
+MeshFlatNodeBoundary<Mesh<3>>::_getFarestNode(const Mesh<3>& mesh, const Rd& x0, const Rd& x1)
 {
   const NodeValue<const Rd>& xr = mesh.xr();
   const auto node_number        = mesh.connectivity().nodeNumber();
@@ -139,7 +139,7 @@ MeshFlatNodeBoundary<const Mesh<3>>::_getFarestNode(const Mesh<3>& mesh, const R
 
 template <>
 TinyVector<3, double>
-MeshFlatNodeBoundary<const Mesh<3>>::_getNormal(const Mesh<3>& mesh)
+MeshFlatNodeBoundary<Mesh<3>>::_getNormal(const Mesh<3>& mesh)
 {
   using R3 = TinyVector<3, double>;
 
@@ -196,7 +196,7 @@ MeshFlatNodeBoundary<const Mesh<3>>::_getNormal(const Mesh<3>& mesh)
 
 template <>
 TinyVector<1, double>
-MeshFlatNodeBoundary<const Mesh<1>>::_getOutgoingNormal(const Mesh<1>& mesh)
+MeshFlatNodeBoundary<Mesh<1>>::_getOutgoingNormal(const Mesh<1>& mesh)
 {
   using R = TinyVector<1, double>;
 
@@ -239,7 +239,7 @@ MeshFlatNodeBoundary<const Mesh<1>>::_getOutgoingNormal(const Mesh<1>& mesh)
 
 template <>
 TinyVector<2, double>
-MeshFlatNodeBoundary<const Mesh<2>>::_getOutgoingNormal(const Mesh<2>& mesh)
+MeshFlatNodeBoundary<Mesh<2>>::_getOutgoingNormal(const Mesh<2>& mesh)
 {
   using R2 = TinyVector<2, double>;
 
@@ -282,7 +282,7 @@ MeshFlatNodeBoundary<const Mesh<2>>::_getOutgoingNormal(const Mesh<2>& mesh)
 
 template <>
 TinyVector<3, double>
-MeshFlatNodeBoundary<const Mesh<3>>::_getOutgoingNormal(const Mesh<3>& mesh)
+MeshFlatNodeBoundary<Mesh<3>>::_getOutgoingNormal(const Mesh<3>& mesh)
 {
   using R3 = TinyVector<3, double>;
 
@@ -324,7 +324,7 @@ MeshFlatNodeBoundary<const Mesh<3>>::_getOutgoingNormal(const Mesh<3>& mesh)
 }
 
 template <typename MeshType>
-MeshFlatNodeBoundary<const MeshType>
+MeshFlatNodeBoundary<MeshType>
 getMeshFlatNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
@@ -332,7 +332,7 @@ getMeshFlatNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundar
     const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
     const RefId& ref          = ref_node_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshFlatNodeBoundary<const MeshType>{mesh, ref_node_list};
+      return MeshFlatNodeBoundary<MeshType>{mesh, ref_node_list};
     }
   }
   for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
@@ -340,7 +340,7 @@ getMeshFlatNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundar
     const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
     const RefId& ref          = ref_face_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshFlatNodeBoundary<const MeshType>{mesh, ref_face_list};
+      return MeshFlatNodeBoundary<MeshType>{mesh, ref_face_list};
     }
   }
 
@@ -350,6 +350,6 @@ getMeshFlatNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundar
   throw NormalError(ost.str());
 }
 
-template MeshFlatNodeBoundary<const Mesh<1>> getMeshFlatNodeBoundary(const Mesh<1>&, const IBoundaryDescriptor&);
-template MeshFlatNodeBoundary<const Mesh<2>> getMeshFlatNodeBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
-template MeshFlatNodeBoundary<const Mesh<3>> getMeshFlatNodeBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<Mesh<1>> getMeshFlatNodeBoundary(const Mesh<1>&, const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<Mesh<2>> getMeshFlatNodeBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<Mesh<3>> getMeshFlatNodeBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFlatNodeBoundary.hpp b/src/mesh/MeshFlatNodeBoundary.hpp
index df3c9d288..453804724 100644
--- a/src/mesh/MeshFlatNodeBoundary.hpp
+++ b/src/mesh/MeshFlatNodeBoundary.hpp
@@ -7,6 +7,8 @@ template <typename MeshType>
 class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
+  static_assert(not std::is_const_v<MeshType>, "MeshType must be non-const");
+
   using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
@@ -34,8 +36,8 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
   MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&)      = default;
 
   template <typename MeshTypeT>
-  friend MeshFlatNodeBoundary<const MeshTypeT> getMeshFlatNodeBoundary(const MeshTypeT& mesh,
-                                                                       const IBoundaryDescriptor& boundary_descriptor);
+  friend MeshFlatNodeBoundary<MeshTypeT> getMeshFlatNodeBoundary(const MeshTypeT& mesh,
+                                                                 const IBoundaryDescriptor& boundary_descriptor);
 
  private:
   MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
@@ -54,7 +56,7 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
 };
 
 template <typename MeshType>
-MeshFlatNodeBoundary<const MeshType> getMeshFlatNodeBoundary(const MeshType& mesh,
-                                                             const IBoundaryDescriptor& boundary_descriptor);
+MeshFlatNodeBoundary<MeshType> getMeshFlatNodeBoundary(const MeshType& mesh,
+                                                       const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_FLAT_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshLineEdgeBoundary.cpp b/src/mesh/MeshLineEdgeBoundary.cpp
index 1d0c37526..f1e0333a9 100644
--- a/src/mesh/MeshLineEdgeBoundary.cpp
+++ b/src/mesh/MeshLineEdgeBoundary.cpp
@@ -6,15 +6,14 @@
 #include <utils/Messenger.hpp>
 
 template <typename MeshType>
-MeshLineEdgeBoundary<const MeshType>
+MeshLineEdgeBoundary<MeshType>
 getMeshLineEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   MeshEdgeBoundary mesh_edge_boundary          = getMeshEdgeBoundary(mesh, boundary_descriptor);
   MeshLineNodeBoundary mesh_line_node_boundary = getMeshLineNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshLineEdgeBoundary<const MeshType>{mesh, mesh_edge_boundary.refEdgeList(),
-                                              mesh_line_node_boundary.direction()};
+  return MeshLineEdgeBoundary<MeshType>{mesh, mesh_edge_boundary.refEdgeList(), mesh_line_node_boundary.direction()};
 }
 
-template MeshLineEdgeBoundary<const Mesh<2>> getMeshLineEdgeBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
-template MeshLineEdgeBoundary<const Mesh<3>> getMeshLineEdgeBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
+template MeshLineEdgeBoundary<Mesh<2>> getMeshLineEdgeBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
+template MeshLineEdgeBoundary<Mesh<3>> getMeshLineEdgeBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshLineEdgeBoundary.hpp b/src/mesh/MeshLineEdgeBoundary.hpp
index b120a6882..39c97a194 100644
--- a/src/mesh/MeshLineEdgeBoundary.hpp
+++ b/src/mesh/MeshLineEdgeBoundary.hpp
@@ -8,6 +8,8 @@ template <typename MeshType>
 class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
+  static_assert(not std::is_const_v<MeshType>, "MeshType must be non-const");
+
   static_assert(MeshType::Dimension > 1, "MeshLineEdgeBoundary makes only sense in dimension 1");
 
   using Rd = TinyVector<MeshType::Dimension, double>;
@@ -17,8 +19,8 @@ class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // cl
 
  public:
   template <typename MeshTypeT>
-  friend MeshLineEdgeBoundary<const MeshTypeT> getMeshLineEdgeBoundary(const MeshTypeT& mesh,
-                                                                       const IBoundaryDescriptor& boundary_descriptor);
+  friend MeshLineEdgeBoundary<MeshTypeT> getMeshLineEdgeBoundary(const MeshTypeT& mesh,
+                                                                 const IBoundaryDescriptor& boundary_descriptor);
 
   PUGS_INLINE
   const Rd&
@@ -43,7 +45,7 @@ class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // cl
 };
 
 template <typename MeshType>
-MeshLineEdgeBoundary<const MeshType> getMeshLineEdgeBoundary(const MeshType& mesh,
-                                                             const IBoundaryDescriptor& boundary_descriptor);
+MeshLineEdgeBoundary<MeshType> getMeshLineEdgeBoundary(const MeshType& mesh,
+                                                       const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_LINE_EDGE_BOUNDARY_HPP
diff --git a/src/mesh/MeshLineFaceBoundary.cpp b/src/mesh/MeshLineFaceBoundary.cpp
index 83bfdf8cf..81ea991ed 100644
--- a/src/mesh/MeshLineFaceBoundary.cpp
+++ b/src/mesh/MeshLineFaceBoundary.cpp
@@ -6,14 +6,13 @@
 #include <utils/Messenger.hpp>
 
 template <typename MeshType>
-MeshLineFaceBoundary<const MeshType>
+MeshLineFaceBoundary<MeshType>
 getMeshLineFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   MeshFaceBoundary mesh_face_boundary          = getMeshFaceBoundary(mesh, boundary_descriptor);
   MeshLineNodeBoundary mesh_line_node_boundary = getMeshLineNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshLineFaceBoundary<const MeshType>{mesh, mesh_face_boundary.refFaceList(),
-                                              mesh_line_node_boundary.direction()};
+  return MeshLineFaceBoundary<MeshType>{mesh, mesh_face_boundary.refFaceList(), mesh_line_node_boundary.direction()};
 }
 
-template MeshLineFaceBoundary<const Mesh<2>> getMeshLineFaceBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
+template MeshLineFaceBoundary<Mesh<2>> getMeshLineFaceBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshLineFaceBoundary.hpp b/src/mesh/MeshLineFaceBoundary.hpp
index f95611bb8..203faac8e 100644
--- a/src/mesh/MeshLineFaceBoundary.hpp
+++ b/src/mesh/MeshLineFaceBoundary.hpp
@@ -8,6 +8,7 @@ template <typename MeshType>
 class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
+  static_assert(not std::is_const_v<MeshType>, "MeshType must be non-const");
   static_assert(MeshType::Dimension == 2, "MeshLineFaceBoundary makes only sense in dimension 2");
 
   using Rd = TinyVector<MeshType::Dimension, double>;
@@ -17,8 +18,8 @@ class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // cl
 
  public:
   template <typename MeshTypeT>
-  friend MeshLineFaceBoundary<const MeshTypeT> getMeshLineFaceBoundary(const MeshTypeT& mesh,
-                                                                       const IBoundaryDescriptor& boundary_descriptor);
+  friend MeshLineFaceBoundary<MeshTypeT> getMeshLineFaceBoundary(const MeshTypeT& mesh,
+                                                                 const IBoundaryDescriptor& boundary_descriptor);
 
   PUGS_INLINE
   const Rd&
@@ -43,7 +44,7 @@ class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // cl
 };
 
 template <typename MeshType>
-MeshLineFaceBoundary<const MeshType> getMeshLineFaceBoundary(const MeshType& mesh,
-                                                             const IBoundaryDescriptor& boundary_descriptor);
+MeshLineFaceBoundary<MeshType> getMeshLineFaceBoundary(const MeshType& mesh,
+                                                       const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_LINE_FACE_BOUNDARY_HPP
diff --git a/src/mesh/MeshLineNodeBoundary.cpp b/src/mesh/MeshLineNodeBoundary.cpp
index dd145d889..1228fa5aa 100644
--- a/src/mesh/MeshLineNodeBoundary.cpp
+++ b/src/mesh/MeshLineNodeBoundary.cpp
@@ -38,7 +38,7 @@ MeshLineNodeBoundary<MeshType>::_checkBoundaryIsLine(const TinyVector<MeshType::
 
 template <>
 TinyVector<2>
-MeshLineNodeBoundary<const Mesh<2>>::_getDirection(const Mesh<2>& mesh)
+MeshLineNodeBoundary<Mesh<2>>::_getDirection(const Mesh<2>& mesh)
 {
   using R2 = TinyVector<2, double>;
 
@@ -65,7 +65,7 @@ MeshLineNodeBoundary<const Mesh<2>>::_getDirection(const Mesh<2>& mesh)
 
 template <>
 TinyVector<3>
-MeshLineNodeBoundary<const Mesh<3>>::_getDirection(const Mesh<3>& mesh)
+MeshLineNodeBoundary<Mesh<3>>::_getDirection(const Mesh<3>& mesh)
 {
   using R3 = TinyVector<3, double>;
 
@@ -106,7 +106,7 @@ MeshLineNodeBoundary<const Mesh<3>>::_getDirection(const Mesh<3>& mesh)
 }
 
 template <typename MeshType>
-MeshLineNodeBoundary<const MeshType>
+MeshLineNodeBoundary<MeshType>
 getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
@@ -114,7 +114,7 @@ getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundar
     const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
     const RefId& ref          = ref_node_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshLineNodeBoundary<const MeshType>{mesh, ref_node_list};
+      return MeshLineNodeBoundary<MeshType>{mesh, ref_node_list};
     }
   }
   for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
@@ -122,7 +122,7 @@ getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundar
     const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
     const RefId& ref          = ref_edge_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshLineNodeBoundary<const MeshType>{mesh, ref_edge_list};
+      return MeshLineNodeBoundary<MeshType>{mesh, ref_edge_list};
     }
   }
   for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
@@ -130,7 +130,7 @@ getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundar
     const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
     const RefId& ref          = ref_face_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshLineNodeBoundary<const MeshType>{mesh, ref_face_list};
+      return MeshLineNodeBoundary<MeshType>{mesh, ref_face_list};
     }
   }
 
@@ -140,5 +140,5 @@ getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundar
   throw NormalError(ost.str());
 }
 
-template MeshLineNodeBoundary<const Mesh<2>> getMeshLineNodeBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
-template MeshLineNodeBoundary<const Mesh<3>> getMeshLineNodeBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
+template MeshLineNodeBoundary<Mesh<2>> getMeshLineNodeBoundary(const Mesh<2>&, const IBoundaryDescriptor&);
+template MeshLineNodeBoundary<Mesh<3>> getMeshLineNodeBoundary(const Mesh<3>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshLineNodeBoundary.hpp b/src/mesh/MeshLineNodeBoundary.hpp
index 50dbfa00b..bb6fc4b04 100644
--- a/src/mesh/MeshLineNodeBoundary.hpp
+++ b/src/mesh/MeshLineNodeBoundary.hpp
@@ -8,6 +8,7 @@ template <typename MeshType>
 class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
+  static_assert(not std::is_const_v<MeshType>, "MeshType must be non-const");
   static_assert(MeshType::Dimension > 1, "MeshLineNodeBoundary makes only sense in dimension greater than 1");
 
   using Rd = TinyVector<MeshType::Dimension, double>;
@@ -24,8 +25,8 @@ class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // cl
 
  public:
   template <typename MeshTypeT>
-  friend MeshLineNodeBoundary<const MeshTypeT> getMeshLineNodeBoundary(const MeshTypeT& mesh,
-                                                                       const IBoundaryDescriptor& boundary_descriptor);
+  friend MeshLineNodeBoundary<MeshTypeT> getMeshLineNodeBoundary(const MeshTypeT& mesh,
+                                                                 const IBoundaryDescriptor& boundary_descriptor);
 
   PUGS_INLINE
   const Rd&
@@ -58,7 +59,7 @@ class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // cl
 };
 
 template <typename MeshType>
-MeshLineNodeBoundary<const MeshType> getMeshLineNodeBoundary(const MeshType& mesh,
-                                                             const IBoundaryDescriptor& boundary_descriptor);
+MeshLineNodeBoundary<MeshType> getMeshLineNodeBoundary(const MeshType& mesh,
+                                                       const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_LINE_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshRandomizer.cpp b/src/mesh/MeshRandomizer.cpp
index 6b7d6f4e4..d6b8c2f1b 100644
--- a/src/mesh/MeshRandomizer.cpp
+++ b/src/mesh/MeshRandomizer.cpp
@@ -281,7 +281,7 @@ class MeshRandomizerHandler::MeshRandomizer<MeshType>::AxisBoundaryCondition
   using Rd = TinyVector<Dimension, double>;
 
  private:
-  const MeshLineNodeBoundary<const MeshType> m_mesh_line_node_boundary;
+  const MeshLineNodeBoundary<MeshType> m_mesh_line_node_boundary;
 
  public:
   const Rd&
@@ -296,7 +296,7 @@ class MeshRandomizerHandler::MeshRandomizer<MeshType>::AxisBoundaryCondition
     return m_mesh_line_node_boundary.nodeList();
   }
 
-  AxisBoundaryCondition(MeshLineNodeBoundary<const MeshType>&& mesh_line_node_boundary)
+  AxisBoundaryCondition(MeshLineNodeBoundary<MeshType>&& mesh_line_node_boundary)
     : m_mesh_line_node_boundary(mesh_line_node_boundary)
   {
     ;
@@ -338,7 +338,7 @@ class MeshRandomizerHandler::MeshRandomizer<MeshType>::SymmetryBoundaryCondition
   using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
-  const MeshFlatNodeBoundary<const MeshType> m_mesh_flat_node_boundary;
+  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
 
  public:
   const Rd&
@@ -353,7 +353,7 @@ class MeshRandomizerHandler::MeshRandomizer<MeshType>::SymmetryBoundaryCondition
     return m_mesh_flat_node_boundary.nodeList();
   }
 
-  SymmetryBoundaryCondition(MeshFlatNodeBoundary<const MeshType>&& mesh_flat_node_boundary)
+  SymmetryBoundaryCondition(MeshFlatNodeBoundary<MeshType>&& mesh_flat_node_boundary)
     : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
   {
     ;
@@ -369,7 +369,7 @@ MeshRandomizerHandler::getRandomizedMesh(
 {
   return std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(mesh)>;
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         MeshRandomizer randomizer(*mesh, bc_descriptor_list);
         return randomizer.getRandomizedMesh();
@@ -388,7 +388,7 @@ MeshRandomizerHandler::getRandomizedMesh(
 {
   return std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(mesh)>;
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         MeshRandomizer randomizer(*mesh, bc_descriptor_list);
         return randomizer.getRandomizedMesh(function_symbol_id);
diff --git a/src/mesh/MeshRelaxer.cpp b/src/mesh/MeshRelaxer.cpp
index 1192c2aed..7dcc86bc0 100644
--- a/src/mesh/MeshRelaxer.cpp
+++ b/src/mesh/MeshRelaxer.cpp
@@ -38,8 +38,8 @@ MeshRelaxer::relax(const std::shared_ptr<const MeshVariant>& p_source_mesh,
 {
   return std::visit(
     [&](auto&& source_mesh, auto&& destination_mesh) -> std::shared_ptr<const MeshVariant> {
-      using SourceMeshType      = typename std::decay_t<decltype(source_mesh)>::element_type;
-      using DestinationMeshType = typename std::decay_t<decltype(destination_mesh)>::element_type;
+      using SourceMeshType      = mesh_type_t<decltype(source_mesh)>;
+      using DestinationMeshType = mesh_type_t<decltype(destination_mesh)>;
       if constexpr (std::is_same_v<SourceMeshType, DestinationMeshType>) {
         return this->_relax(*source_mesh, *destination_mesh, theta);
       } else {
diff --git a/src/mesh/MeshSmoother.cpp b/src/mesh/MeshSmoother.cpp
index 22c8b41b6..eafc9b94b 100644
--- a/src/mesh/MeshSmoother.cpp
+++ b/src/mesh/MeshSmoother.cpp
@@ -328,7 +328,7 @@ class MeshSmootherHandler::MeshSmoother<MeshType>::AxisBoundaryCondition
   using Rd = TinyVector<Dimension, double>;
 
  private:
-  const MeshLineNodeBoundary<const MeshType> m_mesh_line_node_boundary;
+  const MeshLineNodeBoundary<MeshType> m_mesh_line_node_boundary;
 
  public:
   const Rd&
@@ -343,7 +343,7 @@ class MeshSmootherHandler::MeshSmoother<MeshType>::AxisBoundaryCondition
     return m_mesh_line_node_boundary.nodeList();
   }
 
-  AxisBoundaryCondition(MeshLineNodeBoundary<const MeshType>&& mesh_line_node_boundary)
+  AxisBoundaryCondition(MeshLineNodeBoundary<MeshType>&& mesh_line_node_boundary)
     : m_mesh_line_node_boundary(mesh_line_node_boundary)
   {
     ;
@@ -385,7 +385,7 @@ class MeshSmootherHandler::MeshSmoother<MeshType>::SymmetryBoundaryCondition
   using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
-  const MeshFlatNodeBoundary<const MeshType> m_mesh_flat_node_boundary;
+  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
 
  public:
   const Rd&
@@ -400,7 +400,7 @@ class MeshSmootherHandler::MeshSmoother<MeshType>::SymmetryBoundaryCondition
     return m_mesh_flat_node_boundary.nodeList();
   }
 
-  SymmetryBoundaryCondition(MeshFlatNodeBoundary<const MeshType>&& mesh_flat_node_boundary)
+  SymmetryBoundaryCondition(MeshFlatNodeBoundary<MeshType>&& mesh_flat_node_boundary)
     : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
   {
     ;
diff --git a/src/mesh/MeshTraits.hpp b/src/mesh/MeshTraits.hpp
index 976954ed1..f55b53934 100644
--- a/src/mesh/MeshTraits.hpp
+++ b/src/mesh/MeshTraits.hpp
@@ -1,7 +1,9 @@
 #ifndef MESH_TRAITS_HPP
 #define MESH_TRAITS_HPP
 
-#include <cstddef>
+#include <utils/PugsTraits.hpp>
+
+#include <memory>
 #include <type_traits>
 
 template <size_t Dimension>
@@ -22,6 +24,25 @@ constexpr bool is_mesh_v<T&> = is_mesh_v<std::remove_cvref_t<T>>;
 template <typename T>
 concept MeshConcept = is_mesh_v<T>;
 
+// Get mesh type (from std::shared_ptr for instance)
+
+template <typename T>
+struct mesh_type_t_from
+{
+  static_assert(is_mesh_v<T>, "Requires mesh type");
+  using type = std::remove_const_t<T>;
+};
+
+template <MeshConcept T, template <typename T1> typename C>
+struct mesh_type_t_from<C<T>>
+{
+  static_assert(is_shared_ptr_v<std::decay_t<C<T>>>, "Requires std::shared_ptr container");
+  using type = std::remove_const_t<T>;
+};
+
+template <typename T>
+using mesh_type_t = typename mesh_type_t_from<std::decay_t<T>>::type;
+
 // Check mesh kind
 template <MeshConcept MeshType>
 constexpr inline bool is_polygonal_mesh_v = false;
diff --git a/src/mesh/MeshTransformer.cpp b/src/mesh/MeshTransformer.cpp
index 74911b312..d44d3251c 100644
--- a/src/mesh/MeshTransformer.cpp
+++ b/src/mesh/MeshTransformer.cpp
@@ -2,6 +2,7 @@
 
 #include <mesh/Connectivity.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
 
 #include <language/utils/EvaluateAtPoints.hpp>
@@ -30,7 +31,7 @@ MeshTransformer::transform(const FunctionSymbolId& function_id, std::shared_ptr<
 {
   return std::visit(
     [&](auto&& mesh) {
-      using MeshType             = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType             = mesh_type_t<decltype(mesh)>;
       constexpr size_t Dimension = MeshType::Dimension;
       using TransformT           = TinyVector<Dimension>(TinyVector<Dimension>);
 
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index f846dd2eb..5415d59f8 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -5,6 +5,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
 #include <utils/Filesystem.hpp>
 #include <utils/Messenger.hpp>
@@ -320,7 +321,7 @@ GnuplotWriter::_writeAtTime(const MeshVariant& mesh_v,
 
   std::visit(
     [&](auto&& p_mesh) {
-      using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
       if constexpr ((MeshType::Dimension == 1) or (MeshType::Dimension == 2)) {
         this->_write(*p_mesh, output_named_item_data_set, time);
       } else {
@@ -337,7 +338,7 @@ GnuplotWriter::_writeMesh(const MeshVariant& mesh_v) const
 
   std::visit(
     [&](auto&& p_mesh) {
-      using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
       if constexpr ((MeshType::Dimension == 1) or (MeshType::Dimension == 2)) {
         this->_write(*p_mesh, output_named_item_data_set, {});
       } else {
@@ -355,7 +356,7 @@ GnuplotWriter::_write(const MeshVariant& mesh_v,
 
   std::visit(
     [&](auto&& p_mesh) {
-      using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
       if constexpr ((MeshType::Dimension == 1) or (MeshType::Dimension == 2)) {
         this->_write(*p_mesh, output_named_item_data_set, {});
       } else {
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index 745903222..a31d5432f 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -5,6 +5,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
 #include <utils/Filesystem.hpp>
 #include <utils/Messenger.hpp>
@@ -367,7 +368,7 @@ GnuplotWriter1D::_writeAtTime(const MeshVariant& mesh_v,
 
   std::visit(
     [&](auto&& p_mesh) {
-      using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
       if constexpr (MeshType::Dimension == 1) {
         this->_write(*p_mesh, output_named_item_data_set, time);
       } else if constexpr (MeshType::Dimension == 2) {
@@ -391,7 +392,7 @@ GnuplotWriter1D::_write(const MeshVariant& mesh_v,
 
   std::visit(
     [&](auto&& p_mesh) {
-      using MeshType = typename std::decay_t<decltype(p_mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
       if constexpr (MeshType::Dimension == 1) {
         this->_write(*p_mesh, output_named_item_data_set, {});
       } else if constexpr (MeshType::Dimension == 2) {
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 131854ec8..9f00bf12b 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -722,7 +722,7 @@ class AcousticSolverHandler::AcousticSolver<MeshType>::PressureBoundaryCondition
 };
 
 template <>
-class AcousticSolverHandler::AcousticSolver<const Mesh<1>>::PressureBoundaryCondition
+class AcousticSolverHandler::AcousticSolver<Mesh<1>>::PressureBoundaryCondition
 {
  private:
   const MeshNodeBoundary m_mesh_node_boundary;
@@ -871,7 +871,7 @@ AcousticSolverHandler::AcousticSolverHandler(const std::shared_ptr<const MeshVar
 
   std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(mesh)>;
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         m_acoustic_solver = std::make_unique<AcousticSolver<MeshType>>();
       } else {
diff --git a/src/scheme/DiscreteFunctionIntegrator.cpp b/src/scheme/DiscreteFunctionIntegrator.cpp
index 8b2913fe2..92dd526d2 100644
--- a/src/scheme/DiscreteFunctionIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionIntegrator.cpp
@@ -174,10 +174,11 @@ DiscreteFunctionIntegrator::integrate() const
 {
   return std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
-
+      using MeshType = mesh_type_t<decltype(mesh)>;
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         return this->_integrate<MeshType>();
+      } else {
+        throw NormalError("unexpected mesh type");
       }
     },
     m_mesh->variant());
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index 3894eaed8..6c5228a91 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -192,10 +192,12 @@ DiscreteFunctionInterpoler::interpolate() const
 {
   return std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(mesh)>;
 
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         return this->_interpolate<MeshType>();
+      } else {
+        throw NormalError("unexpected mesh type");
       }
     },
     m_mesh->variant());
diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.cpp b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
index 25b460a2e..74bb43d66 100644
--- a/src/scheme/DiscreteFunctionVectorIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
@@ -130,10 +130,12 @@ DiscreteFunctionVectorIntegrator::integrate() const
 
   return std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(mesh)>;
 
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         return this->_integrate<MeshType>();
+      } else {
+        throw NormalError("unexpected mesh type");
       }
     },
     m_mesh->variant());
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.cpp b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
index b30fe107b..547709c45 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
@@ -151,10 +151,12 @@ DiscreteFunctionVectorInterpoler::interpolate() const
 
   return std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(mesh)>;
 
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         return this->_interpolate<MeshType>();
+      } else {
+        throw NormalError("unexpected mesh type");
       }
     },
     m_mesh->variant());
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index e7ec1994b..e1bf11a8b 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -15,6 +15,7 @@
 #include <mesh/MeshFaceBoundary.hpp>
 #include <mesh/MeshFlatFaceBoundary.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/SubItemValuePerItem.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
@@ -240,7 +241,7 @@ FluxingAdvectionSolver<MeshType>::_computeDonorCells(FaceValue<const double> alg
 
 template <>
 void
-FluxingAdvectionSolver<const Mesh<1>>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
+FluxingAdvectionSolver<Mesh<1>>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
 {
   m_donnor_cell = [&] {
     const auto face_to_cell_matrix = m_new_mesh->connectivity().faceToCellMatrix();
@@ -274,7 +275,7 @@ FluxingAdvectionSolver<const Mesh<1>>::_computeDonorCells(FaceValue<const double
 
 template <>
 FaceValue<double>
-FluxingAdvectionSolver<const Mesh<1>>::_computeAlgebraicFluxingVolume()
+FluxingAdvectionSolver<Mesh<1>>::_computeAlgebraicFluxingVolume()
 {
   Array<double> fluxing_volumes{m_new_mesh->numberOfNodes()};
   NodeValue<double> nodal_fluxing_volume(m_new_mesh->connectivity(), fluxing_volumes);
@@ -293,7 +294,7 @@ FluxingAdvectionSolver<const Mesh<1>>::_computeAlgebraicFluxingVolume()
 
 template <>
 FaceValue<double>
-FluxingAdvectionSolver<const Mesh<2>>::_computeAlgebraicFluxingVolume()
+FluxingAdvectionSolver<Mesh<2>>::_computeAlgebraicFluxingVolume()
 {
   const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix();
   FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity());
@@ -320,7 +321,7 @@ FluxingAdvectionSolver<const Mesh<2>>::_computeAlgebraicFluxingVolume()
 
 template <>
 FaceValue<double>
-FluxingAdvectionSolver<const Mesh<3>>::_computeAlgebraicFluxingVolume()
+FluxingAdvectionSolver<Mesh<3>>::_computeAlgebraicFluxingVolume()
 {
   // due to the z component of the jacobian determinant, degree 3
   // polynomials must be exactly integrated
@@ -869,8 +870,8 @@ advectByFluxing(const std::shared_ptr<const MeshVariant> new_mesh_v,
 
   return std::visit(
     [&](auto&& old_mesh, auto&& new_mesh) -> std::vector<std::shared_ptr<const DiscreteFunctionVariant>> {
-      using OldMeshType = typename std::decay_t<decltype(old_mesh)>::element_type;
-      using NewMeshType = typename std::decay_t<decltype(new_mesh)>::element_type;
+      using OldMeshType = mesh_type_t<decltype(old_mesh)>;
+      using NewMeshType = mesh_type_t<decltype(new_mesh)>;
 
       if constexpr (std::is_same_v<OldMeshType, NewMeshType>) {
         return FluxingAdvectionSolver<OldMeshType>{old_mesh, new_mesh}.remap(remapped_variables_with_bc);
diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index c0d011fed..4fc391885 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -793,7 +793,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::PressureBoundaryC
 };
 
 template <>
-class HyperelasticSolverHandler::HyperelasticSolver<const Mesh<1>>::PressureBoundaryCondition
+class HyperelasticSolverHandler::HyperelasticSolver<Mesh<1>>::PressureBoundaryCondition
 {
  private:
   const MeshNodeBoundary m_mesh_node_boundary;
@@ -847,7 +847,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::NormalStressBound
 };
 
 template <>
-class HyperelasticSolverHandler::HyperelasticSolver<const Mesh<1>>::NormalStressBoundaryCondition
+class HyperelasticSolverHandler::HyperelasticSolver<Mesh<1>>::NormalStressBoundaryCondition
 {
  private:
   const MeshNodeBoundary m_mesh_node_boundary;
@@ -947,7 +947,7 @@ HyperelasticSolverHandler::HyperelasticSolverHandler(const std::shared_ptr<const
 
   std::visit(
     [&](auto&& mesh) {
-      using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
+      using MeshType = mesh_type_t<decltype(mesh)>;
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         m_hyperelastic_solver = std::make_unique<HyperelasticSolver<MeshType>>();
       } else {
-- 
GitLab