diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 86c9a0c106e9fe6eece190c93365fc42a47bf864..5d2c577cbff51ec66e5c92c41ec1e757b94a0313 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -12,6 +12,7 @@
 #include <mesh/DiamondDualMeshManager.hpp>
 #include <mesh/GmshReader.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshInterpoler.hpp>
 #include <utils/Exceptions.hpp>
 
 #include <Kokkos_Core.hpp>
@@ -99,6 +100,19 @@ MeshModule::MeshModule()
 
                               ));
 
+  this->_addBuiltinFunction("interpolate",
+                            std::make_shared<BuiltinFunctionEmbedder<
+                              std::shared_ptr<const IMesh>(const std::shared_ptr<const IMesh>&,
+                                                           const std::shared_ptr<const IMesh>&, const double&)>>(
+
+                              [](const std::shared_ptr<const IMesh>& source_mesh,
+                                 const std::shared_ptr<const IMesh>& destination_mesh,
+                                 const double& theta) -> std::shared_ptr<const IMesh> {
+                                return MeshInterpoler{}.interpolate(source_mesh, destination_mesh, theta);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("cartesian1dMesh",
                             std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                               const IMesh>(const TinyVector<1>, const TinyVector<1>, const std::vector<uint64_t>&)>>(
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index c315ced7e506310c08c414d78b4df77aab646d4b..a6ededf56bbb4b0ad28fed7cf0d09eab97490350 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -20,6 +20,7 @@ add_library(
   MeshFaceBoundary.cpp
   MeshFlatFaceBoundary.cpp
   MeshFlatNodeBoundary.cpp
+  MeshInterpoler.cpp
   MeshLineNodeBoundary.cpp
   MeshNodeBoundary.cpp
   MeshRandomizer.cpp
diff --git a/src/mesh/MeshInterpoler.cpp b/src/mesh/MeshInterpoler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..26cba4bff5eeca4b75690b47e30126bac070646b
--- /dev/null
+++ b/src/mesh/MeshInterpoler.cpp
@@ -0,0 +1,59 @@
+#include <mesh/MeshInterpoler.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+
+template <typename ConnectivityType>
+std::shared_ptr<const Mesh<ConnectivityType>>
+MeshInterpoler::_interpolate(const Mesh<ConnectivityType>& source_mesh,
+                             const Mesh<ConnectivityType>& destination_mesh,
+                             const double& theta) const
+{
+  if (source_mesh.shared_connectivity() == destination_mesh.shared_connectivity()) {
+    const ConnectivityType& connectivity = source_mesh.connectivity();
+    NodeValue<TinyVector<ConnectivityType::Dimension>> theta_xr{connectivity};
+    const NodeValue<const TinyVector<ConnectivityType::Dimension>> source_xr      = source_mesh.xr();
+    const NodeValue<const TinyVector<ConnectivityType::Dimension>> destination_xr = destination_mesh.xr();
+
+    const double one_minus_theta = 1 - theta;
+    parallel_for(
+      connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        theta_xr[node_id] = one_minus_theta * source_xr[node_id] + theta * destination_xr[node_id];
+      });
+
+    return std::make_shared<Mesh<ConnectivityType>>(source_mesh.shared_connectivity(), theta_xr);
+  } else {
+    throw NormalError("interpolated meshes must share the same connectivity");
+  }
+}
+
+std::shared_ptr<const IMesh>
+MeshInterpoler::interpolate(const std::shared_ptr<const IMesh>& p_source_mesh,
+                            const std::shared_ptr<const IMesh>& p_destination_mesh,
+                            const double& theta) const
+{
+  if (p_source_mesh->dimension() != p_destination_mesh->dimension()) {
+    throw NormalError("incompatible mesh dimensions");
+  } else {
+    switch (p_source_mesh->dimension()) {
+    case 1: {
+      using MeshType = Mesh<Connectivity<1>>;
+      return this->_interpolate(dynamic_cast<const MeshType&>(*p_source_mesh),
+                                dynamic_cast<const MeshType&>(*p_destination_mesh), theta);
+    }
+    case 2: {
+      using MeshType = Mesh<Connectivity<2>>;
+      return this->_interpolate(dynamic_cast<const MeshType&>(*p_source_mesh),
+                                dynamic_cast<const MeshType&>(*p_destination_mesh), theta);
+    }
+    case 3: {
+      using MeshType = Mesh<Connectivity<3>>;
+      return this->_interpolate(dynamic_cast<const MeshType&>(*p_source_mesh),
+                                dynamic_cast<const MeshType&>(*p_destination_mesh), theta);
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  }
+}
diff --git a/src/mesh/MeshInterpoler.hpp b/src/mesh/MeshInterpoler.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f8a79eea4918bf71b6b7586683a14f2a67194df9
--- /dev/null
+++ b/src/mesh/MeshInterpoler.hpp
@@ -0,0 +1,27 @@
+#ifndef MESH_INTERPOLER_HPP
+#define MESH_INTERPOLER_HPP
+
+class IMesh;
+
+template <typename ConnectivityType>
+class Mesh;
+
+#include <memory>
+
+class MeshInterpoler
+{
+ private:
+  template <typename ConnectivityType>
+  std::shared_ptr<const Mesh<ConnectivityType>> _interpolate(const Mesh<ConnectivityType>& source_mesh,
+                                                             const Mesh<ConnectivityType>& destination_mesh,
+                                                             const double& theta) const;
+
+ public:
+  std::shared_ptr<const IMesh> interpolate(const std::shared_ptr<const IMesh>& p_source_mesh,
+                                           const std::shared_ptr<const IMesh>& p_destination_mesh,
+                                           const double& theta) const;
+  MeshInterpoler()  = default;
+  ~MeshInterpoler() = default;
+};
+
+#endif   // MESH_INTERPOLER_HPP