diff --git a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_conditions.cpp b/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_conditions.cpp index 4548fb2f46e..064c330dab2 100644 --- a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_conditions.cpp +++ b/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(shape, other_mortar_shape_function...); + case BeamToSolid::BeamToSolidMortarShapefunctions::dual_hermite: + { + if constexpr (!BeamInteraction::is_rotation_pair_v) + { + return create_beam_to_solid_volume_pair_mortar(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( shape, mortar_shape_function, mortar_shape_function_rotation); } diff --git a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_input.hpp b/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_input.hpp index 7ec009dd727..f89f43ebafd 100644 --- a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_input.hpp +++ b/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, }; /** diff --git a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_mortar_manager.cpp b/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_mortar_manager.cpp index 371330c53af..b50e5883dcd 100644 --- a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_mortar_manager.cpp +++ b/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_); + } } /** diff --git a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_mortar_manager.hpp b/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_mortar_manager.hpp index 32d39f5c59a..c3078644662 100644 --- a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_mortar_manager.hpp +++ b/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& D_matrix) const; + /** * \brief Assembles the LM constraint into the global RHS vector */ diff --git a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_mortar_shape_functions_dual_hermite.hpp b/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_mortar_shape_functions_dual_hermite.hpp new file mode 100644 index 00000000000..807e58eacff --- /dev/null +++ b/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 + { + }; +} // namespace BeamInteraction + +namespace GeometryPair +{ + template <> + struct ShapeFunctionData + { + double ref_length_; + }; + + template <> + struct SetShapeFunctionData + { + static void set(ShapeFunctionData& shape_function_data, + const Core::Elements::Element* element) + { + const auto* beam_element = dynamic_cast(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 + { + template + static void evaluate(V& N, const T& xi, + const ShapeFunctionData& 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 + { + template + static void print(const ElementData& element_data, + std::ostream& out) + { + constexpr auto max_precision{std::numeric_limits::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 \ No newline at end of file diff --git a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_utils.cpp b/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_utils.cpp index dd286be8deb..04de95e5ff2 100644 --- a/src/beaminteraction/src/contact/beam_to_solid/4C_beaminteraction_contact_beam_to_solid_utils.cpp +++ b/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& 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 diff --git a/src/beaminteraction/src/contact/beam_to_solid/volume/4C_beaminteraction_contact_beam_to_solid_volume_meshtying_pair_mortar.cpp b/src/beaminteraction/src/contact/beam_to_solid/volume/4C_beaminteraction_contact_beam_to_solid_volume_meshtying_pair_mortar.cpp index 15206675b22..b8f6d4b6426 100644 --- a/src/beaminteraction/src/contact/beam_to_solid/volume/4C_beaminteraction_contact_beam_to_solid_volume_meshtying_pair_mortar.cpp +++ b/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 element_data_lambda; + auto element_data_lambda = + GeometryPair::InitializeElementData::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: N_mortar.clear(); N_beam.clear(); N_solid.clear(); + GeometryPair::ShapeFunctionData shape_function_data; + GeometryPair::SetShapeFunctionData::set(shape_function_data, this->element1()); GeometryPair::EvaluateShapeFunction::evaluate( - N_mortar, projected_gauss_point.get_eta()); + N_mortar, projected_gauss_point.get_eta(), shape_function_data); GeometryPair::EvaluateShapeFunction::evaluate( N_beam, projected_gauss_point.get_eta(), this->ele1pos_.shape_function_data_); GeometryPair::EvaluateShapeFunction::evaluate( @@ -425,6 +429,13 @@ namespace BeamInteraction template class BeamToSolidVolumeMeshtyingPairMortar; template class BeamToSolidVolumeMeshtyingPairMortar; template class BeamToSolidVolumeMeshtyingPairMortar; + + template class BeamToSolidVolumeMeshtyingPairMortar; + template class BeamToSolidVolumeMeshtyingPairMortar; + template class BeamToSolidVolumeMeshtyingPairMortar; + template class BeamToSolidVolumeMeshtyingPairMortar; + template class BeamToSolidVolumeMeshtyingPairMortar; + template class BeamToSolidVolumeMeshtyingPairMortar; } // namespace BeamInteraction FOUR_C_NAMESPACE_CLOSE diff --git a/src/beaminteraction/src/contact/beam_to_solid/volume/4C_beaminteraction_contact_beam_to_solid_volume_meshtying_pair_mortar_rotation.hpp b/src/beaminteraction/src/contact/beam_to_solid/volume/4C_beaminteraction_contact_beam_to_solid_volume_meshtying_pair_mortar_rotation.hpp index 94e820ad425..872d4e9bee1 100644 --- a/src/beaminteraction/src/contact/beam_to_solid/volume/4C_beaminteraction_contact_beam_to_solid_volume_meshtying_pair_mortar_rotation.hpp +++ b/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& local_stiff_SB, Core::LinAlg::Matrix& local_stiff_SS) const; }; + + /** + * \brief Compile time flag to check whether a beam-to-solid volume mortar pair includes + * rotational coupling + */ + template