Neumann for elasticity

src/language/algorithms/CMakeLists.txt

@@ -3,7 +3,6 @@
 add_library(PugsLanguageAlgorithms
   AcousticSolverAlgorithm.cpp
   ElasticityDiamondAlgorithm.cpp
-  ElasticityDiamondAlgorithm2.cpp
   Heat5PointsAlgorithm.cpp
   HeatDiamondAlgorithm.cpp
 )

src/language/algorithms/ElasticityDiamondAlgorithm.cpp

@@ -82,50 +82,17 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
           const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
           const RefId& ref          = ref_face_list.refId();
           if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
-            MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
-
             const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
 
-            const auto& node_list = mesh_node_boundary.nodeList();
-
-            Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-              TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(), node_list);
-
-            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-              NodeId node_id = node_list[i_node];
-              if (not is_dirichlet[node_id]) {
-                is_dirichlet[node_id]    = true;
-                dirichlet_value[node_id] = value_list[i_node];
-              }
-            }
-
-            boundary_condition_list.push_back(DirichletBoundaryCondition{node_list, value_list});
-          }
-        }
-
             if constexpr (Dimension > 1) {
-          for (size_t i_ref_node_list = 0;
-               i_ref_node_list < mesh->connectivity().template numberOfRefItemList<ItemType::node>();
-               ++i_ref_node_list) {
-            const auto& ref_node_list = mesh->connectivity().template refItemList<ItemType::node>(i_ref_node_list);
-            const RefId& ref          = ref_node_list.refId();
-            if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
-              const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
-
-              const auto& node_list = ref_node_list.list();
+              MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
 
+              Array<const FaceId> face_list                 = ref_face_list.list();
               Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-                TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(), node_list);
-
-              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-                NodeId node_id = node_list[i_node];
-                if (not is_dirichlet[node_id]) {
-                  is_dirichlet[node_id]    = true;
-                  dirichlet_value[node_id] = value_list[i_node];
-                }
-              }
-
-              boundary_condition_list.push_back(DirichletBoundaryCondition{node_list, value_list});
+                TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(), face_list);
+              boundary_condition_list.push_back(DirichletBoundaryCondition{face_list, value_list});
+            } else {
+              throw NotImplementedError("Neumann conditions are not supported in 1d");
             }
           }
         }
@@ -182,7 +149,8 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
           [&](auto&& bc) {
             using T = std::decay_t<decltype(bc)>;
             if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>) or
-                          (std::is_same_v<T, SymmetryBoundaryCondition>)) {
+                          (std::is_same_v<T, SymmetryBoundaryCondition>) or
+                          (std::is_same_v<T, DirichletBoundaryCondition>)) {
               const auto& face_list = bc.faceList();
 
               for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
@@ -205,6 +173,28 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
 
     const auto& primal_face_to_node_matrix             = mesh->connectivity().faceToNodeMatrix();
     const auto& face_to_cell_matrix                    = mesh->connectivity().faceToCellMatrix();
+    const FaceValue<const bool> primal_face_is_neumann = [&] {
+      FaceValue<bool> face_is_neumann{mesh->connectivity()};
+      face_is_neumann.fill(false);
+      for (const auto& boundary_condition : boundary_condition_list) {
+        std::visit(
+          [&](auto&& bc) {
+            using T = std::decay_t<decltype(bc)>;
+            if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
+              const auto& face_list = bc.faceList();
+
+              for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+                const FaceId face_id     = face_list[i_face];
+                face_is_neumann[face_id] = true;
+              }
+            }
+          },
+          boundary_condition);
+      }
+
+      return face_is_neumann;
+    }();
+
     NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
     if (parallel::size() > 1) {
       throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
@@ -232,21 +222,14 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
       }
     }
 
-    FaceValue<bool> primal_face_is_neumann(mesh->connectivity());
+    FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
     if (parallel::size() > 1) {
       throw NotImplementedError("Calculation of face_is_neumann is incorrect");
     }
 
-    primal_face_is_neumann.fill(false);
+    primal_face_is_dirichlet.fill(false);
     for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-      if (primal_face_is_on_boundary[face_id]) {
-        for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-          NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
-          if (not is_dirichlet[node_id]) {
-            primal_face_is_neumann[face_id] = true;
-          }
-        }
-      }
+      primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
     }
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
 
