diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 1eec3fb02bb419f5ef23f733afe35e0697787982..d4d26979451127cacf44b0719b62588d7ee108e7 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -246,7 +246,7 @@ SchemeModule::SchemeModule()
                                                const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                                  bc_descriptor_list) -> std::shared_ptr<const IMesh> {
                                               MeshSmootherHandler handler;
-                                              return handler.getSmoothedMesh(*p_mesh, bc_descriptor_list);
+                                              return handler.getSmoothedMesh(p_mesh, bc_descriptor_list);
                                             }
 
                                             ));
@@ -259,7 +259,7 @@ SchemeModule::SchemeModule()
                                    bc_descriptor_list,
                                  const FunctionSymbolId& function_symbol_id) -> std::shared_ptr<const IMesh> {
                                 MeshSmootherHandler handler;
-                                return handler.getSmoothedMesh(*p_mesh, bc_descriptor_list, function_symbol_id);
+                                return handler.getSmoothedMesh(p_mesh, bc_descriptor_list, function_symbol_id);
                               }
 
                               ));
@@ -273,7 +273,7 @@ SchemeModule::SchemeModule()
                                  const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list)
                                 -> std::shared_ptr<const IMesh> {
                                 MeshSmootherHandler handler;
-                                return handler.getSmoothedMesh(*p_mesh, bc_descriptor_list, smoothing_zone_list);
+                                return handler.getSmoothedMesh(p_mesh, bc_descriptor_list, smoothing_zone_list);
                               }
 
                               ));
@@ -286,7 +286,7 @@ SchemeModule::SchemeModule()
                                                const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>&
                                                  discrete_function_variant_list) -> std::shared_ptr<const IMesh> {
                                               MeshSmootherHandler handler;
-                                              return handler.getSmoothedMesh(*p_mesh, bc_descriptor_list,
+                                              return handler.getSmoothedMesh(p_mesh, bc_descriptor_list,
                                                                              discrete_function_variant_list);
                                             }
 
diff --git a/src/mesh/MeshSmoother.cpp b/src/mesh/MeshSmoother.cpp
index befa5cdf109991c155cf4d13048e8e702629aa24..a9af80c0129fdb20267b3c62e4430535f00f154f 100644
--- a/src/mesh/MeshSmoother.cpp
+++ b/src/mesh/MeshSmoother.cpp
@@ -11,6 +11,7 @@
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 #include <scheme/AxisBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
@@ -408,26 +409,26 @@ class MeshSmootherHandler::MeshSmoother<Dimension>::SymmetryBoundaryCondition
 
 std::shared_ptr<const IMesh>
 MeshSmootherHandler::getSmoothedMesh(
-  const IMesh& mesh,
+  const std::shared_ptr<const IMesh>& mesh,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
 {
-  switch (mesh.dimension()) {
+  switch (mesh->dimension()) {
   case 1: {
     constexpr size_t Dimension = 1;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh();
   }
   case 2: {
     constexpr size_t Dimension = 2;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh();
   }
   case 3: {
     constexpr size_t Dimension = 3;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh();
   }
   default: {
@@ -438,27 +439,27 @@ MeshSmootherHandler::getSmoothedMesh(
 
 std::shared_ptr<const IMesh>
 MeshSmootherHandler::getSmoothedMesh(
-  const IMesh& mesh,
+  const std::shared_ptr<const IMesh>& mesh,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
   const FunctionSymbolId& function_symbol_id) const
 {
-  switch (mesh.dimension()) {
+  switch (mesh->dimension()) {
   case 1: {
     constexpr size_t Dimension = 1;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh(function_symbol_id);
   }
   case 2: {
     constexpr size_t Dimension = 2;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh(function_symbol_id);
   }
   case 3: {
     constexpr size_t Dimension = 3;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh(function_symbol_id);
   }
   default: {
@@ -469,27 +470,27 @@ MeshSmootherHandler::getSmoothedMesh(
 
 std::shared_ptr<const IMesh>
 MeshSmootherHandler::getSmoothedMesh(
-  const IMesh& mesh,
+  const std::shared_ptr<const IMesh>& mesh,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
   const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const
 {
-  switch (mesh.dimension()) {
+  switch (mesh->dimension()) {
   case 1: {
     constexpr size_t Dimension = 1;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh(smoothing_zone_list);
   }
   case 2: {
     constexpr size_t Dimension = 2;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh(smoothing_zone_list);
   }
   case 3: {
     constexpr size_t Dimension = 3;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh(smoothing_zone_list);
   }
   default: {
@@ -500,27 +501,37 @@ MeshSmootherHandler::getSmoothedMesh(
 
 std::shared_ptr<const IMesh>
 MeshSmootherHandler::getSmoothedMesh(
-  const IMesh& mesh,
+  const std::shared_ptr<const IMesh>& mesh,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
   const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
 {
-  switch (mesh.dimension()) {
+  if (not hasSameMesh(discrete_function_variant_list)) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  std::shared_ptr<const IMesh> common_mesh = getCommonMesh(discrete_function_variant_list);
+
+  if (common_mesh != mesh) {
+    throw NormalError("discrete functions are not defined on the smoothed mesh");
+  }
+
+  switch (mesh->dimension()) {
   case 1: {
     constexpr size_t Dimension = 1;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh(discrete_function_variant_list);
   }
   case 2: {
     constexpr size_t Dimension = 2;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh(discrete_function_variant_list);
   }
   case 3: {
     constexpr size_t Dimension = 3;
     using MeshType             = Mesh<Connectivity<Dimension>>;
-    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
     return smoother.getSmoothedMesh(discrete_function_variant_list);
   }
   default: {
diff --git a/src/mesh/MeshSmoother.hpp b/src/mesh/MeshSmoother.hpp
index 74fde4855809f2c3610e1632f47bfb463a3f5b43..99db9bf76bf9ea771b1ebbc41c1b948c6d451632 100644
--- a/src/mesh/MeshSmoother.hpp
+++ b/src/mesh/MeshSmoother.hpp
@@ -19,21 +19,21 @@ class MeshSmootherHandler
 
  public:
   std::shared_ptr<const IMesh> getSmoothedMesh(
-    const IMesh& mesh,
+    const std::shared_ptr<const IMesh>& mesh,
     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const;
 
   std::shared_ptr<const IMesh> getSmoothedMesh(
-    const IMesh& mesh,
+    const std::shared_ptr<const IMesh>& mesh,
     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
     const FunctionSymbolId& function_symbol_id) const;
 
   std::shared_ptr<const IMesh> getSmoothedMesh(
-    const IMesh& mesh,
+    const std::shared_ptr<const IMesh>& mesh,
     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
     const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const;
 
   std::shared_ptr<const IMesh> getSmoothedMesh(
-    const IMesh& mesh,
+    const std::shared_ptr<const IMesh>& mesh,
     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
     const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& smoothing_zone_list) const;