Add zone or indicator function for smoothed zone

src/language/modules/SchemeModule.cpp

@@ -264,6 +264,34 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("smoothMesh",
+                            std::function(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list)
+                                -> std::shared_ptr<const IMesh> {
+                                MeshSmootherHandler handler;
+                                return handler.getSmoothedMesh(*p_mesh, bc_descriptor_list, smoothing_zone_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("smoothMesh", std::function(
+
+                                            [](std::shared_ptr<const IMesh> p_mesh,
+                                               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                                 bc_descriptor_list,
+                                               const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>&
+                                                 discrete_function_variant_list) -> std::shared_ptr<const IMesh> {
+                                              MeshSmootherHandler handler;
+                                              return handler.getSmoothedMesh(*p_mesh, bc_descriptor_list,
+                                                                             discrete_function_variant_list);
+                                            }
+
+                                            ));
+
   this->_addBuiltinFunction("fixed", std::function(
 
                                        [](std::shared_ptr<const IBoundaryDescriptor> boundary)

src/mesh/MeshSmoother.cpp

@@ -6,10 +6,12 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 #include <scheme/AxisBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <utils/RandomEngine.hpp>
@@ -183,20 +185,11 @@ class MeshSmootherHandler::MeshSmoother
               const NodeId cell_node_id = cell_node_list[i_node];
               mean_position += given_xr[cell_node_id];
               number_of_neighbours++;
-              // std::cout << node_id << " position " << given_xr[node_id] << " position voisin " <<
-              // given_xr[cell_node_id]
-              //           << " mean_position " << mean_position << "\n";
             }
           }
         }
         mean_position    = 1. / number_of_neighbours * mean_position;
         shift_r[node_id] = mean_position - given_xr[node_id];
-        // std::cout << " mean position " << mean_position << given_xr[node_id] << "\n";
-
-        // double nshift    = sqrt(dot(shift_r[node_id], shift_r[node_id]));
-        // if (nshift > max_delta_xr[node_id]) {
-        //   shift_r[node_id] = max_delta_xr[node_id] / nshift * shift_r[node_id];
-        // }
       });
 
     this->_applyBC(shift_r);
@@ -237,6 +230,83 @@ class MeshSmootherHandler::MeshSmoother
     return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
   }
 
+  std::shared_ptr<const IMesh>
+  getSmoothedMesh(const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_descriptor_list) const
+  {
+    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
+
+    NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+    is_displaced.fill(false);
+
+    for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) {
+      MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]);
+      const auto cell_list              = cell_zone.cellList();
+      CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
+      is_zone_cell.fill(false);
+      parallel_for(
+        cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { is_zone_cell[cell_list[i_cell]] = true; });
+      parallel_for(
+        m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+          auto node_cell_list = node_to_cell_matrix[node_id];
+          bool displace       = true;
+          for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+            const CellId cell_id = node_cell_list[i_node_cell];
+            displace &= is_zone_cell[cell_id];
+          }
+          if (displace) {
+            is_displaced[node_id] = true;
+          }
+        });
+    }
+    synchronize(is_displaced);
+    NodeValue<Rd> xr = this->_getDisplacement();
+
+    parallel_for(
+      m_given_mesh.numberOfNodes(),
+      PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+  }
+
+  std::shared_ptr<const IMesh>
+  getSmoothedMesh(
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+  {
+    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
+
+    NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+    is_displaced.fill(false);
+
+    for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) {
+      auto is_zone_cell =
+        discrete_function_variant_list[i_zone]->get<const DiscreteFunctionP0<Dimension, const double>>();
+      parallel_for(
+        m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+          auto node_cell_list = node_to_cell_matrix[node_id];
+          bool displace       = true;
+          for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+            const CellId cell_id = node_cell_list[i_node_cell];
+            displace &= (is_zone_cell[cell_id] != 0);
+          }
+          if (displace) {
+            is_displaced[node_id] = true;
+          }
+        });
+    }
+    synchronize(is_displaced);
+    NodeValue<Rd> xr = this->_getDisplacement();
+
+    parallel_for(
+      m_given_mesh.numberOfNodes(),
+      PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+  }
+
   MeshSmoother(const MeshSmoother&) = delete;
   MeshSmoother(MeshSmoother&&)      = delete;
 
@@ -396,3 +466,65 @@ MeshSmootherHandler::getSmoothedMesh(
   }
   }
 }
+
+std::shared_ptr<const IMesh>
+MeshSmootherHandler::getSmoothedMesh(
+  const IMesh& mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const
+{
+  switch (mesh.dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(smoothing_zone_list);
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(smoothing_zone_list);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(smoothing_zone_list);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IMesh>
+MeshSmootherHandler::getSmoothedMesh(
+  const IMesh& mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+  switch (mesh.dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(discrete_function_variant_list);
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(discrete_function_variant_list);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(discrete_function_variant_list);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}

src/mesh/MeshSmoother.hpp

@@ -8,6 +8,8 @@
 #include <vector>
 
 class FunctionSymbolId;
+class IZoneDescriptor;
+class DiscreteFunctionVariant;
 
 class MeshSmootherHandler
 {
@@ -25,6 +27,16 @@ class MeshSmootherHandler
     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
     const FunctionSymbolId& function_symbol_id) const;
 
+  std::shared_ptr<const IMesh> getSmoothedMesh(
+    const IMesh& mesh,
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+    const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const;
+
+  std::shared_ptr<const IMesh> getSmoothedMesh(
+    const IMesh& mesh,
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& smoothing_zone_list) const;
+
   MeshSmootherHandler()                      = default;
   MeshSmootherHandler(MeshSmootherHandler&&) = default;
   ~MeshSmootherHandler()                     = default;