Add polynomial mesh builder from a polygonal mesh and a degree

Also add simple gnuplot output

src/language/modules/MeshModule.cpp

@@ -31,6 +31,7 @@
 #include <mesh/NumberedBoundaryDescriptor.hpp>
 #include <mesh/NumberedInterfaceDescriptor.hpp>
 #include <mesh/NumberedZoneDescriptor.hpp>
+#include <mesh/PolynomialMeshBuilder.hpp>
 #include <mesh/SubItemArrayPerItemVariant.hpp>
 #include <mesh/SubItemValuePerItemVariant.hpp>
 #include <utils/Exceptions.hpp>
@@ -287,6 +288,16 @@ MeshModule::MeshModule()
                                               return DualMeshManager::instance().getMedianDualMesh(mesh_v);
                                             }
 
+                                            ));
+
+  this->_addBuiltinFunction("polynomialMesh", std::function(
+
+                                                [](const std::shared_ptr<const MeshVariant>& mesh_v,
+                                                   const uint64_t degree) -> std::shared_ptr<const MeshVariant> {
+                                                  PolynomialMeshBuilder builder{mesh_v, degree};
+                                                  return builder.mesh();
+                                                }
+
                                                 ));
 }
 

src/language/utils/EvaluateAtPoints.hpp

@@ -22,7 +22,6 @@ class EvaluateAtPoints<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
                   "invalid input data type");
     static_assert(std::is_same_v<std::remove_const_t<typename OutputArrayT::data_type>, OutputType>,
                   "invalid output data type");
-    Assert(size(value) == size(position));
 
     auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
     auto convert_result = Adapter::getResultConverter(expression.m_data_type);
@@ -39,7 +38,12 @@ class EvaluateAtPoints<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     } else {
       static_assert(std::is_same_v<OutputType, double>, "unexpected output type");
     }
-
+    static_assert(is_item_value_v<InputArrayT> or is_array_v<InputArrayT> or is_item_array_v<InputArrayT> or
+                    is_table_v<InputArrayT>,
+                  "TEMPORARY TEST");
+    if constexpr (is_item_value_v<InputArrayT> or is_array_v<InputArrayT>) {
+      static_assert(is_item_value_v<OutputArrayT> or is_array_v<OutputArrayT>, "incompatible intput and output types");
+      Assert(size(value) == size(position));
       parallel_for(size(position), [=, &expression, &tokens](typename InputArrayT::index_type i) {
         const int32_t t = tokens.acquire();
 
@@ -51,6 +55,26 @@ class EvaluateAtPoints<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
 
         tokens.release(t);
       });
+    } else if constexpr (is_item_array_v<InputArrayT> or is_table_v<InputArrayT>) {
+      static_assert(is_item_array_v<OutputArrayT> or is_table_v<OutputArrayT>, "incompatible intput and output types");
+      Assert(numberOfRows(position) == numberOfRows(value));
+      Assert(numberOfColumns(position) == numberOfColumns(value));
+      parallel_for(numberOfRows(position), [=, &expression, &tokens](typename InputArrayT::index_type i) {
+        const int32_t t = tokens.acquire();
+
+        auto& execution_policy = context_list[t];
+
+        const auto positions = position[i];
+        auto values          = value[i];
+
+        for (size_t j = 0; j < positions.size(); ++j) {
+          Adapter::convertArgs(execution_policy.currentContext(), positions[j]);
+          auto result = expression.execute(execution_policy);
+          values[j]   = convert_result(std::move(result));
+        }
+        tokens.release(t);
+      });
+    }
   }
 
   template <class InputArrayT>

src/mesh/CMakeLists.txt

@@ -44,6 +44,7 @@ add_library(
   MeshUtils.cpp
   MeshVariant.cpp
   PolynomialMesh.cpp
+  PolynomialMeshBuilder.cpp
   SynchronizerManager.cpp
 )
 

src/mesh/MeshTransformer.cpp

@@ -25,6 +25,26 @@ class MeshTransformer::MeshTransformation<OutputType(InputType)>
 
     return std::make_shared<const Mesh<Dimension>>(mesh.shared_connectivity(), xr);
   }
+
+  template <size_t Dimension>
+  static std::shared_ptr<const PolynomialMesh<Dimension>>
+  transform(const FunctionSymbolId& function_symbol_id, const PolynomialMesh<Dimension>& mesh)
+  {
+    if constexpr (Dimension == 2) {
+      NodeValue<OutputType> xr(mesh.connectivity());
+      NodeValue<const InputType> given_xr = mesh.xr();
+      EvaluateAtPoints<OutputType(InputType)>::evaluateTo(function_symbol_id, given_xr, xr);
+
+      FaceArray<const InputType> given_xl = mesh.xl();
+      FaceArray<OutputType> xl(mesh.connectivity(), given_xl.sizeOfArrays());
+
+      EvaluateAtPoints<OutputType(InputType)>::evaluateTo(function_symbol_id, given_xl, xl);
+
+      return std::make_shared<const PolynomialMesh<Dimension>>(mesh.shared_connectivity(), xr, xl);
+    } else {
+      throw NotImplementedError("unexpected dimension for polynomial mesh");
+    }
+  }
 };
 
 std::shared_ptr<const MeshVariant>
