before implementing symmetry

src/language/modules/SchemeModule.cpp

@@ -268,6 +268,34 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("implicitSmoothMesh",
+                            std::function(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const FunctionSymbolId& function_symbol_id) -> std::shared_ptr<const IMesh> {
+                                ImplicitMeshSmootherHandler handler;
+                                return handler.getSmoothedMesh(p_mesh, bc_descriptor_list, function_symbol_id);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("implicitSmoothMesh",
+                            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> {
+                                ImplicitMeshSmootherHandler handler;
+                                return handler.getSmoothedMesh(p_mesh, bc_descriptor_list,
+                                                               discrete_function_variant_list);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("smoothMesh", std::function(
 
                                             [](std::shared_ptr<const IMesh> p_mesh,

src/mesh/ImplicitMeshSmoother.cpp

@@ -32,10 +32,10 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
 
   //  class AxisBoundaryCondition;
   class FixedBoundaryCondition;
-  //  class SymmetryBoundaryCondition;
+  class SymmetryBoundaryCondition;
 
   //  using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
-  using BoundaryCondition = std::variant<FixedBoundaryCondition>;
+  using BoundaryCondition = std::variant<FixedBoundaryCondition, SymmetryBoundaryCondition>;
 
   using BoundaryConditionList = std::vector<BoundaryCondition>;
   BoundaryConditionList m_boundary_condition_list;
@@ -58,11 +58,11 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
         //   }
         //   break;
         // }
-      // case IBoundaryConditionDescriptor::Type::symmetry: {
-      //   bc_list.emplace_back(
-      //     SymmetryBoundaryCondition{getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
-      //   break;
-      // }
+      case IBoundaryConditionDescriptor::Type::symmetry: {
+        bc_list.emplace_back(
+          SymmetryBoundaryCondition{getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        break;
+      }
       case IBoundaryConditionDescriptor::Type::fixed: {
         bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
         break;
@@ -240,7 +240,7 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
   }
 
   NodeValue<Rd>
-  _getPosition() const
+  _getPosition(NodeValue<const bool> is_displaced) const
   {
     const ConnectivityType& connectivity = m_given_mesh.connectivity();
     NodeValue<const Rd> given_xr         = m_given_mesh.xr();
@@ -249,7 +249,12 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
     NodeValue<bool> is_fixed{connectivity};
     is_fixed.fill(false);
     parallel_for(
-      m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { pos_r[node_id] = given_xr[node_id]; });
+      m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+        pos_r[node_id] = given_xr[node_id];
+        if (not is_displaced[node_id]) {
+          is_fixed[node_id] = true;
+        }
+      });
     _browseBC(is_fixed);
 
     NodeValue<int> node_dof_id{connectivity};
@@ -257,51 +262,17 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
     NodeValue<int> id_free{connectivity};
     NodeValue<int> id_fixed{connectivity};
     _computeMatrixSize(is_fixed, id_fixed, id_free, nb_free, nb_fixed);
-    // std::cout << " is_fixed: "
-    //           << "\n";
-    // std::cout << is_fixed << "\n";
-    // std::cout << " id_fixed: "
-    //           << "\n";
-    // std::cout << id_fixed << "\n";
-    // std::cout << " is_free: "
-    //           << "\n";
-    // std::cout << id_free << "\n";
-    //    std::cout << "nb_free " << nb_free << " nb_fixed " << nb_fixed << "\n";
     Array<int> non_zeros_free{nb_free};
     Array<int> non_zeros_fixed{nb_free};
     Array<NodeId> gid_node_free{nb_free};
     Array<NodeId> gid_node_fixed{nb_fixed};
-    // size_t local_free_id  = 0;
-    // size_t local_fixed_id = 0;
     _findLocalIds(is_fixed, gid_node_free, gid_node_fixed, non_zeros_free, non_zeros_fixed);
-    // for (NodeId n_id = 0; n_id < m_given_mesh.numberOfNodes(); n_id++) {
-    //   if (is_fixed[n_id]) {
-    //     gid_node_fixed[local_fixed_id] = n_id;
-    //     local_fixed_id++;
-    //   } else {
-    //     gid_node_free[local_free_id] = n_id;
-    //     local_free_id++;
-    //   }
-    // }
-    // std::cout << " gid_node_fixed: "
-    //           << "\n";
-    // std::cout << gid_node_fixed << "\n";
-    // std::cout << " gid_node_free: "
-    //           << "\n";
-    // std::cout << gid_node_free << "\n";
     std::cout << " nb_free " << nb_free << " nb_fixed " << nb_fixed << "\n";
 
     CRSMatrixDescriptor<double> Afree(nb_free, nb_free, non_zeros_free);
     CRSMatrixDescriptor<double> Afixed(nb_free, nb_fixed, non_zeros_fixed);
     LinearSolver solver;
     _fillMatrix(Afree, Afixed, is_fixed, id_free, id_fixed);
-    // std::cout << " Afree :"
-    //           << "\n";
-    // std::cout << Afree << "\n";
-    // std::cout << " Afixed :"
-    //           << "\n";
-    // std::cout << Afixed << "\n";
-
     Vector<double> F{nb_fixed};
     CRSMatrix Mfree{Afree.getCRSMatrix()};
     CRSMatrix Mfixed{Afixed.getCRSMatrix()};
@@ -312,8 +283,6 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
       Vector<double> X{nb_free};
       Vector<double> b{nb_free};
       b = Mfixed * F;
-      // std::cout << " F " << F << "\n";
-      // std::cout << " b " << b << "\n";
       solver.solveLocalSystem(Mfree, X, b);
       parallel_for(
         m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
@@ -322,16 +291,62 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
           }
         });
     }
+    //    synchronize(pos_r);
     return pos_r;
   }
 
  public:
   std::shared_ptr<const IMesh>
   getSmoothedMesh() const
+  {
+    NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+    is_displaced.fill(true);
+    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    NodeValue<Rd> xr = this->_getPosition(is_displaced);
+
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+  }
+
+  std::shared_ptr<const IMesh>
+  getSmoothedMesh(const FunctionSymbolId& function_symbol_id) const
   {
     NodeValue<const Rd> given_xr = m_given_mesh.xr();
+    NodeValue<const bool> is_displaced =
+      InterpolateItemValue<bool(const Rd)>::interpolate(function_symbol_id, given_xr);
 
-    NodeValue<Rd> xr = this->_getPosition();
+    NodeValue<Rd> xr = this->_getPosition(is_displaced);
+
+    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<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;
+          }
+        });
+    }
+
+    NodeValue<Rd> xr = this->_getPosition(is_displaced);
 
     return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
   }
