[beaminteraction]Dual Hermite shapefunctions for Lagrange Multipliers

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

@@ -13,6 +13,7 @@
 #include "4C_beaminteraction_contact_beam_to_solid_input.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_mortar_manager.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_mortar_manager_contact.hpp"
+#include "4C_beaminteraction_contact_beam_to_solid_mortar_shape_functions_dual_hermite.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_surface_meshtying_pair_gauss_point.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_surface_meshtying_pair_gauss_point_FAD.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_surface_meshtying_pair_mortar.hpp"
@@ -163,6 +164,10 @@ BeamInteraction::BeamToSolidCondition::create_indirect_assembly_manager(
       mortar_manager_params.n_lambda_node_translational = n_lambda_node_pos;
       mortar_manager_params.n_lambda_element_translational = n_lambda_element_pos;
 
+      mortar_manager_params.check_diagonal_d_matrix =
+          beam_to_solid_params_->get_mortar_shape_function_type() ==
+          BeamToSolid::BeamToSolidMortarShapefunctions::dual_hermite;
+
       if (beam_to_solid_params_->is_rotational_coupling())
       {
         // Get the mortar shape functions for rotational coupling
@@ -360,6 +365,19 @@ BeamInteraction::create_beam_to_solid_volume_pair_mortar(const Core::FE::CellTyp
     case BeamToSolid::BeamToSolidMortarShapefunctions::line4:
       return create_beam_to_solid_volume_pair_mortar<BtsClass, BtsMortarTemplateArguments...,
           GeometryPair::t_line4>(shape, other_mortar_shape_function...);
+    case BeamToSolid::BeamToSolidMortarShapefunctions::dual_hermite:
+    {
+      if constexpr (!BeamInteraction::is_rotation_pair_v<BtsClass>)
+      {
+        return create_beam_to_solid_volume_pair_mortar<BtsClass, BtsMortarTemplateArguments...,
+            BeamInteraction::HermiteDual>(shape, other_mortar_shape_function...);
+      }
+
+      FOUR_C_THROW(
+          "Dual Hermite mortar shape functions are currently not implemented for rotational "
+          "beam-to-solid volume meshtying.");
+      return nullptr;
+    }
     default:
       FOUR_C_THROW("Wrong mortar shape function.");
       return nullptr;
@@ -410,7 +428,7 @@ BeamInteraction::BeamToSolidConditionVolumeMeshtying::create_contact_pair_intern
     }
     else
     {
-      // Create the rotational mortart pairs.
+      // Create the rotational mortar pairs.
       return create_beam_to_solid_volume_pair_mortar<BeamToSolidVolumeMeshtyingPairMortarRotation>(
           shape, mortar_shape_function, mortar_shape_function_rotation);
     }

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

@@ -153,7 +153,9 @@ namespace BeamToSolid
     //! Quadratic.
     line3,
     //! Cubic.
-    line4
+    line4,
+    //! Dual shape functions for cubic Hermite interpolation
+    dual_hermite,
   };
 
   /**

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

@@ -548,6 +548,12 @@ void BeamInteraction::BeamToSolidMortarManager::evaluate_and_assemble_global_cou
   kappa_->complete();
   lambda_active_->complete();
   constraint_->complete();
+
+  if (parameters_.check_diagonal_d_matrix)
+  {
+    check_diagonal_like_structure(constraint_lin_beam_);
+    check_diagonal_like_structure(force_beam_lin_lambda_);
+  }
 }
 
 /**

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

@@ -77,6 +77,9 @@ namespace BeamInteraction
     BeamToSolid::BeamToSolidLagrangeFormulation lagrange_formulation =
         BeamToSolid::BeamToSolidLagrangeFormulation::none;
 
+    //! Flag to check whether the relevant D matrix has diagonal-like structure.
+    bool check_diagonal_d_matrix = false;
+
     //! Penalty parameter for positional coupling.
     double penalty_parameter_translational = 0.0;
 
@@ -284,6 +287,12 @@ namespace BeamInteraction
      */
     double get_energy() const;
 
+    /**
+     * \brief Check if there is only one nonzero entry per row of a matrix.
+     */
+    void check_diagonal_like_structure(
+        const std::shared_ptr<Core::LinAlg::SparseMatrix>& D_matrix) const;
+
     /**
      * \brief  Assembles the LM constraint into the global RHS vector
      */

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

@@ -0,0 +1,93 @@
+// This file is part of 4C multiphysics licensed under the
+// GNU Lesser General Public License v3.0 or later.
+//
+// See the LICENSE.md file in the top-level for license information.
+//
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+#ifndef FOUR_C_BEAMINTERACTION_CONTACT_BEAM_TO_SOLID_MORTAR_SHAPE_FUNCTIONS_DUAL_HERMITE_HPP
+#define FOUR_C_BEAMINTERACTION_CONTACT_BEAM_TO_SOLID_MORTAR_SHAPE_FUNCTIONS_DUAL_HERMITE_HPP
+
+#include "4C_config.hpp"
+
+#include "4C_beam3_base.hpp"
+#include "4C_fem_general_utils_fem_shapefunctions.hpp"
+#include "4C_geometry_pair_element.hpp"
+#include "4C_geometry_pair_element_shape_functions.hpp"
+
+FOUR_C_NAMESPACE_OPEN
+
+namespace BeamInteraction
+{
+  struct HermiteDual : public GeometryPair::ElementDiscretizationBase<Core::FE::CellType::line2, 2>
+  {
+  };
+}  // namespace BeamInteraction
+
+namespace GeometryPair
+{
+  template <>
+  struct ShapeFunctionData<BeamInteraction::HermiteDual>
+  {
+    double ref_length_;
+  };
+
+  template <>
+  struct SetShapeFunctionData<BeamInteraction::HermiteDual>
+  {
+    static void set(ShapeFunctionData<BeamInteraction::HermiteDual>& shape_function_data,
+        const Core::Elements::Element* element)
+    {
+      const auto* beam_element = dynamic_cast<const Discret::Elements::Beam3Base*>(element);
+      if (beam_element == nullptr)
+        FOUR_C_THROW(
+            "The element pointer has to point to a valid beam element when evaluating the shape "
+            "function data of a dual Hermite mortar shape function, as we need to get "
+            "RefLength()!");
+
+      shape_function_data.ref_length_ = beam_element->ref_length();
+    }
+  };
+
+
+  template <>
+  struct EvaluateShapeFunction<BeamInteraction::HermiteDual>
+  {
+    template <typename V, typename T>
+    static void evaluate(V& N, const T& xi,
+        const ShapeFunctionData<BeamInteraction::HermiteDual>& shape_function_data)
+    {
+      const T xi2 = xi * xi;
+      const T xi3 = xi2 * xi;
+
+      N(0) = -35.0 / 4.0 * xi3 + 15.0 / 4.0 * xi2 + 15.0 / 4.0 * xi - 3.0 / 4.0;
+
+      N(1) = (105.0 * xi3 - 45.0 / 2.0 * xi2 - 60.0 * xi + 15.0 / 2.0) /
+             shape_function_data.ref_length_;
+
+      N(2) = 35.0 / 4.0 * xi3 + 15.0 / 4.0 * xi2 - 15.0 / 4.0 * xi - 3.0 / 4.0;
+
+      N(3) = (105.0 * xi3 + 45.0 / 2.0 * xi2 - 60.0 * xi - 15.0 / 2.0) /
+             shape_function_data.ref_length_;
+    }
+  };
+
+  template <>
+  struct PrintElementData<BeamInteraction::HermiteDual>
+  {
+    template <typename ScalarType>
+    static void print(const ElementData<BeamInteraction::HermiteDual, ScalarType>& element_data,
+        std::ostream& out)
+    {
+      constexpr auto max_precision{std::numeric_limits<double>::digits10 + 1};
+      out << std::setprecision(max_precision);
+      out << "\nElement reference length: " << element_data.shape_function_data_.ref_length_;
+      out << "\nElement state vector: ";
+      element_data.element_position_.print(out);
+    }
+  };
+}  // namespace GeometryPair
+
+FOUR_C_NAMESPACE_CLOSE
+
+#endif

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

@@ -12,6 +12,7 @@
 #include "4C_beaminteraction_calc_utils.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_input.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_mortar_manager.hpp"
+#include "4C_beaminteraction_contact_beam_to_solid_mortar_shape_functions_dual_hermite.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_surface_params.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_volume_meshtying_params.hpp"
 #include "4C_beaminteraction_contact_pair.hpp"
@@ -130,6 +131,10 @@ BeamInteraction::mortar_shape_functions_to_number_of_lagrange_values(
       const int n_fourier_modes = beam_to_volume_parameters->get_number_of_fourier_modes();
       return (n_fourier_modes * 2 + 1);
     }
+    else if (shape_function == BeamToSolid::BeamToSolidMortarShapefunctions::dual_hermite)
+    {
+      return 2;
+    }
     else
     {
       return 1;
@@ -164,6 +169,12 @@ BeamInteraction::mortar_shape_functions_to_number_of_lagrange_values(
       const unsigned int n_lambda_element = 2 * n_dim * number_of_values_per_entry;
       return {n_lambda_node, n_lambda_element};
     }
+    case BeamToSolid::BeamToSolidMortarShapefunctions::dual_hermite:
+    {
+      const unsigned int n_lambda_node = n_dim * number_of_values_per_entry;
+      const unsigned int n_lambda_element = 0;
+      return {n_lambda_node, n_lambda_element};
+    }
     default:
       FOUR_C_THROW("Mortar shape function not implemented!");
   }
@@ -819,6 +830,51 @@ void BeamInteraction::assemble_local_mortar_contributions(
       Mortar::n_dof_, lambda_gid_pos.data(), local_kappa_active.data());
 }
 
