Evaluate BTB penalty parameters by a function

src/beaminteraction/src/4C_beaminteraction_beam_to_beam_contact_pair.cpp

@@ -27,6 +27,7 @@
 #include "4C_linalg_utils_sparse_algebra_math.hpp"
 #include "4C_structure_timint_impl.hpp"
 #include "4C_utils_exceptions.hpp"
+#include "4C_utils_function_of_time.hpp"
 
 #include <Teuchos_TimeMonitor.hpp>
 
@@ -130,7 +131,6 @@ void BeamInteraction::BeamToBeamContactPair<numnodes, numnodalvalues>::setup()
         "contact search!",
         element1()->id(), element2()->id());
   }
-
   cpvariables_.resize(0);
   gpvariables_.resize(0);
   epvariables_.resize(0);
@@ -242,8 +242,27 @@ bool BeamInteraction::BeamToBeamContactPair<numnodes, numnodalvalues>::evaluate(
   // since most of the elements leave directly after the closest point projection!
   clear_class_variables();
 
+  // initialize the penalty parameter with default one
   double pp = params()->beam_to_beam_contact_params()->beam_to_beam_point_penalty_param();
 
+  // check if a function is specified
+  auto function_id =
+      params()->beam_to_beam_contact_params()->beam_to_beam_function_id_for_point_penalty();
+  if (function_id > -1)
+  {
+    // get penalty the value
+    auto& funct =
+        Global::Problem::instance()->function_by_id<Core::Utils::FunctionOfTime>(function_id);
+    pp = funct.evaluate(time_);
+    if (pp < 0.0)
+    {
+      FOUR_C_THROW("Penalty value is smaller than 0, Please check Function {} at time {}",
+          function_id, time_);
+    }
+  }
+
+
+
   // Subdevide the two elements in segments with linear approximation
   std::vector<Core::LinAlg::Matrix<3, 1, double>> endpoints1;
   std::vector<Core::LinAlg::Matrix<3, 1, double>> endpoints2;
@@ -691,10 +710,29 @@ void BeamInteraction::BeamToBeamContactPair<numnodes, numnodalvalues>::get_activ
                 std::pair<int, int> leftpoint_ids = std::make_pair(leftpoint_id1, leftpoint_id2);
                 TYPE jacobi = get_jacobi_at_xi(element1(), eta1_slave) * jacobi_interval;
 
-                const double parallel_pp =
+                // initialize the penalty parameter with parameter
+                double parallel_pp =
                     params()->beam_to_beam_contact_params()->beam_to_beam_line_penalty_param();
 
-                if (parallel_pp < 0.0) FOUR_C_THROW("BEAMS_BTBLINEPENALTYPARAM not set!");
+                // check if a function is specified
+                auto function_id = params()
+                                       ->beam_to_beam_contact_params()
+                                       ->beam_to_beam_function_id_for_line_penalty();
+                if (function_id > -1)
+                {
+                  // get penalty the value
+                  auto& funct =
+                      Global::Problem::instance()->function_by_id<Core::Utils::FunctionOfTime>(
+                          function_id);
+                  parallel_pp = funct.evaluate(time_);
+                }
+
+
+                if (parallel_pp < 0.0)
+                  FOUR_C_THROW(
+                      "The peanlty parameter from BEAMS_BTBLINEPENALTYPARAM or "
+                      "BEAMS_BTBLINEPENALTYPARAM_BY_FUNCT may not be smaller than 0!");
+
 
                 // Create data container for each Gauss point (in case of small-angle contact the
                 // number of the Gauss point [numgp] and the number of the integration interval

src/beaminteraction/src/4C_beaminteraction_beam_to_beam_contact_params.cpp

@@ -26,7 +26,9 @@ BeamInteraction::BeamToBeamContactParams::BeamToBeamContactParams()
       btb_penalty_law_regularization_c0_(-1.0),
       gap_shift_(0.0),
       btb_point_penalty_param_(-1.0),
+      btb_function_id_for_point_penalty_(-1),
       btb_line_penalty_param_(-1.0),
+      btb_function_id_for_line_penalty_(-1),
       btb_perp_shifting_angle1_(-1.0),
       btb_perp_shifting_angle2_(-1.0),
       btb_parallel_shifting_angle1_(-1.0),
@@ -83,20 +85,41 @@ void BeamInteraction::BeamToBeamContactParams::init()
   /****************************************************************************/
   btb_point_penalty_param_ = beam_contact_params_list.get<double>("BEAMS_BTBPENALTYPARAM");
 
-  if (btb_point_penalty_param_ < 0.0)
-    FOUR_C_THROW("beam-to-beam point penalty parameter must not be negative!");
+  btb_function_id_for_point_penalty_ =
+      beam_contact_params_list.get<int>("BEAMS_BTBPENALTYPARAM_BY_FUNCT");
+
+  if (btb_point_penalty_param_ < 0.0 and btb_function_id_for_point_penalty_ < 0.0)
+    FOUR_C_THROW(
+        "beam-to-beam point penalty parameter or beam-to-beam function id for point penalty must "
+        "be specified!");
+
+  if (btb_point_penalty_param_ > 0.0 and btb_function_id_for_point_penalty_ > 0.0)
+    FOUR_C_THROW(
+        "Please specify the beam-to-beam point penalty parameter either with BEAMS_BTBPENALTYPARAM "
+        "or BEAMS_BTBPENALTYPARAM_BY_FUNCT. Not both!");
 
 
   // input parameters required for all-angle-beam contact formulation ...
   if (beam_contact_params_list.get<bool>("BEAMS_SEGCON"))
   {
     /****************************************************************************/
+
+
+    btb_function_id_for_line_penalty_ =
+        beam_contact_params_list.get<int>("BEAMS_BTBLINEPENALTYPARAM_BY_FUNCT");
+
     btb_line_penalty_param_ = beam_contact_params_list.get<double>("BEAMS_BTBLINEPENALTYPARAM");
 
-    if (btb_line_penalty_param_ < 0.0)
+    if (btb_line_penalty_param_ < 0.0 and btb_function_id_for_line_penalty_ < 0.0)
       FOUR_C_THROW(
-          "You chose all-angle-beam contact algorithm: thus, beam-to-beam line"
-          " penalty parameter must not be negative!");
+          " beam-to-beam line penalty parameter or beam-to-beam function id for line penalty must "
+          "be specified non negative!");
+
+    if (btb_line_penalty_param_ > 0.0 and btb_function_id_for_line_penalty_ > 0.0)
+      FOUR_C_THROW(
+          "Please specify the beam-to-beam line penalty parameter either with "
+          "BEAMS_BTBLINEPENALTYPARAM or BEAMS_BTBLINEPENALTYPARAM_BY_FUNCT. Not both!");
+
 
     /****************************************************************************/
     // Todo find more verbose and expressive naming

src/beaminteraction/src/4C_beaminteraction_beam_to_beam_contact_params.hpp

@@ -74,6 +74,18 @@ namespace BeamInteraction
 
     inline double beam_to_beam_point_penalty_param() const { return btb_point_penalty_param_; }
 
+    inline double beam_to_beam_function_id_for_point_penalty() const
+    {
+      return btb_function_id_for_point_penalty_;
+    }
+
+    inline double beam_to_beam_function_id_for_line_penalty() const
+    {
+      return btb_function_id_for_line_penalty_;
+    }
+
+
+
     inline double beam_to_beam_line_penalty_param() const { return btb_line_penalty_param_; }
 
     inline double beam_to_beam_perp_shifting_angle1() const { return btb_perp_shifting_angle1_; }
@@ -119,9 +131,15 @@ namespace BeamInteraction
     //! beam-to-beam point penalty parameter
     double btb_point_penalty_param_;
 
+    //! function id for the beam-to-beam point penalty evaluated from a function
+    int btb_function_id_for_point_penalty_;
+
     //! beam-to-beam line penalty parameter
     double btb_line_penalty_param_;
 
+    //! function id for beam-to-beam line penalty parameter
+    int btb_function_id_for_line_penalty_;
+
     //! shifting angles [radians] for point contact (near perpendicular configurations) fade
     double btb_perp_shifting_angle1_;
     double btb_perp_shifting_angle2_;

src/beaminteraction/src/4C_beaminteraction_contact_pair.cpp

@@ -42,7 +42,7 @@ BeamInteraction::BeamContactPair::BeamContactPair()
  *----------------------------------------------------------------------------*/
 void BeamInteraction::BeamContactPair::init(
     const std::shared_ptr<BeamInteraction::BeamContactParams> params_ptr,
-    std::vector<Core::Elements::Element const*> elements)
+    std::vector<Core::Elements::Element const*> elements, double time)
 {
   issetup_ = false;
 
@@ -52,6 +52,8 @@ void BeamInteraction::BeamContactPair::init(
   element2_ = elements[1];
 
   isinit_ = true;
+
+  time_ = time;
 }
 
 /*----------------------------------------------------------------------------*

src/beaminteraction/src/4C_beaminteraction_contact_pair.hpp

@@ -81,7 +81,7 @@ namespace BeamInteraction
     virtual ~BeamContactPair() = default;
     //! Initialization
     virtual void init(const std::shared_ptr<BeamInteraction::BeamContactParams> params_ptr,
-        std::vector<Core::Elements::Element const*> elements);
+        std::vector<Core::Elements::Element const*> elements, double time = 0);
 
     //! Setup
     virtual void setup();
@@ -419,6 +419,9 @@ namespace BeamInteraction
     //! indicates if the init() function has been called
     bool isinit_;
 
+    //! time for evaluation of function
+    double time_;
+
     //! indicates if the setup() function has been called
     bool issetup_;
 

src/beaminteraction/src/4C_beaminteraction_submodel_evaluator_beamcontact.cpp

@@ -1103,7 +1103,7 @@ void BeamInteraction::SubmodelEvaluator::BeamContact::create_beam_contact_elemen
 
       for (auto& pair : newbeaminteractionpairs)
       {
-        pair->init(beam_contact_params_ptr_, ele_ptrs);
+        pair->init(beam_contact_params_ptr_, ele_ptrs, g_state().get_time_n());
         pair->setup();
 
         // add to list of current contact pairs

src/beaminteraction/src/contact/4C_beaminteraction_contact_beam_to_beam_input.cpp

@@ -42,6 +42,12 @@ std::vector<Core::IO::InputSpec> BeamInteraction::Contact::BeamToBeam::valid_par
           parameter<double>("BEAMS_BTBLINEPENALTYPARAM",
               {.description = "Penalty parameter per unit length for beam-to-beam line contact",
                   .default_value = -1.0}),
+          parameter<int>("BEAMS_BTBPENALTYPARAM_BY_FUNCT",
+              {.description = "Penalty parameter for beam-to-beam point contact by function id",
+                  .default_value = -1}),
+          parameter<int>("BEAMS_BTBLINEPENALTYPARAM_BY_FUNCT",
+              {.description = "Penalty parameter for beam-to-beam line contact by function id",
+                  .default_value = -1}),
           parameter<double>("BEAMS_PERPSHIFTANGLE1",
               {.description =
                       "Lower shift angle (in degrees) for penalty scaling of large-angle-contact",

tests/input_files/beam3r_line2_genalpha_liegroup_contact_line_penalty_by_function_linpen_twocrossedbeams.4C.yaml

@@ -0,0 +1,140 @@
+TITLE:
+  - "Checks if two overlapping beams"
+  - "can be displaced due to beam-to-beam-line-contact by increasing the beam-penalty"
+  - "value with the function specified from BEAMS_BTBPENALTYPARAM_BY_FUNCT."
+PROBLEM SIZE:
+  ELEMENTS: 2
+  NODES: 4
+  MATERIALS: 1
+  NUMDF: 6
+IO/RUNTIME VTK OUTPUT:
+  INTERVAL_STEPS: 1
+IO/RUNTIME VTK OUTPUT/STRUCTURE:
+  OUTPUT_STRUCTURE: true
+  DISPLACEMENT: true
+  ELEMENT_OWNER: true
+  ELEMENT_GID: true
+IO/RUNTIME VTK OUTPUT/BEAMS:
+  OUTPUT_BEAMS: true
+  DISPLACEMENT: true
+  STRAINS_GAUSSPOINT: true
+  ELEMENT_GID: true
+PROBLEM TYPE:
+  PROBLEMTYPE: "Structure"
+DISCRETISATION:
+  NUMFLUIDDIS: 0
+  NUMALEDIS: 0
+  NUMTHERMDIS: 0
+STRUCTURAL DYNAMIC:
+  INT_STRATEGY: "Standard"
+  DYNAMICTYPE: "GenAlphaLieGroup"
+  TIMESTEP: 0.1
+  NUMSTEP: 10
+  MAXTIME: 1
+  M_DAMP: 0
+  K_DAMP: 0
+  TOLDISP: 1e-08
+  MAXITER: 25
+  MASSLIN: "rotations"
+  LINEAR_SOLVER: 1
+STRUCTURAL DYNAMIC/GENALPHA:
+  BETA: 0.5
+  GAMMA: 1
+  ALPHA_M: 0.5
+  ALPHA_F: 0.5
+  RHO_INF: -1
+BEAM INTERACTION:
+  REPARTITIONSTRATEGY: "Everydt"
+  SEARCH_STRATEGY: bounding_volume_hierarchy
+BOUNDINGVOLUME STRATEGY:
+  WRITE_GEOMETRIC_SEARCH_VISUALIZATION: true
+  BEAM_RADIUS_EXTENSION_FACTOR: 11
+BEAM INTERACTION/BEAM TO BEAM CONTACT:
+  BEAMS_NEWGAP: false
+  BEAMS_BTBLINEPENALTYPARAM_BY_FUNCT: 2
+  BEAMS_BTBPENALTYPARAM: 0.0001
+  BEAMS_SEGCON: true
+  BEAMS_PERPSHIFTANGLE1: 70
+  BEAMS_PERPSHIFTANGLE2: 80
+  BEAMS_PARSHIFTANGLE1: 70
+  BEAMS_PARSHIFTANGLE2: 80
+  BEAMS_SEGANGLE: 10
+BINNING STRATEGY:
+  BIN_SIZE_LOWER_BOUND: 10
+  DOMAINBOUNDINGBOX: "-2 -2 -2 2 2 2"
+SOLVER 1:
+  SOLVER: "UMFPACK"
+  NAME: "Structure_Solver"
+DESIGN POINT DIRICH CONDITIONS:
+  - E: 1
+    NUMDOF: 6
+    ONOFF: [1, 1, 1, 1, 1, 1]
+    VAL: [0, 0, 0, 0, 0, 0]
+    FUNCT: [0, 0, 0, 0, 0, 0]
+  - E: 2
+    NUMDOF: 6
+    ONOFF: [1, 1, 1, 1, 1, 1]
+    VAL: [0, 0, 0, 0, 0, 0]
+    FUNCT: [0, 0, 0, 0, 0, 0]
+  - E: 3
+    NUMDOF: 6
+    ONOFF: [1, 1, 1, 1, 1, 1]
+    VAL: [0, 0, 0, 0, 0, 0]
+    FUNCT: [0, 0, 0, 0, 0, 0]
+BEAM INTERACTION/BEAM TO BEAM CONTACT CONDITIONS:
+  - E: 1
+    COUPLING_ID: 1
+  - E: 2
+    COUPLING_ID: 1
+DNODE-NODE TOPOLOGY:
+  - "NODE 1 DNODE 1"
+  - "NODE 2 DNODE 2"
+  - "NODE 3 DNODE 3"
+  - "NODE 4 DNODE 4"
+DLINE-NODE TOPOLOGY:
+  - "NODE 1 DLINE 1"
+  - "NODE 2 DLINE 1"
+  - "NODE 3 DLINE 2"
+  - "NODE 4 DLINE 2"
+NODE COORDS:
+  - "NODE 1 COORD -1.00000000000000e+00 0.000000000000000e+00 0.000000000000000e+00"
+  - "NODE 2 COORD 1.00000000000000e+00 0.0000000000000000e+00 0.000000000000000e+00"
+  - "NODE 3 COORD -1.00000000000000e+00 0.000000000000000e+00 0.1000000000000e+00"
+  - "NODE 4 COORD 1.00000000000000e+00 0.0000000000000000e+00 0.01000000000000e+00"
+STRUCTURE ELEMENTS:
+  - "1 BEAM3R LINE2 1 2 MAT 1 TRIADS 0 0 0 0 0 0"
+  - "2 BEAM3R LINE2 3 4 MAT 1 TRIADS 0 0 0 0 0 0"
+MATERIALS:
+  - MAT: 1
+    MAT_BeamReissnerElastHyper:
+      YOUNG: 100000
+      SHEARMOD: 50000
+      DENS: 0.001
+      CROSSAREA: 0.00282743
+      SHEARCORR: 1.1
+      MOMINPOL: 1.27235e-06
+      MOMIN2: 6.36173e-07
+      MOMIN3: 6.36173e-07
+FUNCT1:
+  - SYMBOLIC_FUNCTION_OF_TIME: "t"
+FUNCT2:
+  - SYMBOLIC_FUNCTION_OF_TIME: "10*t"
+RESULT DESCRIPTION:
+  - STRUCTURE:
+      DIS: "structure"
+      NODE: 4
+      QUANTITY: "dispx"
+      VALUE: 1.56653329972578407e-03
+      TOLERANCE: 1e-12
+  - STRUCTURE:
+      DIS: "structure"
+      NODE: 4
+      QUANTITY: "dispy"
+      VALUE: 0
+      TOLERANCE: 1e-12
+  - STRUCTURE:
+      DIS: "structure"
+      NODE: 4
+      QUANTITY: "dispz"
+      VALUE: 4.72233623507744282e-02
+      TOLERANCE: 1e-12