@@ -447,10 +430,10 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
         return computed_dual_node_primal_node_id;
       }();
 
-      TinyMatrix<Dimension, double> eye = zero;
-      for (size_t i = 0; i < Dimension; ++i) {
-        eye(i, i) = 1;
-      }
+      TinyMatrix<Dimension, double> eye = identity;
+      // for (size_t i = 0; i < Dimension; ++i) {
+      //    eye(i, i) = 1;
+      // }
 
       const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
       const auto& primal_face_cell_is_reversed                       = mesh->connectivity().cellFaceIsReversed();
@@ -479,6 +462,9 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
               for (size_t i = 0; i < Dimension; ++i) {
                 for (size_t j = 0; j < Dimension; ++j) {
                   S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) += M(i, j);
+                  if (primal_face_is_neumann[face_id]) {
+                    S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) -= M(i, j);
+                  }
                 }
               }
             } else {
@@ -507,7 +493,7 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
           for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
             NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
 
-            if (not primal_node_is_on_boundary[node_id]) {
+            //            if (not primal_node_is_on_boundary[node_id]) {
             const TinyVector<Dimension> nil = [&] {
               if (is_face_reversed_for_cell_i) {
                 return -primal_nlr(face_id, i_node);
@@ -533,6 +519,29 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
                     for (size_t j = 0; j < Dimension; ++j) {
                       S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -=
                         w_rj(node_id, j_cell) * M(i, j);
+                      if (primal_face_is_neumann[face_id]) {
+                        S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
+                          w_rj(node_id, j_cell) * M(i, j);
+                      }
+                    }
+                  }
+                }
+                if (primal_node_is_on_boundary[node_id]) {
+                  for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
+                    FaceId l_id = node_to_face_matrix[node_id][l_face];
+                    if (primal_face_is_on_boundary[l_id]) {
+                      for (size_t i = 0; i < Dimension; ++i) {
+                        for (size_t j = 0; j < Dimension; ++j) {
+                          S(cell_dof_number[i_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) -=
+                            w_rl(node_id, l_face) * M(i, j);
+                        }
+                      }
+                      if (primal_face_is_neumann[face_id]) {
+                        for (size_t i = 0; i < Dimension; ++i) {
+                          for (size_t j = 0; j < Dimension; ++j) {
+                            S(face_dof_number[face_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) +=
+                              w_rl(node_id, l_face) * M(i, j);
+                          }
                         }
                       }
                     }
@@ -540,6 +549,15 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
                 }
               }
             }
+            //            }
+          }
+        }
+      }
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        if (primal_face_is_dirichlet[face_id]) {
+          for (size_t i = 0; i < Dimension; ++i) {
+            S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + i) += 1;
+          }
         }
       }
 
@@ -551,7 +569,6 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
         }
       }
 
-      const auto& primal_cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
       // Dirichlet
       NodeValue<bool> node_tag{mesh->connectivity()};
       node_tag.fill(false);