+/**
+ *
+ */
+void BeamInteraction::BeamToSolidMortarManager::check_diagonal_like_structure(
+    const std::shared_ptr<Core::LinAlg::SparseMatrix>& D_matrix) const
+{
+  const double atol = 1.0e-14;
+  const double rtol = 1.0e-12;
+
+
+  for (int lid = 0; lid < D_matrix->row_map().num_my_elements(); ++lid)
+  {
+    const int row_gid = D_matrix->row_map().gid(lid);
+
+    int num_extracted = 0;
+    double* values = nullptr;
+    int* col_lids = nullptr;
+
+    D_matrix->extract_my_row_view(lid, num_extracted, values, col_lids);
+
+    if (num_extracted == 0 || num_extracted == 1)
+    {
+      continue;
+    }
+
+    int nnz_in_row = 0;
+    double row_scale = 0.0;
+
+    for (int k = 0; k < num_extracted; ++k)
+    {
+      row_scale = std::max(row_scale, std::abs(values[k]));
+      const double tol = atol + rtol * row_scale;
+      if (std::abs(values[k]) > tol)
+      {
+        ++nnz_in_row;
+      }
+    }
+
+    if (nnz_in_row != 1)
+    {
+      FOUR_C_THROW("Row {} has {} nonzero entries, but expected exactly one.", row_gid, nnz_in_row);
+    }
+  }
+}
+
 
 /**
  * Explicit template initialization of template functions.
@@ -920,50 +976,62 @@ namespace BeamInteraction
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex8, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex8, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex8, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex8, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex20, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex20, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex20, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex20, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex27, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex27, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex27, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_hex27, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tet4, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tet4, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tet4, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_tet4, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tet10, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tet10, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tet10, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_tet10, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_nurbs27, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_nurbs27, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_nurbs27, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_nurbs27, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad4, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad4, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad4, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad4, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad8, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad8, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad8, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad8, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad9, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad9, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad9, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_quad9, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tri3, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tri3, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tri3, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_tri3, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tri6, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tri6, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_tri6, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_tri6, HermiteDual);
 
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_nurbs9, t_line2);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_nurbs9, t_line3);
   initialize_template_assemble_local_mortar_contributions(t_hermite, t_nurbs9, t_line4);
+  initialize_template_assemble_local_mortar_contributions(t_hermite, t_nurbs9, HermiteDual);
 }  // namespace BeamInteraction
 
 FOUR_C_NAMESPACE_CLOSE

src/beaminteraction/src/contact/beam_to_solid/volume/4C_beaminteraction_contact_beam_to_solid_volume_meshtying_pair_mortar.cpp

@@ -8,6 +8,7 @@
 #include "4C_beaminteraction_contact_beam_to_solid_volume_meshtying_pair_mortar.hpp"
 
 #include "4C_beaminteraction_contact_beam_to_solid_mortar_manager.hpp"
+#include "4C_beaminteraction_contact_beam_to_solid_mortar_shape_functions_dual_hermite.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_utils.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_visualization_output_writer_base.hpp"
 #include "4C_beaminteraction_contact_beam_to_solid_visualization_output_writer_visualization.hpp"
@@ -112,7 +113,8 @@ void BeamInteraction::BeamToSolidVolumeMeshtyingPairMortar<Beam, Solid,
   if (visualization_discret != nullptr || visualization_continuous != nullptr)
   {
     // Setup variables.
-    GeometryPair::ElementData<Mortar, double> element_data_lambda;
+    auto element_data_lambda =
+        GeometryPair::InitializeElementData<Mortar, double>::initialize(this->element1());
     Core::LinAlg::Matrix<3, 1, scalar_type> X;
     Core::LinAlg::Matrix<3, 1, scalar_type> r;
     Core::LinAlg::Matrix<3, 1, scalar_type> u;
@@ -330,8 +332,10 @@ void BeamInteraction::BeamToSolidVolumeMeshtyingPairMortar<Beam, Solid, Mortar>:
       N_mortar.clear();
       N_beam.clear();
       N_solid.clear();
+      GeometryPair::ShapeFunctionData<Mortar> shape_function_data;
+      GeometryPair::SetShapeFunctionData<Mortar>::set(shape_function_data, this->element1());
       GeometryPair::EvaluateShapeFunction<Mortar>::evaluate(
-          N_mortar, projected_gauss_point.get_eta());
+          N_mortar, projected_gauss_point.get_eta(), shape_function_data);
       GeometryPair::EvaluateShapeFunction<Beam>::evaluate(
           N_beam, projected_gauss_point.get_eta(), this->ele1pos_.shape_function_data_);
       GeometryPair::EvaluateShapeFunction<Solid>::evaluate(
@@ -425,6 +429,13 @@ namespace BeamInteraction
   template class BeamToSolidVolumeMeshtyingPairMortar<t_hermite, t_tet4, t_line4>;
   template class BeamToSolidVolumeMeshtyingPairMortar<t_hermite, t_tet10, t_line4>;
   template class BeamToSolidVolumeMeshtyingPairMortar<t_hermite, t_nurbs27, t_line4>;
+
+  template class BeamToSolidVolumeMeshtyingPairMortar<t_hermite, t_hex8, HermiteDual>;
+  template class BeamToSolidVolumeMeshtyingPairMortar<t_hermite, t_hex20, HermiteDual>;
+  template class BeamToSolidVolumeMeshtyingPairMortar<t_hermite, t_hex27, HermiteDual>;
+  template class BeamToSolidVolumeMeshtyingPairMortar<t_hermite, t_tet4, HermiteDual>;
+  template class BeamToSolidVolumeMeshtyingPairMortar<t_hermite, t_tet10, HermiteDual>;
+  template class BeamToSolidVolumeMeshtyingPairMortar<t_hermite, t_nurbs27, HermiteDual>;
 }  // namespace BeamInteraction
 
 FOUR_C_NAMESPACE_CLOSE

src/beaminteraction/src/contact/beam_to_solid/volume/4C_beaminteraction_contact_beam_to_solid_volume_meshtying_pair_mortar_rotation.hpp

@@ -127,6 +127,17 @@ namespace BeamInteraction
         Core::LinAlg::Matrix<Solid::n_dof_, n_dof_rot_, double>& local_stiff_SB,
         Core::LinAlg::Matrix<Solid::n_dof_, Solid::n_dof_, double>& local_stiff_SS) const;
   };
+
+  /**
+   * \brief Compile time flag to check whether a beam-to-solid volume mortar pair includes
+   * rotational coupling
+   */
+  template <template <typename...> class T>
+  inline constexpr bool is_rotation_pair_v = false;
+  template <>
+  inline constexpr bool
+      is_rotation_pair_v<BeamInteraction::BeamToSolidVolumeMeshtyingPairMortarRotation> = true;
+
 }  // namespace BeamInteraction
 
 FOUR_C_NAMESPACE_CLOSE