add normal stress BC

src/language/modules/SchemeModule.cpp

@@ -265,6 +265,18 @@ SchemeModule::SchemeModule()
 
                                           ));
 
+  this->_addBuiltinFunction("normal-stress",
+                            std::function(
+
+                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+                                 const FunctionSymbolId& normal_stress_id)
+                                -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                                return std::make_shared<DirichletBoundaryConditionDescriptor>("normal-stress", boundary,
+                                                                                              normal_stress_id);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("velocity", std::function(
 
                                           [](std::shared_ptr<const IBoundaryDescriptor> boundary,

src/scheme/HyperelasticSolver.cpp

@@ -82,8 +82,11 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
   class SymmetryBoundaryCondition;
   class VelocityBoundaryCondition;
 
-  using BoundaryCondition = std::
-    variant<FixedBoundaryCondition, PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
+  using BoundaryCondition = std::variant<FixedBoundaryCondition,
+                                         PressureBoundaryCondition,
+                                         NormalStressBoundaryCondition,
+                                         SymmetryBoundaryCondition,
+                                         VelocityBoundaryCondition>;
 
   using BoundaryConditionList = std::vector<BoundaryCondition>;
 
@@ -299,6 +302,29 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
             bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
           }
 
+        } else if (dirichlet_bc_descriptor.name() == "normal-stress") {
+          const FunctionSymbolId normal_stress_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+          if constexpr (Dimension == 1) {
+            MeshNodeBoundary<Dimension> mesh_node_boundary =
+              getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+            Array<const Rdxd> node_values =
+              InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::node>(normal_stress_id, mesh.xr(),
+                                                                                   mesh_node_boundary.nodeList());
+
+            bc_list.emplace_back(NormalStressBoundaryCondition{mesh_node_boundary, node_values});
+          } else {
+            MeshFaceBoundary<Dimension> mesh_face_boundary =
+              getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+            Array<const Rdxd> face_values =
+              InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::face>(normal_stress_id, mesh_data.xl(),
+                                                                                   mesh_face_boundary.faceList());
+            bc_list.emplace_back(NormalStressBoundaryCondition{mesh_face_boundary, face_values});
+          }
+
         } else {
           is_valid_boundary_condition = false;
         }
@@ -319,6 +345,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
   }
 
   void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
+  void _applyNormalStressBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
   void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
   void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
   void
@@ -328,6 +355,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
                            NodeValue<Rd>& br) const
   {
     this->_applyPressureBC(bc_list, mesh, br);
+    this->_applyNormalStressBC(bc_list, mesh, br);
     this->_applySymmetryBC(bc_list, Ar, br);
     this->_applyVelocityBC(bc_list, Ar, br);
   }
@@ -602,6 +630,71 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyPressureBC(const
   }
 }
 
+template <size_t Dimension>
+void
+HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyNormalStressBC(const BoundaryConditionList& bc_list,
+                                                                               const MeshType& mesh,
+                                                                               NodeValue<Rd>& br) const
+{
+  for (const auto& boundary_condition : bc_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using T = std::decay_t<decltype(bc)>;
+        if constexpr (std::is_same_v<NormalStressBoundaryCondition, T>) {
+          MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+          if constexpr (Dimension == 1) {
+            const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+            const auto& node_list  = bc.nodeList();
+            const auto& value_list = bc.valueList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(size_t i_node) {
+                const NodeId node_id       = node_list[i_node];
+                const auto& node_cell_list = node_to_cell_matrix[node_id];
+                Assert(node_cell_list.size() == 1);
+
+                CellId node_cell_id              = node_cell_list[0];
+                size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+
+                br[node_id] += value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
+              });
+          } else {
+            const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+            const auto& face_to_cell_matrix               = mesh.connectivity().faceToCellMatrix();
+            const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
+            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+            const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();
+
+            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) {
+              const FaceId face_id       = face_list[i_face];
+              const auto& face_cell_list = face_to_cell_matrix[face_id];
+              Assert(face_cell_list.size() == 1);
+
+              CellId face_cell_id              = face_cell_list[0];
+              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+              const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+              const auto& face_nodes = face_to_node_matrix[face_id];
+
+              for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+                NodeId node_id = face_nodes[i_node];
+                br[node_id] += sign * value_list[i_face] * Nlr(face_id, i_node);
+              }
+            }
+          }
+        }
+      },
+      boundary_condition);
+  }
+}
+
 template <size_t Dimension>
 void
 HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
@@ -745,6 +838,61 @@ class HyperelasticSolverHandler::HyperelasticSolver<1>::PressureBoundaryConditio
   ~PressureBoundaryCondition() = default;
 };
 
+template <size_t Dimension>
+class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::NormalStressBoundaryCondition
+{
+ private:
+  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+  const Array<const Rdxd> m_value_list;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_mesh_face_boundary.faceList();
+  }
+
+  const Array<const Rdxd>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  NormalStressBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
+                                const Array<const Rdxd>& value_list)
+    : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
+  {}
+
+  ~NormalStressBoundaryCondition() = default;
+};
+
+template <>
+class HyperelasticSolverHandler::HyperelasticSolver<1>::NormalStressBoundaryCondition
+{
+ private:
+  const MeshNodeBoundary<1> m_mesh_node_boundary;
+  const Array<const Rdxd> m_value_list;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  const Array<const Rdxd>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  NormalStressBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const Rdxd>& value_list)
+    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+  {}
+
+  ~NormalStressBoundaryCondition() = default;
+};
+
 template <size_t Dimension>
 class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::VelocityBoundaryCondition
 {