[beaminteraction] Individual augmentation scaling parameters

src/beaminteraction/src/contact/4C_beaminteraction_contact_submodel_evaluator.cpp

@@ -360,7 +360,8 @@ bool BeamInteraction::SubmodelEvaluator::BeamContact::evaluate_force()
   for (auto& assembly_manager : assembly_managers_)
   {
     assembly_manager->evaluate_force_stiff(discret_ptr(), beam_interaction_data_state_ptr(),
-        beam_interaction_data_state_ptr()->get_force_np(), nullptr);
+        beam_interaction_data_state_ptr()->get_force_np(), nullptr,
+        beam_contact_params_ptr_->beam_to_solid_volume_meshtying_params());
   }
 
   return true;
@@ -378,7 +379,8 @@ bool BeamInteraction::SubmodelEvaluator::BeamContact::evaluate_stiff()
   for (auto& assembly_manager : assembly_managers_)
   {
     assembly_manager->evaluate_force_stiff(discret_ptr(), beam_interaction_data_state_ptr(),
-        nullptr, beam_interaction_data_state_ptr()->get_stiff());
+        nullptr, beam_interaction_data_state_ptr()->get_stiff(),
+        beam_contact_params_ptr_->beam_to_solid_volume_meshtying_params());
   }
 
   return true;
@@ -397,7 +399,8 @@ bool BeamInteraction::SubmodelEvaluator::BeamContact::evaluate_force_stiff()
   for (auto& assembly_manager : assembly_managers_)
     assembly_manager->evaluate_force_stiff(discret_ptr(), beam_interaction_data_state_ptr(),
         beam_interaction_data_state_ptr()->get_force_np(),
-        beam_interaction_data_state_ptr()->get_stiff());
+        beam_interaction_data_state_ptr()->get_stiff(),
+        beam_contact_params_ptr_->beam_to_solid_volume_meshtying_params());
 
   print_active_beam_contact_set(Core::IO::cout.os(Core::IO::debug));
 

src/beaminteraction/src/contact/4C_beaminteraction_contact_submodel_evaluator_assembly_manager.hpp

@@ -11,6 +11,7 @@
 
 #include "4C_config.hpp"
 
+#include "4C_beaminteraction_contact_beam_to_solid_volume_meshtying_params.hpp"
 #include "4C_linalg_fevector.hpp"
 #include "4C_linalg_vector.hpp"
 #include "4C_utils_exceptions.hpp"
@@ -73,7 +74,9 @@ namespace BeamInteraction
       virtual void evaluate_force_stiff(std::shared_ptr<Core::FE::Discretization> discret,
           const std::shared_ptr<const Solid::ModelEvaluator::BeamInteractionDataState>& data_state,
           std::shared_ptr<Core::LinAlg::FEVector<double>> fe_sysvec,
-          std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat)
+          std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat,
+          const std::shared_ptr<BeamInteraction::BeamToSolidVolumeMeshtyingParams>
+              beam_to_solid_volume_meshtying_params)
       {
         FOUR_C_THROW("Not implemented!");
       }

src/beaminteraction/src/contact/4C_beaminteraction_contact_submodel_evaluator_assembly_manager_direct.cpp

@@ -36,7 +36,9 @@ void BeamInteraction::SubmodelEvaluator::BeamContactAssemblyManagerDirect::evalu
     std::shared_ptr<Core::FE::Discretization> discret,
     const std::shared_ptr<const Solid::ModelEvaluator::BeamInteractionDataState>& data_state,
     std::shared_ptr<Core::LinAlg::FEVector<double>> fe_sysvec,
