symmetry for elastoplacity

src/language/algorithms/ElasticityDiamondAlgorithm.cpp

@@ -70,9 +70,8 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
     case IBoundaryConditionDescriptor::Type::symmetry: {
       const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
         dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-      if (sym_bc_descriptor.name() == "symmetry") {
-        for (size_t i_ref_face_list = 0;
-             i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+      for (size_t i_ref_face_list = 0; i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>();
+           ++i_ref_face_list) {
         const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
         const RefId& ref          = ref_face_list.refId();
         if (ref == sym_bc_descriptor.boundaryDescriptor()) {
@@ -81,11 +80,10 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
             boundary_condition_list.push_back(SymmetryBoundaryCondition{face_list});
           } else {
             throw NotImplementedError("Symmetry conditions are not supported in 1d");
-              break;
-            }
           }
         }
       }
+      break;
     }
     case IBoundaryConditionDescriptor::Type::dirichlet: {
       const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
@@ -265,10 +263,10 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
 
     primal_face_is_dirichlet.fill(false);
     for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-      primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
+      primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]) &&
+                                           (!primal_face_is_symmetry[face_id]));
     }
     MeshDataType& mesh_data                          = MeshDataManager::instance().getMeshData(*mesh);
-
     const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
 
     const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
@@ -278,6 +276,14 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
     CellValuePerNode<double> w_rj{mesh->connectivity()};
     FaceValuePerNode<double> w_rl{mesh->connectivity()};
 
+    const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
+    auto project_to_face = [&](const TinyVector<Dimension>& x, const FaceId face_id) -> const TinyVector<Dimension> {
+      TinyVector<Dimension> proj;
+      const TinyVector<Dimension> nil = primal_nlr(face_id, 0);
+      proj                            = x - ((x - xl[face_id]), nil) * nil;
+      return proj;
+    };
+
     for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
       w_rl[i] = std::numeric_limits<double>::signaling_NaN();
     }
@@ -334,7 +340,15 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
           for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
             FaceId face_id = node_to_face_matrix[i_node][i_face];
             if (primal_face_is_on_boundary[face_id]) {
+              if (primal_face_is_symmetry[face_id]) {
+                for (size_t j = 0; j < face_to_cell_matrix[face_id].size(); ++j) {
+                  const CellId cell_id                 = face_to_cell_matrix[face_id][j];
+                  TinyVector<Dimension> xproj          = project_to_face(xj[cell_id], face_id);
+                  A(i, node_to_cell.size() + cpt_face) = xproj[i - 1];
+                }
+              } else {
                 A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
+              }
               cpt_face++;
             }
           }
@@ -467,11 +481,7 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
       }();
 
       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();
       const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells();
       SparseMatrixDescriptor S{number_of_dof * Dimension};
@@ -494,6 +504,8 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
             const CellId j_id = primal_face_to_cell[j_cell];
             TinyMatrix<Dimension, double> M =
               beta_mu_l * eye + beta_mu_l * tensorProduct(nil, nil) + beta_lambda_l * tensorProduct(nil, nil);
+            TinyMatrix<Dimension, double> N = tensorProduct(nil, nil);
+
             if (i_cell == j_cell) {
               for (size_t i = 0; i < Dimension; ++i) {
                 for (size_t j = 0; j < Dimension; ++j) {
@@ -501,6 +513,12 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
                   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);
                   }
+                  if (primal_face_is_symmetry[face_id]) {
+                    S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
+                      -((i == j) ? 1 : 0) + N(i, j);
+                    S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + j) +=
+                      (i == j) ? 1 : 0;
+                  }
                 }
               }
             } else {
@@ -653,26 +671,13 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
         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];
-          }
-        }
-      }
       Vector r = A * U - b;
-
-      std::cout << "initial residu = " << std::sqrt((r, r)) << '\n';
+      std::cout << "initial (real) residu = " << std::sqrt((r, r)) << '\n';
       BiCGStab{b, A, U, 10000, 1e-15};
 
       r = A * U - b;
 
-      std::cout << "real residu = " << std::sqrt((r, r)) << '\n';
+      std::cout << "final (real) residu = " << std::sqrt((r, r)) << '\n';
 
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
         for (size_t i = 0; i < Dimension; ++i) {

src/language/algorithms/ElasticityDiamondAlgorithm.hpp

 #ifndef ELASTICITY_DIAMOND_ALGORITHM_HPP
 #define ELASTICITY_DIAMOND_ALGORITHM_HPP
 
+#include <algebra/TinyVector.hpp>
 #include <language/utils/FunctionSymbolId.hpp>
+#include <memory>
 #include <mesh/IMesh.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-
-#include <memory>
 #include <variant>
 #include <vector>
 

src/scheme/SymmetryBoundaryConditionDescriptor.hpp

@@ -12,7 +12,7 @@ class SymmetryBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
   std::ostream&
   _write(std::ostream& os) const final
   {
-    os << "symmetry(" << *m_boundary_descriptor << ")";
+    os << "symmetry(" << m_name << ',' << *m_boundary_descriptor << ")";
     return os;
   }