plumbing for symmetry boundary conditions for elasticity

src/language/algorithms/ElasticityDiamondAlgorithm.cpp

@@ -68,11 +68,25 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
 
     switch (bc_descriptor->type()) {
     case IBoundaryConditionDescriptor::Type::symmetry: {
-      // const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-      //   dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-      throw NotImplementedError("NIY");
+      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) {
+          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()) {
+            if constexpr (Dimension > 1) {
+              Array<const FaceId> face_list = ref_face_list.list();
+              boundary_condition_list.push_back(SymmetryBoundaryCondition{face_list});
+            } else {
+              throw NotImplementedError("Symmetry conditions are not supported in 1d");
               break;
             }
+          }
+        }
+      }
+    }
     case IBoundaryConditionDescriptor::Type::dirichlet: {
       const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
         dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
@@ -195,6 +209,28 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
       return face_is_neumann;
     }();
 
+    const FaceValue<const bool> primal_face_is_symmetry = [&] {
+      FaceValue<bool> face_is_symmetry{mesh->connectivity()};
+      face_is_symmetry.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, SymmetryBoundaryCondition>)) {
+              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_symmetry[face_id] = true;
+              }
+            }
+          },
+          boundary_condition);
+      }
+
+      return face_is_symmetry;
+    }();
+
     NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
     if (parallel::size() > 1) {
       throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
@@ -754,6 +790,7 @@ template <size_t Dimension>
 class ElasticityDiamondScheme<Dimension>::SymmetryBoundaryCondition
 {
  private:
+  const Array<const TinyVector<Dimension>> m_value_list;
   const Array<const FaceId> m_face_list;
 
  public:

src/scheme/SymmetryBoundaryConditionDescriptor.hpp

@@ -16,9 +16,16 @@ class SymmetryBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
     return os;
   }
 
+  const std::string_view m_name;
+
   std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
 
  public:
+  std::string_view
+  name() const
+  {
+    return m_name;
+  }
   const IBoundaryDescriptor&
   boundaryDescriptor() const
   {