Fix normal pressure Neumann conditions on 2D solid lines

src/solid_ele/4C_solid_ele.hpp

@@ -188,6 +188,12 @@ namespace Discret::Elements
              solid_ele_property_.element_technology == ElementTechnology::shell_eas_ans;
     }
 
+    [[nodiscard]] double reference_thickness() const
+      requires(dim == 2)
+    {
+      return solid_ele_property_.reference_thickness;
+    }
+
     void vis_names(std::map<std::string, int>& names) const override;
 
     void set_integration_rule(const Core::FE::GaussIntegration& integration_rule);

src/solid_ele/4C_solid_ele_line.cpp

@@ -13,6 +13,7 @@
 #include "4C_legacy_enum_definitions_element_actions.hpp"
 #include "4C_linalg_serialdensematrix.hpp"
 #include "4C_linalg_tensor.hpp"
+#include "4C_solid_ele.hpp"
 #include "4C_solid_ele_neumann_evaluator.hpp"
 #include "4C_utils_exceptions.hpp"
 
@@ -235,6 +236,24 @@ int Discret::Elements::SolidLine<dim>::evaluate_neumann(Teuchos::ParameterList&
           return params.get("total time", -1.0);
       });
 
+  const auto& load_type = condition.parameters().get<std::string>("TYPE");
+  if (load_type == "pseudo_orthopressure" || load_type == "orthopressure")
+  {
+    if constexpr (dim == 2)
+    {
+      const auto* parent_solid = dynamic_cast<const Discret::Elements::Solid<2>*>(parent_element());
+      FOUR_C_ASSERT_ALWAYS(
+          parent_solid != nullptr, "Solid line has no two-dimensional solid parent.");
+      Discret::Elements::evaluate_normal_pressure_by_element(*this, discretization, condition, lm,
+          elevec1, elemat1, total_time, parent_solid->reference_thickness());
+      return 0;
+    }
+    else
+    {
+      FOUR_C_THROW("Normal pressure on a solid line is only supported in two dimensions.");
+    }
+  }
+
   Discret::Elements::evaluate_neumann_by_element<dim>(
       *this, discretization, condition, elevec1, total_time);
   return 0;

src/solid_ele/4C_solid_ele_neumann_evaluator.cpp

@@ -11,6 +11,7 @@
 #include "4C_fem_general_cell_type_traits.hpp"
 #include "4C_fem_general_element.hpp"
 #include "4C_fem_general_element_integration.hpp"
+#include "4C_fem_general_extract_values.hpp"
 #include "4C_fem_general_utils_gausspoints.hpp"
 #include "4C_fem_general_utils_local_connectivity_matrices.hpp"
 #include "4C_global_data.hpp"
@@ -21,6 +22,99 @@
 
 FOUR_C_NAMESPACE_OPEN
 