@@ -36,11 +56,11 @@ MeshTransformer::transform(const FunctionSymbolId& function_id, std::shared_ptr<
       using MeshType             = mesh_type_t<decltype(mesh)>;
       constexpr size_t Dimension = MeshType::Dimension;
       using TransformT           = TinyVector<Dimension>(TinyVector<Dimension>);
-      if constexpr (is_polygonal_mesh_v<MeshType>) {
+      // if constexpr (is_polygonal_mesh_v<MeshType>) {
       return std::make_shared<MeshVariant>(MeshTransformation<TransformT>::transform(function_id, *mesh));
-      } else {
-        throw NotImplementedError("unexpected mesh type " + demangle<MeshType>());
-      }
+      // } else {
+      //   throw NotImplementedError("unexpected mesh type " + demangle<MeshType>());
+      // }
     },
     mesh_v->variant());
 }

src/mesh/PolynomialMesh.hpp

@@ -24,7 +24,7 @@ class PolynomialMesh : public std::enable_shared_from_this<PolynomialMesh<MeshDi
 
   const std::shared_ptr<const Connectivity> m_connectivity;
   const NodeValue<const Rd> m_xr;
-  const EdgeArray<const Rd> m_xl;
+  const FaceArray<const Rd> m_xl;
 
   std::ostream&
   _write(std::ostream& os) const
@@ -131,7 +131,7 @@ class PolynomialMesh : public std::enable_shared_from_this<PolynomialMesh<MeshDi
   }
 
   PUGS_INLINE
-  const EdgeArray<const Rd>&
+  const FaceArray<const Rd>&
   xl() const
   {
     return m_xl;
@@ -142,7 +142,7 @@ class PolynomialMesh : public std::enable_shared_from_this<PolynomialMesh<MeshDi
   PUGS_INLINE
   PolynomialMesh(const std::shared_ptr<const Connectivity>& connectivity,
                  const NodeValue<const Rd>& xr,
-                 const EdgeArray<const Rd>& xl)
+                 const FaceArray<const Rd>& xl)
     : m_id{GlobalVariableManager::instance().getAndIncrementMeshId()}, m_connectivity{connectivity}, m_xr{xr}, m_xl{xl}
   {
     ;

src/mesh/PolynomialMeshBuilder.cpp

+#include <mesh/PolynomialMeshBuilder.hpp>
+
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <mesh/PolynomialMesh.hpp>
+
+#include <string>
+
+template <MeshConcept MeshType>
+void
+PolynomialMeshBuilder::_build(const MeshType& mesh, const uint64_t& degree)
+{
+  constexpr size_t Dimension = MeshType::Dimension;
+  static_assert(Dimension == 2);
+
+  using Rd = TinyVector<Dimension>;
+
+  if (degree < 2) {
+    throw NormalError("provided degree must be at least 2");
+  }
+
+  const Connectivity<Dimension>& connectivity = mesh.connectivity();
+
+  const NodeValue<const Rd> xr = mesh.xr();
+  FaceArray<Rd> xl(connectivity, degree - 1);
+
+  auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+
+  const size_t number_of_face_nodes = degree - 1;
+  const double ratio                = 1. / degree;
+
+  parallel_for(
+    connectivity.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+      const Rd x0 = xr[face_to_node_matrix[face_id][0]];
+      const Rd x1 = xr[face_to_node_matrix[face_id][1]];
+
+      for (size_t i_face_node = 0; i_face_node < number_of_face_nodes; ++i_face_node) {
+        xl[face_id][i_face_node] = (i_face_node + 1) * ratio * x1 + (number_of_face_nodes - i_face_node) * ratio * x0;
+      }
+    });
+
+  std::shared_ptr<const PolynomialMesh<Dimension>> polynomial_mesh =
+    std::make_shared<PolynomialMesh<Dimension>>(mesh.shared_connectivity(), xr, xl);
+
+  m_mesh = std::make_shared<MeshVariant>(polynomial_mesh);
+}
+
+PolynomialMeshBuilder::PolynomialMeshBuilder(const std::shared_ptr<const MeshVariant>& mesh_v, const uint64_t& degree)
+{
+  std::visit(
+    [&](auto&& mesh) {
+      using MeshType = mesh_type_t<decltype(mesh)>;
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+        if constexpr (MeshType::Dimension == 2) {
+          this->_build(*mesh, degree);
+        } else {
+          throw NotImplementedError("cannot build polynomial mesh in dimension " + std::to_string(MeshType::Dimension));
+        }
+      } else {
+        throw NormalError("input mesh must be a polygonal mesh. Provided: " + demangle<MeshType>());
+      }
+    },
+    mesh_v->variant());
+}