@@ -560,102 +577,109 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
           [&](auto&& bc) {
             using T = std::decay_t<decltype(bc)>;
             if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
-              const auto& node_list  = bc.nodeList();
+              const auto& face_list  = bc.faceList();
               const auto& value_list = bc.valueList();
-              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-                const NodeId node_id = node_list[i_node];
-                if (not node_tag[node_id]) {
-                  node_tag[node_id] = true;
-
-                  const auto& node_cells = primal_node_to_cell_matrix[node_id];
-                  for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
-                    const CellId cell_id = node_cells[i_cell];
-
-                    const auto& primal_cell_to_face = primal_cell_to_face_matrix[cell_id];
-
-                    for (size_t i_cell_face = 0; i_cell_face < primal_cell_to_face.size(); ++i_cell_face) {
-                      FaceId face_id = primal_cell_to_face_matrix[cell_id][i_cell_face];
+              for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+                const FaceId face_id = face_list[i_face];
 
-                      const bool is_face_reversed_for_cell_i = [=] {
-                        for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) {
-                          FaceId primal_face_id = primal_cell_to_face[i_face];
-                          if (primal_face_id == face_id) {
-                            return primal_face_cell_is_reversed(cell_id, i_face);
+                for (size_t i = 0; i < Dimension; ++i) {
+                  b[(face_dof_number[face_id] * Dimension) + i] += value_list[i_face][i];
                 }
               }
-                        Assert(false, "cannot find cell's face");
-                        return false;
-                      }();
-
-                      const auto& primal_face_to_node = primal_face_to_node_matrix[face_id];
-
-                      for (size_t i_face_node = 0; i_face_node < primal_face_to_node.size(); ++i_face_node) {
-                        if (primal_face_to_node[i_face_node] == node_id) {
-                          const TinyVector<Dimension> nil = [=] {
-                            if (is_face_reversed_for_cell_i) {
-                              return -primal_nlr(face_id, i_face_node);
-                            } else {
-                              return primal_nlr(face_id, i_face_node);
             }
-                          }();
-
-                          CellId dual_cell_id = face_dual_cell_id[face_id];
+          },
+          boundary_condition);
+      }
 
-                          for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
-                               ++i_dual_node) {
-                            const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                            if (dual_node_primal_node_id[dual_node_id] == node_id) {
-                              const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-                              TinyMatrix<Dimension, double> M = alpha_mu_l[face_id] * (Clr, nil) * eye +
-                                                                alpha_mu_l[face_id] * tensorProduct(Clr, nil) +
-                                                                alpha_lambda_l[face_id] * tensorProduct(nil, Clr);
+      for (const auto& boundary_condition : boundary_condition_list) {
+        std::visit(
+          [&](auto&& bc) {
+            using T = std::decay_t<decltype(bc)>;
+            if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
+              const auto& face_list  = bc.faceList();
+              const auto& value_list = bc.valueList();
+              for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+                FaceId face_id = face_list[i_face];
+                Assert(face_to_cell_matrix[face_id].size() == 1);
                 for (size_t i = 0; i < Dimension; ++i) {
-                                b[(cell_dof_number[cell_id] * Dimension) + i] += (M * value_list[i_node])[i];
-                              }
-                            }
+                  b[face_dof_number[face_id] * Dimension + i] += mes_l[face_id] * value_list[i_face][i];   // sign
                 }
               }
             }
+          },
+          boundary_condition);
       }
+
+      Vector<double> U{number_of_dof * Dimension};
+      U = 1;
+      CellValue<TinyVector<Dimension>> Speed{mesh->connectivity()};
+      FaceValue<TinyVector<Dimension>> Ul = InterpolateItemValue<TinyVector<Dimension>(
+        TinyVector<Dimension>)>::template interpolate<ItemType::face>(U_id, mesh_data.xl());
+      FaceValue<TinyVector<Dimension>> Speed_face{mesh->connectivity()};
+
+      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+        for (size_t i = 0; i < Dimension; ++i) {
+          U[(cell_dof_number[cell_id] * Dimension) + i] = Uj[cell_id][i];
         }
       }
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        for (size_t i = 0; i < Dimension; ++i) {
+          if (primal_face_is_on_boundary[face_id]) {
+            U[(face_dof_number[face_id] * Dimension) + i] = Ul[face_id][i];
           }
         }
-          },
-          boundary_condition);
       }
+      Vector r = A * U - b;
 
-      Vector<double> U{number_of_dof * Dimension};
-      U = 0;
-
+      std::cout << "initial residu = " << std::sqrt((r, r)) << '\n';
       BiCGStab{b, A, U, 10000, 1e-15};
 
-      Vector r = A * U - b;
+      r = A * U - b;
 
       std::cout << "real residu = " << std::sqrt((r, r)) << '\n';
 
-      CellValue<TinyVector<Dimension>> Speed{mesh->connectivity()};
-
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
         for (size_t i = 0; i < Dimension; ++i) {
           Speed[cell_id][i] = U[(cell_dof_number[cell_id] * Dimension) + i];
         }
       }
-      Vector<double> Uexacte{number_of_dof * Dimension};
-      for (CellId j = 0; j < number_of_dof; ++j) {
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        for (size_t i = 0; i < Dimension; ++i) {
+          if (primal_face_is_on_boundary[face_id]) {
+            Speed_face[face_id][i] = U[(face_dof_number[face_id] * Dimension) + i];
+          } else {
+            Speed_face[face_id][i] = Ul[face_id][i];
+          }
+        }
+      }
+      Vector<double> Uexacte{mesh->numberOfCells() * Dimension};
+      for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
         for (size_t l = 0; l < Dimension; ++l) {
           Uexacte[(cell_dof_number[j] * Dimension) + l] = Uj[j][l];
         }
       }
 
