Merge pull request #1691 from isteinbrecher/beam-to-beam-point-couplings_with_lagrange_multipliers

Beam to beam point couplings with Lagrange multipliers

Author: isteinbrecher

src/beaminteraction/src/4C_beaminteraction_beam_to_beam_point_coupling_pair.hpp

@@ -55,8 +55,10 @@ namespace BeamInteraction
     //! Flag on type of constraint enforcement.
     enum class ConstraintEnforcement
     {
-      penalty_direct,
-      penalty_indirect
+      penalty_direct,      //! Direct assembly via penalty method.
+      penalty_indirect,    //! Set up full Lagrange-multiplier system and regularize
+      lagrange_multiplier  //! Set up full Lagrange-multiplier system and solve with Lagrange
+                           //! multipliers
     };
     ConstraintEnforcement constraint_enforcement = ConstraintEnforcement::penalty_direct;
 
@@ -238,8 +240,11 @@ namespace BeamInteraction
      */
     double get_energy() const override
     {
-      FOUR_C_THROW("get_energy not implemented yet!");
-      return 0.0;
+      if (parameters_.constraint_enforcement ==
+          BeamToBeamPointCouplingPairParameters::ConstraintEnforcement::lagrange_multiplier)
+        return 0.0;
+      else
+        FOUR_C_THROW("get_energy not implemented yet!");
     }
 
     /**

src/beaminteraction/src/4C_beaminteraction_beam_to_beam_point_coupling_pair_condition.cpp

@@ -187,6 +187,15 @@ BeamInteraction::BeamToBeamPointCouplingConditionIndirect::create_indirect_assem
         .penalty_parameter_translational = parameters_.penalty_parameter_pos,
         .penalty_parameter_rotational = parameters_.penalty_parameter_rot};
 
+    if (parameters_.constraint_enforcement ==
+        BeamToBeamPointCouplingPairParameters::ConstraintEnforcement::lagrange_multiplier)
+    {
+      mortar_manager_parameters.constraint_enforcement =
+          Inpar::BeamToSolid::BeamToSolidConstraintEnforcement::lagrange;
+      mortar_manager_parameters.lagrange_formulation =
+          Inpar::BeamToSolid::BeamToSolidLagrangeFormulation::saddlepoint;
+    }
+
     auto mortar_manager = std::make_shared<BeamInteraction::BeamToSolidMortarManager>(
         discret, mortar_manager_parameters);
     mortar_manager->setup();

src/beaminteraction/src/4C_beaminteraction_beam_to_solid_mortar_manager.cpp

@@ -431,14 +431,23 @@ void BeamInteraction::BeamToSolidMortarManager::evaluate_force_stiff_penalty_reg
  *
  */
 void BeamInteraction::BeamToSolidMortarManager::evaluate_coupling_terms_lagrange(
-    const std::shared_ptr<const Solid::ModelEvaluator::BeamInteractionDataState>& data_state)
+    const std::shared_ptr<const Solid::ModelEvaluator::BeamInteractionDataState>& data_state,
+    std::shared_ptr<Core::LinAlg::SparseMatrix> stiff,
+    std::shared_ptr<Core::LinAlg::FEVector<double>> force)
 {
   // Evaluate the global coupling terms
   evaluate_and_assemble_global_coupling_contributions(data_state->get_dis_col_np());
 
   // For now we need to set a link to the global Lagrange multiplier vector in the beam interaction
   // data state.
   global_lambda_ = data_state->get_lambda();
+
+  // Add the force and stiffness contributions that are assembled directly by the pairs.
+  Core::LinAlg::Vector<double> lambda_col(*lambda_dof_colmap_);
+  Core::LinAlg::export_to(*global_lambda_, lambda_col);
+  for (const auto& elepairptr : contact_pairs_)
+    elepairptr->evaluate_and_assemble(
+        *discret_, this, force, stiff, lambda_col, *data_state->get_dis_col_np());
 }
 
 /**
@@ -780,6 +789,23 @@ void BeamInteraction::BeamToSolidMortarManager::assemble_stiff(
   auto block_lm_displ_row_map = jac_block_sparse_matrix_base->matrix(1, 0).row_map();
   auto block_displ_lm_row_map = jac_block_sparse_matrix_base->matrix(0, 1).row_map();
 
+  Core::LinAlg::SparseMatrix lm_lm =
+      Core::LinAlg::SparseMatrix(block_lm_displ_row_map, 81, true, true);
+
+  auto lambda_non_active = Core::LinAlg::Vector<double>(lambda_active_->get_map());
+  for (int lid = 0; lid < lambda_active_->get_map().num_my_elements(); lid++)
+  {
+    if (lambda_active_->get_values()[lid] < 0.1)
+      lambda_non_active.replace_local_value(lid, 1.0);
+    else
+      lambda_non_active.replace_local_value(lid, 0.0);
+  }
+  auto lambda_non_active_structure_map = Core::LinAlg::Vector<double>(block_lm_displ_row_map);
+  Core::LinAlg::export_to(lambda_non_active, lambda_non_active_structure_map);
+  Core::LinAlg::SparseMatrix lambda_non_active_matrix(lambda_non_active_structure_map);
+  lambda_non_active_matrix.complete();
+  lm_lm.add(lambda_non_active_matrix, false, 1.0, 1.0);
+
   if (parameters_.lagrange_formulation ==
       Inpar::BeamToSolid::BeamToSolidLagrangeFormulation::regularized)
   {
@@ -790,10 +816,13 @@ void BeamInteraction::BeamToSolidMortarManager::assemble_stiff(
     Core::LinAlg::SparseMatrix kappa_penalty_inv_mat(kappa_vector);
     kappa_penalty_inv_mat.scale(-1.0 / penalty_translation);
     kappa_penalty_inv_mat.complete();
-    gstate.assign_model_block(jac, kappa_penalty_inv_mat, Inpar::Solid::model_beaminteraction,
-        Solid::MatBlockType::lm_lm);
+    lm_lm.add(kappa_penalty_inv_mat, false, 1.0, 1.0);
   }
 
+  gstate.assign_model_block(
+      jac, lm_lm, Inpar::Solid::model_beaminteraction, Solid::MatBlockType::lm_lm);
+
+
   Core::LinAlg::SparseMatrix lm_displ =
       Core::LinAlg::SparseMatrix(*lambda_dof_rowmap_, 81, true, true);
   Core::LinAlg::matrix_add(*constraint_lin_beam_, false, 1.0, lm_displ, 0.0);

src/beaminteraction/src/4C_beaminteraction_beam_to_solid_mortar_manager.hpp

@@ -214,9 +214,13 @@ namespace BeamInteraction
      * \brief Evaluate the coupling contributions resulting from a lagrange
      * type of constraint enforcement.
      * @param data_state (in) Beam interaction data state.
+     * @param stiff (out) Global stiffness matrix. Mortar terms are added to it.
+     * @param force (out) Global force vector. Mortar terms are added to it.
      */
     void evaluate_coupling_terms_lagrange(
-        const std::shared_ptr<const Solid::ModelEvaluator::BeamInteractionDataState>& data_state);
+        const std::shared_ptr<const Solid::ModelEvaluator::BeamInteractionDataState>& data_state,
+        std::shared_ptr<Core::LinAlg::SparseMatrix> stiff,
+        std::shared_ptr<Core::LinAlg::FEVector<double>> force);
 
     /**
      * \brief Check if the model uses Lagrange multipliers.

src/beaminteraction/src/4C_beaminteraction_submodel_evaluator_beamcontact_assembly_manager_indirect.cpp

@@ -29,7 +29,7 @@ void BeamInteraction::SubmodelEvaluator::BeamContactAssemblyManagerInDirect::eva
 {
   if (mortar_manager_->have_lagrange_dofs())
   {
-    mortar_manager_->evaluate_coupling_terms_lagrange(data_state);
+    mortar_manager_->evaluate_coupling_terms_lagrange(data_state, fe_sysmat, fe_sysvec);
   }
   else
   {

src/inpar/4C_inpar_beaminteraction.cpp

@@ -229,8 +229,8 @@ void Inpar::BeamInteraction::set_valid_conditions(
         Core::Conditions::geometry_type_line);
     penalty_coupling_condition_indirect.add_component(parameter<int>("COUPLING_ID"));
     penalty_coupling_condition_indirect.add_component(group("PARAMETERS",
-        {parameter<double>("POSITIONAL_PENALTY_PARAMETER"),
-            parameter<double>("ROTATIONAL_PENALTY_PARAMETER"),
+        {parameter<double>("POSITIONAL_PENALTY_PARAMETER", {.default_value = 0.0}),
+            parameter<double>("ROTATIONAL_PENALTY_PARAMETER", {.default_value = 0.0}),
             parameter<double>("PROJECTION_VALID_FACTOR",
                 {.description = "Factor multiplied with sum of cross section "
                                 "radii to define valid projection distance",
@@ -240,10 +240,8 @@ void Inpar::BeamInteraction::set_valid_conditions(
                     .default_value = 0}),
             parameter<::FourC::BeamInteraction::BeamToBeamPointCouplingPairParameters::
                     ConstraintEnforcement>("CONSTRAINT_ENFORCEMENT",
-                {.description = "How the constraints for this condition shall be enforced",
-                    .default_value =
-                        ::FourC::BeamInteraction::BeamToBeamPointCouplingPairParameters::
-                            ConstraintEnforcement::penalty_direct})},
+                {.description = "Specifies the constraint enforcement technique for beam-to-beam "
+                                "point couplings."})},
         {.required = false}));
     condlist.push_back(penalty_coupling_condition_indirect);
   }

…eam_penalty_point_coupling_elbow.4C.yaml …eam_to_beam_point_coupling_elbow.4C.yamltests/input_files/beam3r_herm2line3_beam_to_beam_penalty_point_coupling_elbow.4C.yaml renamed to tests/input_files/beam3r_herm2line3_beam_to_beam_point_coupling_elbow.4C.yaml tests/input_files/beam3r_herm2line3_beam_to_beam_penalty_point_coupling_elbow.4C.yaml renamed to tests/input_files/beam3r_herm2line3_beam_to_beam_point_coupling_elbow.4C.yaml

@@ -39,7 +39,7 @@ BINNING STRATEGY:
   BIN_SIZE_LOWER_BOUND: 10
   DOMAINBOUNDINGBOX: "-20 -20 -20 20 20 20"
 STRUCT NOX/Status Test:
-  XML File: "beam3r_herm2line3_beam_to_beam_penalty_point_coupling_elbow.xml"
+  XML File: "beam3r_herm2line3_beam_to_beam_point_coupling_elbow.xml"
 MATERIALS:
   - MAT: 1
     MAT_BeamReissnerElastHyper:

…to_beam_penalty_point_coupling_elbow.xml …e3_beam_to_beam_point_coupling_elbow.xmltests/input_files/beam3r_herm2line3_beam_to_beam_penalty_point_coupling_elbow.xml renamed to tests/input_files/beam3r_herm2line3_beam_to_beam_point_coupling_elbow.xml tests/input_files/beam3r_herm2line3_beam_to_beam_penalty_point_coupling_elbow.xml renamed to tests/input_files/beam3r_herm2line3_beam_to_beam_point_coupling_elbow.xml

@@ -61,7 +61,7 @@ STRUCT NOX/Printing:
   Linear Solver Details: true
   Test Details: true
 STRUCT NOX/Status Test:
-  XML File: "beam3r_herm2line3_beam_to_beam_penalty_point_coupling_elbow.xml"
+  XML File: "beam3r_herm2line3_beam_to_beam_point_coupling_elbow.xml"
 STRUCTURAL DYNAMIC:
   LINEAR_SOLVER: 1
   INT_STRATEGY: "Standard"

…alty_point_coupling_elbow_offset.4C.yaml …beam_point_coupling_elbow_offset.4C.yamltests/input_files/beam3r_herm2line3_beam_to_beam_penalty_point_coupling_elbow_offset.4C.yaml renamed to tests/input_files/beam3r_herm2line3_beam_to_beam_point_coupling_elbow_offset.4C.yaml tests/input_files/beam3r_herm2line3_beam_to_beam_penalty_point_coupling_elbow_offset.4C.yaml renamed to tests/input_files/beam3r_herm2line3_beam_to_beam_point_coupling_elbow_offset.4C.yaml

@@ -71,7 +71,7 @@ STRUCT NOX/Printing:
   Linear Solver Details: true
   Test Details: true
 STRUCT NOX/Status Test:
-  XML File: "beam3r_herm2line3_beam_to_beam_penalty_point_coupling_elbow.xml"
+  XML File: "beam3r_herm2line3_beam_to_beam_point_coupling_elbow.xml"
 STRUCTURAL DYNAMIC:
   LINEAR_SOLVER: 1
   INT_STRATEGY: "Standard"