src/mesh/PolynomialMeshBuilder.hpp

+#ifndef POLYNOMIAL_MESH_BUILDER_HPP
+#define POLYNOMIAL_MESH_BUILDER_HPP
+
+#include <mesh/MeshTraits.hpp>
+
+class MeshVariant;
+
+#include <memory>
+
+class PolynomialMeshBuilder
+{
+ protected:
+  std::shared_ptr<const MeshVariant> m_mesh;
+
+  template <MeshConcept MeshType>
+  void _build(const MeshType& mesh, const uint64_t& degree);
+
+ public:
+  std::shared_ptr<const MeshVariant>
+  mesh() const
+  {
+    return m_mesh;
+  }
+
+  PolynomialMeshBuilder(const std::shared_ptr<const MeshVariant>& mesh_v, const uint64_t& degree);
+  ~PolynomialMeshBuilder() = default;
+};
+
+#endif   // POLYNOMIAL_MESH_BUILDER_HPP

src/output/GnuplotWriter.cpp

@@ -176,9 +176,10 @@ GnuplotWriter::_writeDataAtNodes(const MeshType& mesh,
       }
     }
   } else if constexpr (MeshType::Dimension == 2) {
-    auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
     auto cell_is_owned = mesh.connectivity().cellIsOwned();
 
+    if constexpr (is_polygonal_mesh_v<MeshType>) {
+      auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
       for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
         if (cell_is_owned[cell_id]) {
           const auto& cell_nodes = cell_to_node_matrix[cell_id];
@@ -200,6 +201,43 @@ GnuplotWriter::_writeDataAtNodes(const MeshType& mesh,
           fout << "\n\n\n";
         }
       }
+    } else {
+      static_assert(is_polynomial_mesh_v<MeshType>);
+      auto cell_to_face_matrix   = mesh.connectivity().cellToFaceMatrix();
+      auto face_to_node_matrix   = mesh.connectivity().faceToNodeMatrix();
+      auto cell_face_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+      if (output_named_item_data_set.size() > 0) {
+        throw NotImplementedError("cannot store data on polynomial mesh");
+      }
+
+      for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+        if (cell_is_owned[cell_id]) {
+          const auto& cell_faces = cell_to_face_matrix[cell_id];
+          for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
+            const FaceId& face_id              = cell_faces[i_face];
+            TinyVector<MeshType::Dimension> x0 = mesh.xr()[face_to_node_matrix[face_id][0]];
+            TinyVector<MeshType::Dimension> x1 = mesh.xr()[face_to_node_matrix[face_id][1]];
+            if (not cell_face_is_reversed(cell_id, i_face)) {
+              fout << x0[0] << ' ' << x0[1] << '\n';
+              for (size_t i_face_node = 0; i_face_node < mesh.xl().sizeOfArrays(); ++i_face_node) {
+                fout << mesh.xl()[face_id][i_face_node][0] << ' ' << mesh.xl()[face_id][i_face_node][1] << '\n';
+              }
+              fout << x1[0] << ' ' << x1[1] << '\n';
+            } else {
+              fout << x1[0] << ' ' << x1[1] << '\n';
+              for (size_t i_face_node = 0; i_face_node < mesh.xl().sizeOfArrays(); ++i_face_node) {
+                const size_t i_reversed_face_node = mesh.xl().sizeOfArrays() - 1 - i_face_node;
+                fout << mesh.xl()[face_id][i_reversed_face_node][0] << ' '
+                     << mesh.xl()[face_id][i_reversed_face_node][1] << '\n';
+              }
+              fout << x0[0] << ' ' << x0[1] << '\n';
+            }
+          }
+          fout << "\n\n\n";
+        }
+      }
+    }
   } else {
     throw UnexpectedError("invalid mesh dimension");
   }
@@ -344,15 +382,11 @@ GnuplotWriter::_writeMesh(const MeshVariant& mesh_v) const
   std::visit(
     [&](auto&& p_mesh) {
       using MeshType = mesh_type_t<decltype(p_mesh)>;
-      if constexpr (is_polygonal_mesh_v<MeshType>) {
       if constexpr ((MeshType::Dimension == 1) or (MeshType::Dimension == 2)) {
         this->_write(*p_mesh, output_named_item_data_set, {});
       } else {
         throw NormalError("gnuplot format is not available in dimension " + stringify(MeshType::Dimension));
       }
-      } else {
-        throw NormalError("gnuplot format is not available for this mesh type " + demangle<MeshType>());
-      }
     },
     mesh_v.variant());
 }

src/output/OutputNamedItemValueSet.hpp

@@ -138,6 +138,12 @@ class OutputNamedItemDataSet
   std::vector<std::pair<std::string, ItemDataVariant>> m_name_itemvariant_list;
 
  public:
+  size_t
+  size() const
+  {
+    return m_name_itemvariant_list.size();
+  }
+
   auto
   begin() const
   {