Merge branch 'feature/mesh-interpoler' into 'develop'

Feature/mesh interpoler

See merge request !119

src/language/modules/MeshModule.cpp

@@ -12,52 +12,12 @@
 #include <mesh/DiamondDualMeshManager.hpp>
 #include <mesh/GmshReader.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshInterpoler.hpp>
+#include <mesh/MeshTransformer.hpp>
 #include <utils/Exceptions.hpp>
 
 #include <Kokkos_Core.hpp>
 
-template <typename T>
-class MeshTransformation;
-template <typename OutputType, typename InputType>
-class MeshTransformation<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
-{
-  static constexpr size_t Dimension = OutputType::Dimension;
-  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
-
- public:
-  static inline std::shared_ptr<Mesh<Connectivity<Dimension>>>
-  transform(const FunctionSymbolId& function_symbol_id, std::shared_ptr<const IMesh> p_mesh)
-  {
-    using MeshType             = Mesh<Connectivity<Dimension>>;
-    const MeshType& given_mesh = dynamic_cast<const MeshType&>(*p_mesh);
-
-    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
-    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
-
-    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
-
-    NodeValue<const InputType> given_xr = given_mesh.xr();
-    NodeValue<OutputType> xr(given_mesh.connectivity());
-
-    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
-    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
-
-    parallel_for(given_mesh.numberOfNodes(), [=, &expression, &tokens](NodeId r) {
-      const int32_t t = tokens.acquire();
-
-      auto& execution_policy = context_list[t];
-
-      Adapter::convertArgs(execution_policy.currentContext(), given_xr[r]);
-      auto result = expression.execute(execution_policy);
-      xr[r]       = convert_result(std::move(result));
-
-      tokens.release(t);
-    });
-
-    return std::make_shared<MeshType>(given_mesh.shared_connectivity(), xr);
-  }
-};
-
 MeshModule::MeshModule()
 {
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IMesh>>);
@@ -78,23 +38,20 @@ MeshModule::MeshModule()
 
                               [](std::shared_ptr<const IMesh> p_mesh,
                                  const FunctionSymbolId& function_id) -> std::shared_ptr<const IMesh> {
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  using TransformT = TinyVector<1>(TinyVector<1>);
-                                  return MeshTransformation<TransformT>::transform(function_id, p_mesh);
-                                }
-                                case 2: {
-                                  using TransformT = TinyVector<2>(TinyVector<2>);
-                                  return MeshTransformation<TransformT>::transform(function_id, p_mesh);
-                                }
-                                case 3: {
-                                  using TransformT = TinyVector<3>(TinyVector<3>);
-                                  return MeshTransformation<TransformT>::transform(function_id, p_mesh);
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
+                                return MeshTransformer{}.transform(function_id, p_mesh);
                               }
+
+                              ));
+
+  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);
                               }
 
                               ));

src/mesh/CMakeLists.txt

@@ -20,9 +20,11 @@ add_library(
   MeshFaceBoundary.cpp
   MeshFlatFaceBoundary.cpp
   MeshFlatNodeBoundary.cpp
+  MeshInterpoler.cpp
   MeshLineNodeBoundary.cpp
   MeshNodeBoundary.cpp
   MeshRandomizer.cpp
+  MeshTransformer.cpp
   SynchronizerManager.cpp)
 
 # Additional dependencies

src/mesh/MeshInterpoler.cpp

+#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");
+    }
+    }
+  }
+}

src/mesh/MeshInterpoler.hpp

+#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

src/mesh/MeshTransformer.cpp

+#include <mesh/MeshTransformer.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+
+#include <language/utils/FunctionTable.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+template <typename OutputType, typename InputType>
+class MeshTransformer::MeshTransformation<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+{
+  static constexpr size_t Dimension = OutputType::Dimension;
+  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
+
+ public:
+  static inline std::shared_ptr<Mesh<Connectivity<Dimension>>>
+  transform(const FunctionSymbolId& function_symbol_id, std::shared_ptr<const IMesh> p_mesh)
+  {
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    const MeshType& given_mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
+    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+
+    NodeValue<const InputType> given_xr = given_mesh.xr();
+    NodeValue<OutputType> xr(given_mesh.connectivity());
+
+    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+
+    parallel_for(given_mesh.numberOfNodes(), [=, &expression, &tokens](NodeId r) {
+      const int32_t t = tokens.acquire();
+
+      auto& execution_policy = context_list[t];
+
+      Adapter::convertArgs(execution_policy.currentContext(), given_xr[r]);
+      auto result = expression.execute(execution_policy);
+      xr[r]       = convert_result(std::move(result));
+
+      tokens.release(t);
+    });
+
+    return std::make_shared<MeshType>(given_mesh.shared_connectivity(), xr);
+  }
+};
+
+std::shared_ptr<const IMesh>
+MeshTransformer::transform(const FunctionSymbolId& function_id, std::shared_ptr<const IMesh> p_mesh)
+
+{
+  switch (p_mesh->dimension()) {
+  case 1: {
+    using TransformT = TinyVector<1>(TinyVector<1>);
+    return MeshTransformation<TransformT>::transform(function_id, p_mesh);
+  }
+  case 2: {
+    using TransformT = TinyVector<2>(TinyVector<2>);
+    return MeshTransformation<TransformT>::transform(function_id, p_mesh);
+  }
+  case 3: {
+    using TransformT = TinyVector<3>(TinyVector<3>);
+    return MeshTransformation<TransformT>::transform(function_id, p_mesh);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}

src/mesh/MeshTransformer.hpp

+#ifndef MESH_TRANSFORMER_HPP
+#define MESH_TRANSFORMER_HPP
+
+class IMesh;
+
+template <typename ConnectivityType>
+class Mesh;
+
+class FunctionSymbolId;
+
+#include <memory>
+
+class MeshTransformer
+{
+  template <typename T>
+  class MeshTransformation;
+
+ public:
+  std::shared_ptr<const IMesh> transform(const FunctionSymbolId& function_symbol_id,
+                                         std::shared_ptr<const IMesh> p_mesh);
+
+  MeshTransformer()  = default;
+  ~MeshTransformer() = default;
+};
+
+#endif   // MESH_TRANSFORMER_HPP