-      Vector<double> error{number_of_dof * Dimension};
+      Vector<double> error{mesh->numberOfCells() * Dimension};
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
         for (size_t i = 0; i < Dimension; ++i) {
           error[(cell_id * Dimension) + i] = (Speed[cell_id][i] - Uj[cell_id][i]) * sqrt(primal_Vj[cell_id]);
         }
       }
+      Vector<double> error_face{mesh->numberOfFaces() * Dimension};
+      parallel_for(
+        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+          if (primal_face_is_on_boundary[face_id]) {
+            for (size_t i = 0; i < Dimension; ++i) {
+              error_face[face_id * Dimension + i] = (Speed_face[face_id][i] - Ul[face_id][i]) * sqrt(mes_l[face_id]);
+            }
+          } else {
+            error_face[face_id] = 0;
+          }
+        });
 
-      std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n";
+      std::cout << "||Error||_2 (cell)= " << std::sqrt((error, error)) << "\n";
+      std::cout << "||Error||_2 (face)= " << std::sqrt((error_face, error_face)) << "\n";
+      std::cout << "||Error||_2 (total)= " << std::sqrt((error, error)) + std::sqrt((error_face, error_face)) << "\n";
 
       VTKWriter vtk_writer("Speed_" + std::to_string(Dimension), 0.01);
 
@@ -672,13 +696,13 @@ class ElasticityDiamondScheme<Dimension>::DirichletBoundaryCondition
 {
  private:
   const Array<const TinyVector<Dimension>> m_value_list;
-  const Array<const NodeId> m_node_list;
+  const Array<const FaceId> m_face_list;
 
  public:
-  const Array<const NodeId>&
-  nodeList() const
+  const Array<const FaceId>&
+  faceList() const
   {
-    return m_node_list;
+    return m_face_list;
   }
 
   const Array<const TinyVector<Dimension>>&
@@ -687,10 +711,10 @@ class ElasticityDiamondScheme<Dimension>::DirichletBoundaryCondition
     return m_value_list;
   }
 
-  DirichletBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
-    : m_value_list{value_list}, m_node_list{node_list}
+  DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const TinyVector<Dimension>>& value_list)
+    : m_value_list{value_list}, m_face_list{face_list}
   {
-    Assert(m_value_list.size() == m_node_list.size());
+    Assert(m_value_list.size() == m_face_list.size());
   }
 
   ~DirichletBoundaryCondition() = default;

src/language/algorithms/ElasticityDiamondAlgorithm.hpp

@@ -26,4 +26,4 @@ class ElasticityDiamondScheme
                           const FunctionSymbolId& U_id);
 };
 
-#endif   // ELASTICITY_DIAMOND_ALGORITHM_HPP
+#endif   // ELASTICITY_DIAMOND_ALGORITHM2_HPP

src/language/modules/SchemeModule.cpp

@@ -2,7 +2,6 @@
 
 #include <language/algorithms/AcousticSolverAlgorithm.hpp>
 #include <language/algorithms/ElasticityDiamondAlgorithm.hpp>
-#include <language/algorithms/ElasticityDiamondAlgorithm2.hpp>
 #include <language/algorithms/Heat5PointsAlgorithm.hpp>
 #include <language/algorithms/HeatDiamondAlgorithm.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
@@ -294,39 +293,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("elasticity2",
-                            std::make_shared<BuiltinFunctionEmbedder<
-                              void(std::shared_ptr<const IMesh>,
-                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&,
-                                   const FunctionSymbolId&)>>(
-
-                              [](std::shared_ptr<const IMesh> p_mesh,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const FunctionSymbolId& lambda_id, const FunctionSymbolId& mu_id,
-                                 const FunctionSymbolId& f_id, const FunctionSymbolId& U_id) -> void {
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  ElasticityDiamondScheme2<1>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
-                                  break;
-                                }
-                                case 2: {
-                                  ElasticityDiamondScheme2<2>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
-                                  break;
-                                }
-                                case 3: {
-                                  ElasticityDiamondScheme2<3>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
-                                  break;
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
-                                }
-                              }
-
-                              ));
-
   this->_addBuiltinFunction("heat5points",
                             std::make_shared<BuiltinFunctionEmbedder<
                               void(std::shared_ptr<const IMesh>,