diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 0f5de1a9bf388b250e4ce8963fa856e7860b7f22..d7c8755b05e082df8478d755200123e92fd16101 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/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,
diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index 9a94c787f8be7cf4bf64060418a3c0e8ae9eeb0f..81717e62bdb9aab35849a9efb4a2f8fd835766c2 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/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
{