-    std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat)
+    std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat,
+    const std::shared_ptr<BeamInteraction::BeamToSolidVolumeMeshtyingParams>
+        beam_to_solid_volume_meshtying_params)
 {
   // resulting discrete element force vectors of the two interacting elements
   std::vector<Core::LinAlg::SerialDenseVector> eleforce(2);

src/beaminteraction/src/contact/4C_beaminteraction_contact_submodel_evaluator_assembly_manager_direct.hpp

@@ -47,7 +47,9 @@ namespace BeamInteraction
       void evaluate_force_stiff(std::shared_ptr<Core::FE::Discretization> discret,
           const std::shared_ptr<const Solid::ModelEvaluator::BeamInteractionDataState>& data_state,
           std::shared_ptr<Core::LinAlg::FEVector<double>> fe_sysvec,
-          std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat) override;
+          std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat,
+          const std::shared_ptr<BeamInteraction::BeamToSolidVolumeMeshtyingParams>
+              beam_to_solid_volume_meshtying_params) override;
 
       /**
        * \brief Return a const reference to the contact pairs in this assembly manager.

src/beaminteraction/src/contact/4C_beaminteraction_contact_submodel_evaluator_assembly_manager_indirect.cpp

@@ -13,6 +13,9 @@
 #include "4C_beaminteraction_str_model_evaluator_datastate.hpp"
 #include "4C_comm_mpi_utils.hpp"
 #include "4C_fem_discretization.hpp"
+#include "4C_linalg_vector.hpp"
+
+#include <memory>
 
 FOUR_C_NAMESPACE_OPEN
 
@@ -24,7 +27,9 @@ void BeamInteraction::SubmodelEvaluator::BeamContactAssemblyManagerInDirect::eva
     std::shared_ptr<Core::FE::Discretization> discret,
     const std::shared_ptr<const Solid::ModelEvaluator::BeamInteractionDataState>& data_state,
     std::shared_ptr<Core::LinAlg::FEVector<double>> fe_sysvec,
-    std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat)
+    std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat,
+    const std::shared_ptr<BeamInteraction::BeamToSolidVolumeMeshtyingParams>
+        beam_to_solid_volume_meshtying_params)
 {
   if (mortar_manager_->have_lagrange_dofs())
   {
@@ -33,8 +38,67 @@ void BeamInteraction::SubmodelEvaluator::BeamContactAssemblyManagerInDirect::eva
     if (mortar_manager_->get_lagrange_formulation() ==
         BeamToSolid::BeamToSolidLagrangeFormulation::augmented)
     {
+      const auto* dof_row_map = discret->dof_row_map();
+      const auto* solid_map = &mortar_manager_->get_solid_dof_row_map();
+      const auto* beam_map = &mortar_manager_->get_beam_dof_row_map();
+
+      // Create matrix and vector for the penalty regularization contributions.
+      std::shared_ptr<Core::LinAlg::SparseMatrix> sysmat_penalty = nullptr;
+      if (fe_sysmat != nullptr)
+      {
+        sysmat_penalty = std::make_shared<Core::LinAlg::SparseMatrix>(*dof_row_map, 81);
+      }
+      std::shared_ptr<Core::LinAlg::FEVector<double>> sysvec_penalty = nullptr;
+      if (fe_sysvec != nullptr)
+      {
+        sysvec_penalty = std::make_shared<Core::LinAlg::FEVector<double>>(*dof_row_map);
+      }
+
+      // Evaluate the penalty coupling contributions.
       mortar_manager_->evaluate_force_stiff_penalty_regularization(
-          data_state, fe_sysmat, fe_sysvec);
+          data_state, sysmat_penalty, sysvec_penalty);
+      if (sysmat_penalty != nullptr) sysmat_penalty->complete();
+      if (sysvec_penalty != nullptr) sysvec_penalty->complete();
+
+      // Create the scaling vector holding the diagonal entries for the augmentation scaling
+      // matrix.
+      Core::LinAlg::Vector<double> augmentation_scaling_vector(*dof_row_map);
+      augmentation_scaling_vector.put_scalar(1.0);
+
+      auto add_scaling_values_to_vector = [&](const Core::LinAlg::Map& sub_map, double val)
+      {
+        for (int i = 0; i < sub_map.num_my_elements(); ++i)
+        {
+          int gid = sub_map.gid(i);
+          int lid_dof_row_map = dof_row_map->lid(gid);
+          if (lid_dof_row_map != -1)
+          {
+            augmentation_scaling_vector.replace_local_value(lid_dof_row_map, val);
+          }
+          else
+          {
+            FOUR_C_THROW(
+                "Global ID {} from sub map is not found in dof row map. This should not happen.",
+                gid);
+          }
+        }
+      };
+
+      add_scaling_values_to_vector(*solid_map,
+          beam_to_solid_volume_meshtying_params->get_augmentation_scaling_parameter_solid());
+      add_scaling_values_to_vector(*beam_map,
+          beam_to_solid_volume_meshtying_params->get_augmentation_scaling_parameter_beam());
+
+      // Scale the penalty contributions and add them to the system matrix and vector.
+      if (fe_sysmat != nullptr)
+      {
+        sysmat_penalty->left_scale(augmentation_scaling_vector);
+        fe_sysmat->add(*sysmat_penalty, false, 1.0, 1.0);
+      }
+      if (fe_sysvec != nullptr)
+      {
+        fe_sysvec->multiply(1.0, augmentation_scaling_vector, *sysvec_penalty, 1.0);
+      }
     }
   }
   else

src/beaminteraction/src/contact/4C_beaminteraction_contact_submodel_evaluator_assembly_manager_indirect.hpp

@@ -53,7 +53,9 @@ namespace BeamInteraction
       void evaluate_force_stiff(std::shared_ptr<Core::FE::Discretization> discret,
           const std::shared_ptr<const Solid::ModelEvaluator::BeamInteractionDataState>& data_state,
           std::shared_ptr<Core::LinAlg::FEVector<double>> fe_sysvec,
-          std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat) override;
+          std::shared_ptr<Core::LinAlg::SparseMatrix> fe_sysmat,
+          const std::shared_ptr<BeamInteraction::BeamToSolidVolumeMeshtyingParams>
+              beam_to_solid_volume_meshtying_params) override;
 
       /**
        * \brief Return a const pointer to the mortar manager.

src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_input.cpp

@@ -91,6 +91,15 @@ std::vector<Core::IO::InputSpec> BeamToSolid::valid_parameters()
           {.description = "Penalty parameter for beam-to-solid volume meshtying",
               .default_value = 0.0}),
 
+      parameter<double>("AUGMENTATION_SCALING_PARAMETER_BEAM",
+          {.description = "Augmentation scaling parameter for the beam row",
+              .default_value = 1.0,
+              .validator = Validators::positive<double>()}),
+
+      parameter<double>("AUGMENTATION_SCALING_PARAMETER_SOLID",
+          {.description = "Augmentation scaling parameter for the solid row",
+              .default_value = 1.0,
+              .validator = Validators::positive<double>()}),
 
       // This option only has an effect during a restart simulation.
       // - No:  (default) The coupling is treated the same way as during a non restart

src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_mortar_manager.hpp

@@ -262,6 +262,22 @@ namespace BeamInteraction
       return lambda_dof_rowmap_;
     }
 
+    /**
+     * \brief Returns a const reference to the beam dof row map
+     */
+    [[nodiscard]] const Core::LinAlg::Map& get_beam_dof_row_map() const
+    {
+      return *beam_dof_rowmap_;
+    }
+
+    /**
+     * \brief Returns  a const reference to the solid dof row map
+     */
+    [[nodiscard]] const Core::LinAlg::Map& get_solid_dof_row_map() const
+    {
+      return *solid_dof_rowmap_;
+    }
+
     /**
      * \brief Calculate the internal mortar penalty energy.
      * @return Global mortar energy.

src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_params_base.cpp

@@ -25,6 +25,8 @@ BeamInteraction::BeamToSolidParamsBase::BeamToSolidParamsBase()
       contact_discretization_(BeamToSolid::BeamToSolidContactDiscretization::none),
       mortar_shape_function_(BeamToSolid::BeamToSolidMortarShapefunctions::none),
       penalty_parameter_(-1.0),
+      augmentation_scaling_parameter_beam_(-1.0),
+      augmentation_scaling_parameter_solid_(-1.0),
       gauss_rule_(Core::FE::GaussRule1D::undefined),
       rotational_coupling_(false),
       lagrange_formulation_(BeamToSolid::BeamToSolidLagrangeFormulation::none)
@@ -69,6 +71,16 @@ void BeamInteraction::BeamToSolidParamsBase::set_base_params(
     if (penalty_parameter_ < 0.0)
       FOUR_C_THROW("beam-to-volume-meshtying penalty parameter must not be negative!");
 
+    if (beam_to_solid_params_list.isParameter("AUGMENTATION_SCALING_PARAMETER_BEAM") and
+        beam_to_solid_params_list.isParameter("AUGMENTATION_SCALING_PARAMETER_SOLID"))
+    {
+      // Augmentation scaling parameters.
+      augmentation_scaling_parameter_beam_ =
+          beam_to_solid_params_list.get<double>("AUGMENTATION_SCALING_PARAMETER_BEAM");
+      augmentation_scaling_parameter_solid_ =
+          beam_to_solid_params_list.get<double>("AUGMENTATION_SCALING_PARAMETER_SOLID");
+    }
+
     // Gauss rule for integration along the beam (segments).
     gauss_rule_ =
         GeometryPair::int_to_gauss_rule1_d(beam_to_solid_params_list.get<int>("GAUSS_POINTS"));

src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_params_base.hpp

@@ -117,6 +117,22 @@ namespace BeamInteraction
      */
     inline double get_penalty_parameter() const { return penalty_parameter_; }
 
