add symmetry BCs to the implicit smoother

src/mesh/ImplicitMeshSmoother.cpp

@@ -79,27 +79,31 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
   }
 
   void
-  _browseBC(NodeValue<bool>& is_fixed) const
+  _browseBC(NodeValue<bool>& is_fixed, NodeValue<int>& is_symmetric) const
   {
+    int n_sym = 1;
     for (auto&& boundary_condition : m_boundary_condition_list) {
       std::visit(
         [&](auto&& bc) {
           using BCType = std::decay_t<decltype(bc)>;
-          // if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+          if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
             //   const Rd& n = bc.outgoingNormal();
 
             //   const Rdxd I   = identity;
             //   const Rdxd nxn = tensorProduct(n, n);
             //   const Rdxd P   = I - nxn;
 
-          //   const Array<const NodeId>& node_list = bc.nodeList();
-          //   parallel_for(
-          //     node_list.size(), PUGS_LAMBDA(const size_t i_node) {
-          //       const NodeId node_id = node_list[i_node];
-
-          //       shift[node_id] = P * shift[node_id];
-          //     });
-
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+                const NodeId node_id = node_list[i_node];
+                if (is_symmetric[node_id] != 0) {
+                  is_fixed[node_id] = true;
+                } else if (!is_fixed[node_id]) {
+                  is_symmetric[node_id] = n_sym;
+                }
+              });
+            n_sym += 1;
             // } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
             //   if constexpr (Dimension > 1) {
             //     const Rd& t = bc.direction();
@@ -117,8 +121,7 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
             //     throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
             //   }
 
-          // } else
-          if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
+          } else if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
             const Array<const NodeId>& node_list = bc.nodeList();
             parallel_for(
               node_list.size(), PUGS_LAMBDA(const size_t i_node) {
@@ -167,22 +170,23 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
     int local_free_id  = 0;
     int local_fixed_id = 0;
     non_zeros_free.fill(1);
+    // Warning this is not correct
     non_zeros_fixed.fill(0);
     for (NodeId n_id = 0; n_id < m_given_mesh.numberOfNodes(); n_id++) {
-      // std::cout << " n_id " << n_id << "\n";
       if (is_fixed[n_id]) {
         gid_node_fixed[local_fixed_id] = n_id;
-        for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
-          FaceId face_id = node_to_face_matrix[n_id][i_face];
-          for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
-            NodeId node_id = face_to_node_matrix[face_id][i_node];
-            if ((node_id == n_id) or (is_fixed[node_id])) {
-              continue;
-            } else {
-              non_zeros_fixed[local_fixed_id] += 1;
-            }
-          }
-        }
+        // for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
+        //   FaceId face_id = node_to_face_matrix[n_id][i_face];
+
+        //   for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
+        //     NodeId node_id = face_to_node_matrix[face_id][i_node];
+        //     if ((node_id == n_id) or (is_fixed[node_id])) {
+        //       continue;
+        //     } else {
+        //       non_zeros_fixed[local_fixed_id] += 1;
+        //     }
+        //   }
+        // }
         local_fixed_id++;
       } else {
         gid_node_free[local_free_id] = n_id;
@@ -190,8 +194,10 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
           FaceId face_id = node_to_face_matrix[n_id][i_face];
           for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
             NodeId node_id = face_to_node_matrix[face_id][i_node];
-            if ((node_id == n_id) or (is_fixed[node_id])) {
+            if ((node_id == n_id)) {
               continue;
+            } else if ((is_fixed[node_id])) {
+              non_zeros_fixed[local_free_id] += 1;
             } else {
               non_zeros_free[local_free_id] += 1;
             }
@@ -205,6 +211,7 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
   _fillMatrix(CRSMatrixDescriptor<double>& Afree,
               CRSMatrixDescriptor<double>& Afixed,
               const NodeValue<bool>& is_fixed,
+              const NodeValue<int>& is_symmetric,
               const NodeValue<int>& id_free,
               const NodeValue<int>& id_fixed) const
   // ,
@@ -219,6 +226,24 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
         continue;
       } else {
         Assert(id_free[n_id] == local_free_id, "bad matrix definition");
+        if (is_symmetric[n_id]) {
+          for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
+            FaceId face_id = node_to_face_matrix[n_id][i_face];
+            for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
+              NodeId node_id = face_to_node_matrix[face_id][i_node];
+              if (node_id == n_id) {
+                continue;
+              } else if ((is_fixed[node_id])) {
+                Afree(local_free_id, local_free_id) += 1;
+                Afixed(local_free_id, id_fixed[node_id]) -= 1;
+              } else if (is_symmetric[node_id] == is_symmetric[n_id]) {
+                Afree(local_free_id, local_free_id) += 1;
+                Afree(local_free_id, id_free[node_id]) -= 1;
+              }
+            }
+          }
+
+        } else {
           for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
             FaceId face_id = node_to_face_matrix[n_id][i_face];
             for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
@@ -234,11 +259,55 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
               }
             }
           }
