Merge pull request #1968 from amgebauer/move-constr2D-tests-to-new-solids

Move all constr2D tests to new solid elements

Author: amgebauer

src/adapter/4C_adapter_str_structure_new.cpp

@@ -606,9 +606,13 @@ void Adapter::StructureBaseAlgorithmNew::detect_element_technologies(
       if (shell7p->get_ele_tech().find(Inpar::Solid::EleTech::eas) != shell7p->get_ele_tech().end())
         iseas_local = 1;
 
-    Discret::Elements::Solid<3>* solid = dynamic_cast<Discret::Elements::Solid<3>*>(actele);
-    if (solid != nullptr)
-      if (solid->have_eas()) iseas_local = 1;
+    const auto* solid3 = dynamic_cast<Discret::Elements::Solid<3>*>(actele);
+    if (solid3 != nullptr)
+      if (solid3->have_eas()) iseas_local = 1;
+
+    const auto* solid2 = dynamic_cast<Discret::Elements::Solid<2>*>(actele);
+    if (solid2 != nullptr)
+      if (solid2->have_eas()) iseas_local = 1;
 
     // Detect non-additive rotation-vector DOFs --------------------------------
     if (actele->element_type() == Discret::Elements::Beam3rType::instance() or

src/core/legacy_enum_definitions/4C_legacy_enum_definitions_element_actions.hpp

@@ -41,6 +41,8 @@ namespace Core::Elements
     struct_calc_reset_istep,  //!< reset elementwise internal variables, during iteration to last
                               //!< converged state
     struct_calc_energy,       //!< compute internal energy
+    calc_struct_constrarea,
+    calc_struct_areaconstrstiff,
     struct_postprocess_thickness,         //!< postprocess thickness of membrane finite elements
     struct_init_gauss_point_data_output,  //!< initialize quantities for output of gauss point
                                           //!< data
@@ -99,6 +101,10 @@ namespace Core::Elements
       return struct_calc_reset_istep;
     else if (action == "calc_struct_energy")
       return struct_calc_energy;
+    else if (action == "calc_struct_constrarea")
+      return calc_struct_constrarea;
+    else if (action == "calc_struct_areaconstrstiff")
+      return calc_struct_areaconstrstiff;
     else if (action == "struct_init_gauss_point_data_output")
       return struct_init_gauss_point_data_output;
     else if (action == "struct_gauss_point_data_output")
@@ -150,6 +156,10 @@ namespace Core::Elements
         return "struct_calc_stress";
       case struct_calc_thickness:
         return "struct_calc_thickness";
+      case calc_struct_constrarea:
+        return "calc_struct_constrarea";
+      case calc_struct_areaconstrstiff:
+        return "calc_struct_areaconstrstiff";
       case struct_calc_eleload:
         return "struct_calc_eleload";
       case struct_calc_fsiload:

src/solid_3D_ele/4C_solid_3D_ele_calc_lib_plane.cpp

@@ -313,8 +313,10 @@ double Discret::Elements::evaluate_material_strain_energy(Mat::So3Material& mate
     const Core::LinAlg::SymmetricTensor<double, 2, 2>& gl_strain, Teuchos::ParameterList& params,
     const Mat::EvaluationContext<2>& context, const int gp, const int eleGID)
 {
+  constexpr auto dummy_defgrd =
+      Core::LinAlg::get_full(Core::LinAlg::TensorGenerators::identity<double, 2, 2>);
   const Core::LinAlg::Tensor<double, 2, 2> defgrd =
-      Discret::Elements::compute_deformation_gradient_from_gl_strains({}, gl_strain);
+      Discret::Elements::compute_deformation_gradient_from_gl_strains(dummy_defgrd, gl_strain);
 
   double strain_energy = 0.0;
   transform_to_3d(material, element_properties, defgrd, gl_strain, params, context, gp, eleGID,

src/solid_3D_ele/4C_solid_3D_ele_line.cpp

@@ -7,7 +7,12 @@
 
 #include "4C_solid_3D_ele_line.hpp"
 
+#include "4C_fem_discretization.hpp"
+#include "4C_fem_general_cell_type_traits.hpp"
+#include "4C_fem_general_extract_values.hpp"
+#include "4C_legacy_enum_definitions_element_actions.hpp"
 #include "4C_linalg_serialdensematrix.hpp"
+#include "4C_linalg_tensor.hpp"
 #include "4C_solid_3D_ele_neumann_evaluator.hpp"
 #include "4C_utils_exceptions.hpp"
 
@@ -19,6 +24,32 @@
 FOUR_C_NAMESPACE_OPEN
 
 
+namespace
+{
+  template <Core::FE::CellType celltype, unsigned dim>
+  std::array<Core::LinAlg::Tensor<double, dim>, Core::FE::num_nodes(celltype)>
+  evaluate_current_coordinates(const Core::Elements::Element& ele,
+      const Core::FE::Discretization& discretization, const std::vector<int>& lm)
+  {
+    constexpr unsigned num_dofs = dim * Core::FE::num_nodes(celltype);
+    std::array<double, num_dofs> mydisp =
+        Core::FE::extract_values_as_array<num_dofs>(*discretization.get_state("displacement"), lm);
+
+    std::array<Core::LinAlg::Tensor<double, dim>, Core::FE::num_nodes(celltype)>
+        current_nodal_coordinates;
+
+    for (int i = 0; i < Core::FE::num_nodes(celltype); ++i)
+    {
+      for (std::size_t d = 0; d < dim; ++d)
+      {
+        current_nodal_coordinates[i](d) = ele.nodes()[i]->x()[d] + mydisp[i * dim + d];
+      }
+    }
+
+    return current_nodal_coordinates;
+  }
+}  // namespace
+
 template <unsigned dim>
 Discret::Elements::SolidLineType<dim> Discret::Elements::SolidLineType<dim>::instance_;
 
@@ -107,12 +138,92 @@ void Discret::Elements::SolidLine<dim>::print(std::ostream& os) const
   Element::print(os);
 }
 
+template <unsigned dim>
+int Discret::Elements::SolidLine<dim>::evaluate(Teuchos::ParameterList& params,
+    Core::FE::Discretization& discretization, std::vector<int>& lm,
+    Core::LinAlg::SerialDenseMatrix& elematrix1, Core::LinAlg::SerialDenseMatrix& elematrix2,
+    Core::LinAlg::SerialDenseVector& elevector1, Core::LinAlg::SerialDenseVector& elevector2,
+    Core::LinAlg::SerialDenseVector& elevector3)
+{
+  set_params_interface_ptr(params);
+
+  Core::Elements::ActionType act =
+      Core::Elements::string_to_action_type(params.get<std::string>("action"));
+
+  switch (act)
+  {
+    case Core::Elements::ActionType::calc_struct_constrarea:
+    {
+      FOUR_C_ASSERT_ALWAYS(
+          shape() == Core::FE::CellType::line2, "Area Constraint only works for line2 elements!");
+
+      constexpr Core::FE::CellType celltype = Core::FE::CellType::line2;
+
+
+      // We are not interested in volume of ghosted elements
+      if (Core::Communication::my_mpi_rank(discretization.get_comm()) != owner()) return 0;
+
+      auto current_nodal_coordinates =
+          evaluate_current_coordinates<celltype, dim>(*this, discretization, lm);
+
+      elevector3[0] = 0.5 * (current_nodal_coordinates[0](1) + current_nodal_coordinates[1](1)) *
+                      (current_nodal_coordinates[1](0) - current_nodal_coordinates[0](0));
+    }
+    break;
+    case Core::Elements::ActionType::calc_struct_areaconstrstiff:
+    {
+      FOUR_C_ASSERT_ALWAYS(
+          shape() == Core::FE::CellType::line2, "Area Constraint only works for line2 elements!");
+
+      constexpr Core::FE::CellType celltype = Core::FE::CellType::line2;
+
+      // We are not interested in volume of ghosted elements
+      if (Core::Communication::my_mpi_rank(discretization.get_comm()) != owner()) return 0;
+
+      auto current_nodal_coordinates =
+          evaluate_current_coordinates<celltype, dim>(*this, discretization, lm);
+
+      FOUR_C_ASSERT(elevector1.length() == 4,
+          "Expected elevector1 to have size 4 for line2 area constraint! Given vector has size {}",
+          elevector1.length());
+      FOUR_C_ASSERT(elematrix1.num_cols() == 4 && elematrix1.num_rows() == 4,
+          "Expected elematrix1 to have 4 rows and 4 columns for line2 area constraint! Given "
+          "matrix has {} rows and {} columns",
+          elematrix1.num_rows(), elematrix1.num_cols());
+
+      elevector1[0] = 0.5 * (current_nodal_coordinates[0](1) + current_nodal_coordinates[1](1));
+      elevector1[1] = 0.5 * (current_nodal_coordinates[0](0) - current_nodal_coordinates[1](0));
+      elevector1[2] = -0.5 * (current_nodal_coordinates[0](1) + current_nodal_coordinates[1](1));
+      elevector1[3] = 0.5 * (current_nodal_coordinates[0](0) - current_nodal_coordinates[1](0));
+      elevector2 = elevector1;
+
+      elematrix1.put_scalar(0.0);
+      elematrix1(0, 1) = 0.5;
+      elematrix1(0, 3) = 0.5;
+      elematrix1(1, 0) = 0.5;
+      elematrix1(1, 2) = -0.5;
+      elematrix1(2, 1) = -0.5;
+      elematrix1(2, 3) = -0.5;
+      elematrix1(3, 0) = 0.5;
+      elematrix1(3, 2) = -0.5;
+
+      elevector3[0] = 0.5 * (current_nodal_coordinates[0](1) + current_nodal_coordinates[1](1)) *
+                      (current_nodal_coordinates[1](0) - current_nodal_coordinates[0](0));
+    }
+    break;
+    default:
+      FOUR_C_THROW("Unknown type of action {} for SolidLine element", action_type_to_string(act));
+  }
+  return 0;
+}
+
 template <unsigned dim>
 int Discret::Elements::SolidLine<dim>::evaluate_neumann(Teuchos::ParameterList& params,
     Core::FE::Discretization& discretization, const Core::Conditions::Condition& condition,
     std::vector<int>& lm, Core::LinAlg::SerialDenseVector& elevec1,
     Core::LinAlg::SerialDenseMatrix* elemat1)
 {
+  static_assert(dim == 2 || dim == 3, "SolidLine element only implemented for 2D and 3D!");
   set_params_interface_ptr(params);
 
   const double total_time = std::invoke(

src/solid_3D_ele/4C_solid_3D_ele_line.hpp

@@ -97,6 +97,12 @@ namespace Discret::Elements
       return SolidLineType<dim>::instance();
     }
 
+    int evaluate(Teuchos::ParameterList& params, Core::FE::Discretization& discretization,
+        std::vector<int>& lm, Core::LinAlg::SerialDenseMatrix& elematrix1,
+        Core::LinAlg::SerialDenseMatrix& elematrix2, Core::LinAlg::SerialDenseVector& elevector1,
+        Core::LinAlg::SerialDenseVector& elevector2,
+        Core::LinAlg::SerialDenseVector& elevector3) override;
+
     int evaluate_neumann(Teuchos::ParameterList& params, Core::FE::Discretization& discretization,
         const Core::Conditions::Condition& condition, std::vector<int>& lm,
         Core::LinAlg::SerialDenseVector& elevec1,