+namespace
+{
+  template <Core::FE::CellType celltype>
+  void evaluate_normal_pressure(Core::Elements::Element& element,
+      const Core::FE::Discretization& discretization, const Core::Conditions::Condition& condition,
+      const std::vector<int>& dof_index_array, Core::LinAlg::SerialDenseVector& force,
+      Core::LinAlg::SerialDenseMatrix* load_linearization, const double total_time,
+      const double reference_thickness, const std::string& displacement_state)
+  {
+    constexpr int num_nodes = Core::FE::num_nodes(celltype);
+    const int num_dof_per_node = element.num_dof_per_node(*element.nodes()[0]);
+    const auto& onoff = condition.parameters().get<std::vector<int>>("ONOFF");
+    const auto& values = condition.parameters().get<std::vector<double>>("VAL");
+    const auto& function_ids = condition.parameters().get<std::vector<std::optional<int>>>("FUNCT");
+    FOUR_C_ASSERT_ALWAYS(onoff.size() >= 2 && values.size() >= 2 && function_ids.size() >= 2,
+        "A two-dimensional pressure condition requires at least two entries in ONOFF, VAL, and "
+        "FUNCT.");
+    FOUR_C_ASSERT_ALWAYS(
+        onoff[0] == 1 && onoff[1] == 0, "Normal pressure must be activated on the first dof only.");
+
+    const auto displacement = discretization.get_state(displacement_state);
+    FOUR_C_ASSERT_ALWAYS(
+        displacement != nullptr, "Cannot get state vector '{}'.", displacement_state);
+    const std::vector<double> local_displacement =
+        Core::FE::extract_values(*displacement, dof_index_array);
+
+    Core::LinAlg::Matrix<num_nodes, 2> current_coordinates;
+    for (int node = 0; node < num_nodes; ++node)
+      for (int component = 0; component < 2; ++component)
+        current_coordinates(node, component) =
+            element.nodes()[node]->x()[component] +
+            local_displacement[node * num_dof_per_node + component];
+
+    const auto integration = Core::FE::create_gauss_integration<celltype>(
+        Discret::Elements::get_gauss_rule_stiffness_matrix<celltype>());
+    for (int gp = 0; gp < integration.num_points(); ++gp)
+    {
+      const auto xi = Core::Elements::evaluate_parameter_coordinate<celltype>(integration, gp);
+      Core::Elements::ElementNodes<celltype, 2> nodes{.coordinates = current_coordinates};
+      const auto shape = Core::Elements::evaluate_shape_functions_and_derivs<celltype>(xi, nodes);
+
+      Core::LinAlg::Matrix<2, 1> current_coordinate;
+      current_coordinate.multiply_tn(current_coordinates, shape.values);
+      const double function_factor =
+          function_ids[0].has_value() && function_ids[0].value() > 0
+              ? Global::Problem::instance()
+                    ->function_by_id<Core::Utils::FunctionOfSpaceTime>(function_ids[0].value())
+                    .evaluate(current_coordinate.as_span(), total_time, 0)
+              : 1.0;
+
+      Discret::Elements::add_normal_pressure_load<celltype>(shape, current_coordinates,
+          values[0] * function_factor * integration.weight(gp), reference_thickness,
+          num_dof_per_node, force, load_linearization);
+    }
+  }
+}  // namespace
+
+void Discret::Elements::evaluate_normal_pressure_by_element(Core::Elements::Element& element,
+    const Core::FE::Discretization& discretization, const Core::Conditions::Condition& condition,
+    const std::vector<int>& dof_index_array, Core::LinAlg::SerialDenseVector& force,
+    Core::LinAlg::SerialDenseMatrix* load_linearization, const double total_time,
+    const double reference_thickness)
+{
+  const auto& type = condition.parameters().get<std::string>("TYPE");
+  std::string displacement_state;
+  if (type == "pseudo_orthopressure")
+  {
+    displacement_state = "displacement";
+    load_linearization = nullptr;
+  }
+  else if (type == "orthopressure")
+  {
+    FOUR_C_ASSERT_ALWAYS(
+        Global::Problem::instance()->structural_dynamic_params().get<bool>("LOADLIN"),
+        "If you use NEUMANN CONDITIONS with TYPE: \"orthopressure\" you need to set "
+        "'LOADLIN: true' in 'STRUCTURAL DYNAMIC'.");
+    displacement_state = "displacement new";
+  }
+  else
+  {
+    FOUR_C_THROW("Expected a normal-pressure condition, got TYPE '{}'.", type);
+  }
+
+  using supported_celltypes =
+      Core::FE::CelltypeSequence<Core::FE::CellType::line2, Core::FE::CellType::line3>;
+  Core::FE::cell_type_switch<supported_celltypes>(element.shape(),
+      [&](auto celltype_t)
+      {
+        evaluate_normal_pressure<celltype_t()>(element, discretization, condition, dof_index_array,
+            force, load_linearization, total_time, reference_thickness, displacement_state);
+      });
+}
+
 template <int dim>
 void Discret::Elements::evaluate_neumann_by_element(Core::Elements::Element& element,
     const Core::FE::Discretization& discretization, const Core::Conditions::Condition& condition,

src/solid_ele/4C_solid_ele_neumann_evaluator.hpp

@@ -10,8 +10,12 @@
 #include "4C_config.hpp"
 
 #include "4C_fem_general_element.hpp"
+#include "4C_fem_general_element_integration.hpp"
+#include "4C_linalg_fixedsizematrix.hpp"
+#include "4C_linalg_serialdensematrix.hpp"
 #include "4C_linalg_serialdensevector.hpp"
 
+#include <array>
 #include <vector>
 
 FOUR_C_NAMESPACE_OPEN
@@ -23,6 +27,60 @@ namespace Core::FE
 
 namespace Discret::Elements
 {
+  template <Core::FE::CellType celltype>
+    requires(Core::FE::dim<celltype> == 1)
+  /*!
+   * @brief Add pressure normal to a plane solid boundary.
+   *
+   * The pressure is multiplied by the parent solid's reference thickness. If requested, this
+   * function also assembles the consistent negative load linearization used by the structural
+   * residual.
+   */
+  void add_normal_pressure_load(const Core::Elements::ShapeFunctionsAndDerivatives<celltype>& shape,
+      const Core::LinAlg::Matrix<Core::FE::num_nodes(celltype), 2>& current_coordinates,
+      const double pressure_times_weight, const double reference_thickness,
+      const int num_dof_per_node, Core::LinAlg::SerialDenseVector& force,
+      Core::LinAlg::SerialDenseMatrix* load_linearization)
+  {
+    double radial_derivative = 0.0;
+    double axial_derivative = 0.0;
+    for (int node = 0; node < Core::FE::num_nodes(celltype); ++node)
+    {
+      radial_derivative += shape.derivatives(0, node) * current_coordinates(node, 0);
+      axial_derivative += shape.derivatives(0, node) * current_coordinates(node, 1);
+    }
+
+    const std::array<double, 2> normal_measure{
+        reference_thickness * axial_derivative, -reference_thickness * radial_derivative};
+    for (int node = 0; node < Core::FE::num_nodes(celltype); ++node)
+      for (int component = 0; component < 2; ++component)
+        force[node * num_dof_per_node + component] +=
+            shape.values(node) * pressure_times_weight * normal_measure[component];
+
+    if (load_linearization == nullptr) return;
+
+    for (int force_node = 0; force_node < Core::FE::num_nodes(celltype); ++force_node)
+      for (int coordinate_node = 0; coordinate_node < Core::FE::num_nodes(celltype);
+          ++coordinate_node)
+      {
+        const double factor = shape.values(force_node) * pressure_times_weight *
+                              reference_thickness * shape.derivatives(0, coordinate_node);
+        (*load_linearization)(
+            force_node* num_dof_per_node, coordinate_node* num_dof_per_node + 1) -= factor;
+        (*load_linearization)(
+            force_node* num_dof_per_node + 1, coordinate_node * num_dof_per_node) += factor;
+      }
+  }
+
+  /*!
+   * @brief Evaluate pseudo-orthopressure or follower orthopressure on a plane solid boundary line.
+   */
+  void evaluate_normal_pressure_by_element(Core::Elements::Element& element,
+      const Core::FE::Discretization& discretization, const Core::Conditions::Condition& condition,
+      const std::vector<int>& dof_index_array, Core::LinAlg::SerialDenseVector& force,
+      Core::LinAlg::SerialDenseMatrix* load_linearization, double total_time,
+      double reference_thickness);
+
   /*!
    * @brief Evaluates a Neumann condition @p condition for the element @p element.
    *