diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 86c9a0c106e9fe6eece190c93365fc42a47bf864..61767078aeaa0f8362393196a01ce9cde167a064 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/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);
}
));
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index c315ced7e506310c08c414d78b4df77aab646d4b..f536bae7afe8b6d1bc979e4abd9feccaf7c86a7d 100644
--- a/src/mesh/CMakeLists.txt
+++ b/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
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
diff --git a/src/mesh/MeshTransformer.cpp b/src/mesh/MeshTransformer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..67b3688926dcfb765c44799cff4674955515eb8b
--- /dev/null
+++ b/src/mesh/MeshTransformer.cpp
@@ -0,0 +1,71 @@
+#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");
+ }
+ }
+}
diff --git a/src/mesh/MeshTransformer.hpp b/src/mesh/MeshTransformer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ce0a4d47eab0557004234d22ea5331a0d2b0f793
--- /dev/null
+++ b/src/mesh/MeshTransformer.hpp
@@ -0,0 +1,26 @@
+#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