From bc654ee555d1411f53c4326d708608207b20fb1c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 20 Sep 2021 18:16:21 +0200
Subject: [PATCH] Add a mesh interpolation functionality

---
 src/language/modules/MeshModule.cpp | 14 +++++++
 src/mesh/CMakeLists.txt             |  1 +
 src/mesh/MeshInterpoler.cpp         | 59 +++++++++++++++++++++++++++++
 src/mesh/MeshInterpoler.hpp         | 27 +++++++++++++
 4 files changed, 101 insertions(+)
 create mode 100644 src/mesh/MeshInterpoler.cpp
 create mode 100644 src/mesh/MeshInterpoler.hpp

diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 86c9a0c10..5d2c577cb 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 c315ced7e..a6ededf56 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 000000000..26cba4bff
--- /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 000000000..f8a79eea4
--- /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
-- 
GitLab