@@ -347,45 +362,6 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
   ~ImplicitMeshSmoother() = default;
 };
 
-// template <size_t Dimension>
-// class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<Dimension>::AxisBoundaryCondition
-// {
-//  public:
-//   using Rd = TinyVector<Dimension, double>;
-
-//  private:
-//   const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
-
-//  public:
-//   const Rd&
-//   direction() const
-//   {
-//     return m_mesh_line_node_boundary.direction();
-//   }
-
-//   const Array<const NodeId>&
-//   nodeList() const
-//   {
-//     return m_mesh_line_node_boundary.nodeList();
-//   }
-
-//   AxisBoundaryCondition(MeshLineNodeBoundary<Dimension>&& mesh_line_node_boundary)
-//     : m_mesh_line_node_boundary(mesh_line_node_boundary)
-//   {
-//     ;
-//   }
-
-//   ~AxisBoundaryCondition() = default;
-// };
-
-// template <>
-// class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<1>::AxisBoundaryCondition
-// {
-//  public:
-//   AxisBoundaryCondition()  = default;
-//   ~AxisBoundaryCondition() = default;
-// };
-
 template <size_t Dimension>
 class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<Dimension>::FixedBoundaryCondition
 {
@@ -404,36 +380,36 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<Dimension>::FixedBoundar
   ~FixedBoundaryCondition() = default;
 };
 
-// template <size_t Dimension>
-// class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<Dimension>::SymmetryBoundaryCondition
-// {
-//  public:
-//   using Rd = TinyVector<Dimension, double>;
+template <size_t Dimension>
+class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<Dimension>::SymmetryBoundaryCondition
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
 
-//  private:
-//   const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+ private:
+  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
 
-//  public:
-//   const Rd&
-//   outgoingNormal() const
-//   {
-//     return m_mesh_flat_node_boundary.outgoingNormal();
-//   }
+ public:
+  const Rd&
+  outgoingNormal() const
+  {
+    return m_mesh_flat_node_boundary.outgoingNormal();
+  }
 
-//   const Array<const NodeId>&
-//   nodeList() const
-//   {
-//     return m_mesh_flat_node_boundary.nodeList();
-//   }
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_flat_node_boundary.nodeList();
+  }
 
-//   SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
-//     : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
-//   {
-//     ;
-//   }
+  SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
+    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+  {
+    ;
+  }
 
-//   ~SymmetryBoundaryCondition() = default;
-// };
+  ~SymmetryBoundaryCondition() = default;
+};
 
 std::shared_ptr<const IMesh>
 ImplicitMeshSmootherHandler::getSmoothedMesh(
@@ -462,3 +438,70 @@ ImplicitMeshSmootherHandler::getSmoothedMesh(
   }
   }
 }
+std::shared_ptr<const IMesh>
+ImplicitMeshSmootherHandler::getSmoothedMesh(
+  const std::shared_ptr<const IMesh>& mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const FunctionSymbolId& function_symbol_id) const
+{
+  switch (mesh->dimension()) {
+  case 1: {
+    throw NotImplementedError("ImplicitMeshSmoother not implemented in 1D");
+    break;
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    ImplicitMeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(function_symbol_id);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    ImplicitMeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(function_symbol_id);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IMesh>
+ImplicitMeshSmootherHandler::getSmoothedMesh(
+  const std::shared_ptr<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
+{
+  if (not hasSameMesh(discrete_function_variant_list)) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  std::shared_ptr<const IMesh> common_mesh = getCommonMesh(discrete_function_variant_list);
+
+  if (common_mesh != mesh) {
+    throw NormalError("discrete functions are not defined on the smoothed mesh");
+  }
+
+  switch (mesh->dimension()) {
+  case 1: {
+    throw NotImplementedError("ImplicitMeshSmoother not implemented in 1D");
+    break;
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    ImplicitMeshSmoother 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>>;
+    ImplicitMeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(discrete_function_variant_list);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}

src/mesh/ImplicitMeshSmoother.hpp

@@ -41,6 +41,14 @@ class ImplicitMeshSmootherHandler
   std::shared_ptr<const IMesh> getSmoothedMesh(
     const std::shared_ptr<const IMesh>& mesh,
     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const;
+  std::shared_ptr<const IMesh> getSmoothedMesh(
+    const std::shared_ptr<const IMesh>& mesh,
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+    const FunctionSymbolId& function_symbol_id) const;
+  std::shared_ptr<const IMesh> getSmoothedMesh(
+    const std::shared_ptr<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;
 
   // std::shared_ptr<const IMesh> getSmoothedMesh(
   //   const std::shared_ptr<const IMesh>& mesh,