+    /**
+     * \brief Returns the augmentation scaling parameter for the beam DOFs.
+     */
+    [[nodiscard]] inline double get_augmentation_scaling_parameter_beam() const
+    {
+      return augmentation_scaling_parameter_beam_;
+    }
+
+    /**
+     * \brief Returns the augmentation scaling parameter for the solid DOFs.
+     */
+    [[nodiscard]] inline double get_augmentation_scaling_parameter_solid() const
+    {
+      return augmentation_scaling_parameter_solid_;
+    }
+
     /**
      * \brief Returns the Gauss rule.
      * @return gauss rule.
@@ -157,6 +173,10 @@ namespace BeamInteraction
     //! Penalty parameter.
     double penalty_parameter_;
 
+    //! Augmentation scaling parameters
+    double augmentation_scaling_parameter_beam_;
+    double augmentation_scaling_parameter_solid_;
+
     //! Gauss rule to be used.
     Core::FE::GaussRule1D gauss_rule_;
 

src/core/linalg/src/sparse/4C_linalg_fevector.hpp

@@ -201,7 +201,10 @@ namespace Core::LinAlg
 
     //! Multiply a Core::LinAlg::MultiVector<double> with another, element-by-element.
     void multiply(double ScalarAB, const Epetra_MultiVector& A, const Epetra_MultiVector& B,
-        double ScalarThis);
+        double ScalarThis)
+    {
+      ASSERT_EPETRA_CALL(vector_->Multiply(ScalarAB, A, B, ScalarThis));
+    }
 
     //! Imports an Epetra_DistObject using the Core::LinAlg::Import object.
     void import(const Epetra_SrcDistObject& A, const Core::LinAlg::Import& Importer,