+        }
         local_free_id++;
       }
     }
   }
 
+  void
+  _correct_sym(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
+  {
+    for (auto&& boundary_condition : m_boundary_condition_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using BCType = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+            const Rd& n = bc.outgoingNormal();
+
+            const Rdxd I   = identity;
+            const Rdxd nxn = tensorProduct(n, n);
+            const Rdxd P   = I - nxn;
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+                const NodeId node_id = node_list[i_node];
+                Rd depl              = new_xr[node_id] - old_xr[node_id];
+                new_xr[node_id]      = old_xr[node_id] + P * depl;
+              });
+            // } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
+            //   if constexpr (Dimension > 1) {
+            //     const Rd& t = bc.direction();
+
+            //     const Rdxd txt = tensorProduct(t, t);
+
+            //     const Array<const NodeId>& node_list = bc.nodeList();
+            //     parallel_for(
+            //       node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+            //         const NodeId node_id = node_list[i_node];
+
+            //         shift[node_id] = txt * shift[node_id];
+            //       });
+            //   } else {
+            //     throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
+            //   }
+          }
+        },
+        boundary_condition);
+    }
+  }
+
   NodeValue<Rd>
   _getPosition(NodeValue<const bool> is_displaced) const
   {
@@ -247,7 +316,9 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
 
     NodeValue<Rd> pos_r{connectivity};
     NodeValue<bool> is_fixed{connectivity};
+    NodeValue<int> is_symmetric{connectivity};
     is_fixed.fill(false);
+    is_symmetric.fill(0);
     parallel_for(
       m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
         pos_r[node_id] = given_xr[node_id];
@@ -255,24 +326,26 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
           is_fixed[node_id] = true;
         }
       });
-    _browseBC(is_fixed);
 
+    _browseBC(is_fixed, is_symmetric);
+    // std::cout << "is_fixed " << is_fixed << "\n";
+    // std::cout << "is_symmetric " << is_symmetric << "\n";
     NodeValue<int> node_dof_id{connectivity};
     size_t nb_free, nb_fixed;
     NodeValue<int> id_free{connectivity};
     NodeValue<int> id_fixed{connectivity};
     _computeMatrixSize(is_fixed, id_fixed, id_free, nb_free, nb_fixed);
+    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};
     _findLocalIds(is_fixed, gid_node_free, gid_node_fixed, non_zeros_free, non_zeros_fixed);
-    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);
+    _fillMatrix(Afree, Afixed, is_fixed, is_symmetric, id_free, id_fixed);
     Vector<double> F{nb_fixed};
     CRSMatrix Mfree{Afree.getCRSMatrix()};
     CRSMatrix Mfixed{Afixed.getCRSMatrix()};
@@ -283,6 +356,10 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
       Vector<double> X{nb_free};
       Vector<double> b{nb_free};
       b = Mfixed * F;
+      // std::cout << " Mfixed " << Mfixed << "\n";
+      // std::cout << " F " << F << "\n";
+      // std::cout << " b " << b << "\n";
+      // std::cout << " Mfree " << Mfree << "\n";
       solver.solveLocalSystem(Mfree, X, b);
       parallel_for(
         m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
@@ -291,6 +368,7 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
           }
         });
     }
+    _correct_sym(m_given_mesh.xr(), pos_r);
     //    synchronize(pos_r);
     return pos_r;
   }