diff --git a/packages/contact_mechanics.cmake b/packages/contact_mechanics.cmake index 98609cf76..874b27225 100644 --- a/packages/contact_mechanics.cmake +++ b/packages/contact_mechanics.cmake @@ -1,75 +1,76 @@ #=============================================================================== # @file contact_mechanics.cmake # # @author Mohit Pundir # # @date creation: Fri Sep 03 2010 # @date last modification: Wed Jun 23 2021 # # @brief package description for contact mechanics # # # @section LICENSE # # Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License along # with Akantu. If not, see . # #=============================================================================== package_declare(contact_mechanics DEPENDS model_couplers cohesive_element DESCRIPTION "Use Contact Mechanics package of Akantu") package_declare_sources(contact_mechanics model/contact_mechanics/contact_mechanics_model.hh model/contact_mechanics/contact_mechanics_model.cc model/contact_mechanics/contact_detector.hh model/contact_mechanics/contact_detector.cc model/contact_mechanics/contact_detector_inline_impl.cc model/contact_mechanics/contact_element.hh model/contact_mechanics/geometry_utils.hh model/contact_mechanics/geometry_utils.cc model/contact_mechanics/geometry_utils_inline_impl.cc model/contact_mechanics/resolution.hh model/contact_mechanics/resolution.cc model/contact_mechanics/resolution_utils.hh model/contact_mechanics/resolution_utils.cc model/contact_mechanics/resolutions/resolution_penalty.hh + model/contact_mechanics/resolutions/resolution_penalty_tmpl.hh model/contact_mechanics/resolutions/resolution_penalty.cc model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc model/contact_mechanics/surface_selector.hh model/contact_mechanics/surface_selector.cc model/model_couplers/coupler_solid_contact.hh model/model_couplers/coupler_solid_contact_tmpl.hh model/model_couplers/coupler_solid_contact.cc model/model_couplers/coupler_solid_cohesive_contact.hh model/model_couplers/coupler_solid_cohesive_contact.cc model/model_couplers/cohesive_contact_solvercallback.hh model/model_couplers/cohesive_contact_solvercallback.cc ) package_declare_documentation_files(contact_mechanics manual-contactmechanicsmodel.tex manual-contact-detector.tex ) package_declare_documentation(contact_mechanics "This package activates the contact mechanics model") diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc index 688b4db6e..8b40b6c93 100644 --- a/src/model/contact_mechanics/resolution.cc +++ b/src/model/contact_mechanics/resolution.cc @@ -1,222 +1,222 @@ /** * @file resolution.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 - * @date last modification: Wed Apr 07 2021 + * @date last modification: Fri Oct 21 2022 * * @brief Implementation of common part of the contact resolution class * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "resolution.hh" #include "contact_mechanics_model.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Resolution::Resolution(ContactMechanicsModel & model, const ID & id) : Parsable(ParserType::_contact_resolution, id), id(id), fem(model.getFEEngine()), model(model) { AKANTU_DEBUG_IN(); spatial_dimension = model.getMesh().getSpatialDimension(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Resolution::~Resolution() = default; /* -------------------------------------------------------------------------- */ void Resolution::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("mu", mu, Real(0.), _pat_parsable | _pat_modifiable, "Friction Coefficient"); registerParam("is_master_deformable", is_master_deformable, bool(false), _pat_parsable | _pat_readable, "Is master surface deformable"); } /* -------------------------------------------------------------------------- */ void Resolution::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "Contact Resolution " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Resolution::assembleInternalForces(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); this->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleInternalForces() { AKANTU_DEBUG_IN(); for (const auto & element : model.getContactElements()) { auto nb_nodes = element.getNbNodes(); Vector local_fn(nb_nodes * spatial_dimension); computeNormalForce(element, local_fn); Vector local_ft(nb_nodes * spatial_dimension); computeTangentialForce(element, local_ft); Vector local_fc(nb_nodes * spatial_dimension); local_fc = local_fn + local_ft; assembleLocalToGlobalArray(element, local_fn, model.getNormalForce()); assembleLocalToGlobalArray(element, local_ft, model.getTangentialForce()); assembleLocalToGlobalArray(element, local_fc, model.getInternalForce()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleLocalToGlobalArray(const ContactElement & element, Vector & local, Array & global) { auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn = const_cast(model.getMesh()).getConnectivity(master); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i - 1]; } return elem_conn; }; auto & surface_selector = model.getContactDetector().getSurfaceSelector(); auto & slave_list = surface_selector.getSlaveList(); auto & master_list = surface_selector.getMasterList(); auto connectivity = get_connectivity(element.slave, element.master); UInt nb_dofs = global.getNbComponent(); UInt nb_nodes = is_master_deformable ? connectivity.size() : 1; Real alpha = is_master_deformable ? 0.5 : 1.; for (UInt i : arange(nb_nodes)) { UInt n = connectivity[i]; auto slave_result = std::find(slave_list.begin(), slave_list.end(), n); auto master_result = std::find(master_list.begin(), master_list.end(), n); for (UInt j : arange(nb_dofs)) { UInt offset_node = n * nb_dofs + j; global[offset_node] += alpha * local[i * nb_dofs + j]; } } } /* -------------------------------------------------------------------------- */ void Resolution::assembleStiffnessMatrix(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); auto & global_stiffness = const_cast(model.getDOFManager().getMatrix("K")); for (const auto & element : model.getContactElements()) { auto nb_nodes = element.getNbNodes(); Matrix local_kn(nb_nodes * spatial_dimension, nb_nodes * spatial_dimension); computeNormalModuli(element, local_kn); assembleLocalToGlobalMatrix(element, local_kn, global_stiffness); Matrix local_kt(nb_nodes * spatial_dimension, nb_nodes * spatial_dimension); computeTangentialModuli(element, local_kt); assembleLocalToGlobalMatrix(element, local_kt, global_stiffness); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleLocalToGlobalMatrix(const ContactElement & element, const Matrix & local, SparseMatrix & global) { auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn = const_cast(model.getMesh()).getConnectivity(master); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i - 1]; } return elem_conn; }; auto connectivity = get_connectivity(element.slave, element.master); auto nb_dofs = spatial_dimension; UInt nb_nodes = is_master_deformable ? connectivity.size() : 1; UInt total_nb_dofs = nb_dofs * nb_nodes; std::vector equations; for (UInt i : arange(connectivity.size())) { UInt conn = connectivity[i]; for (UInt j : arange(nb_dofs)) { equations.push_back(conn * nb_dofs + j); } } for (UInt i : arange(total_nb_dofs)) { UInt row = equations[i]; for (UInt j : arange(total_nb_dofs)) { UInt col = equations[j]; global.add(row, col, local(i, j)); } } } /* -------------------------------------------------------------------------- */ void Resolution::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void Resolution::afterSolveStep(__attribute__((unused)) bool converged) {} } // namespace akantu diff --git a/src/model/contact_mechanics/resolution.hh b/src/model/contact_mechanics/resolution.hh index eaef3d209..9877adcf7 100644 --- a/src/model/contact_mechanics/resolution.hh +++ b/src/model/contact_mechanics/resolution.hh @@ -1,236 +1,242 @@ /** * @file resolution.hh * * @author Mohit Pundir * * @date creation: Fri Jun 18 2010 - * @date last modification: Wed Apr 07 2021 + * @date last modification: Fri Oct 21 2022 * * @brief Mother class for all contact resolutions * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_factory.hh" #include "contact_element.hh" #include "fe_engine.hh" #include "geometry_utils.hh" #include "parsable.hh" #include "parser.hh" #include "resolution_utils.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RESOLUTION_HH__ #define __AKANTU_RESOLUTION_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class ContactMechanicsModel; } // namespace akantu namespace akantu { /** * Interface of all contact resolutions * Prerequisites for a new resolution * - inherit from this class * - implement the following methods: * \code * * virtual void computeNormalForce(); * virtual void computeTangentialForce(); * virtual void computeNormalModuli(); * virtual void computeTangentialModuli(); * * \endcode * */ class Resolution : public Parsable { /* ------------------------------------------------------------------------ */ /* Constructor/Destructor */ /* ------------------------------------------------------------------------ */ public: /// instantiate contact resolution with defaults Resolution(ContactMechanicsModel & model, const ID & id = ""); /// Destructor ~Resolution() override; protected: void initialize(); /// computes coordinates of a given element void computeCoordinates(const Element &, Matrix &); /* ------------------------------------------------------------------------ */ /* Functions that resolutions should reimplement for force */ /* ------------------------------------------------------------------------ */ public: /// computes the force vector due to normal traction virtual void computeNormalForce(__attribute__((unused)) const ContactElement & /*unused*/, __attribute__((unused)) Vector & /*unused*/) { AKANTU_TO_IMPLEMENT(); } /// computes the tangential force vector due to frictional traction virtual void computeTangentialForce(__attribute__((unused)) const ContactElement & /*unused*/, __attribute__((unused)) Vector & /*unused*/) { AKANTU_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* Functions that resolutions should reimplement for stiffness */ /* ------------------------------------------------------------------------ */ public: /// compute the normal moduli due to normal traction virtual void computeNormalModuli(__attribute__((unused)) const ContactElement & /*unused*/, __attribute__((unused)) Matrix & /*unused*/) { AKANTU_TO_IMPLEMENT(); } /// compute the tangent moduli due to tangential traction virtual void computeTangentialModuli(__attribute__((unused)) const ContactElement & /*unused*/, __attribute__((unused)) Matrix & /*unused*/) { AKANTU_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// assemble the residual for this resolution void assembleInternalForces(GhostType ghost_type); /// assemble the stiffness matrix for this resolution void assembleStiffnessMatrix(GhostType ghost_type); private: /// assemble the residual for this resolution void assembleInternalForces(); /// assemble the local array to global array for a contact element void assembleLocalToGlobalArray(const ContactElement & /*element*/, Vector & /*local*/, Array & /*global*/); /// assemble the local stiffness to global stiffness for a contact element void assembleLocalToGlobalMatrix(const ContactElement & /*element*/, const Matrix & /*local*/, SparseMatrix & /*global*/); public: virtual void beforeSolveStep(); virtual void afterSolveStep(bool converged = true); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const ID &); public: /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: ID id; /// friction coefficient : mu Real mu; /// spatial dimension UInt spatial_dimension; /// is master surface deformable bool is_master_deformable; /// Link to the fe engine object in the model FEEngine & fem; /// resolution name std::string name; /// model to which the resolution belong ContactMechanicsModel & model; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Resolution & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { using ResolutionFactory = Factory; /// macaulay bracket to convert positive gap to zero template T macaulay(T var) { return var < 0 ? 0 : var; } template T heaviside(T var) { return var < 0 ? 0 : 1.0; } } // namespace akantu -#define INSTANTIATE_RESOLUTION_ONLY(res_name) class res_name +#define INSTANTIATE_RESOLUTION_ONLY(res_name, dtype, pen_fun, pen_der) \ + class res_name, pen_der> -#define RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name) \ +#define RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR( \ + id, res_name, dtype, pen_fun, pen_der)\ [](UInt dim, const ID &, ContactMechanicsModel & model, \ const ID & id) -> std::unique_ptr { \ switch (dim) { \ case 1: \ - return std::make_unique(model, id); \ + return std::make_unique, pen_der>>(\ + model, id);\ case 2: \ - return std::make_unique(model, id); \ + return std::make_unique, pen_der>>(\ + model, id);\ case 3: \ - return std::make_unique(model, id); \ + return std::make_unique, pen_der>>(\ + model, id);\ default: \ AKANTU_EXCEPTION( \ "The dimension " \ << dim << "is not a valid dimension for the contact resolution " \ << #id); \ } \ } -#define INSTANTIATE_RESOLUTION(id, res_name) \ - INSTANTIATE_RESOLUTION_ONLY(res_name); \ +#define INSTANTIATE_RESOLUTION(id, res_name, dtype, pen_fun, pen_der) \ + INSTANTIATE_RESOLUTION_ONLY(res_name, dtype, pen_fun, pen_der); \ static bool resolution_is_alocated_##id [[gnu::unused]] = \ ResolutionFactory::getInstance().registerAllocator( \ - #id, RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name)) + #id, RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR( \ + id, res_name, dtype,pen_fun, pen_der)) #endif /* __AKANTU_RESOLUTION_HH__ */ diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc index 35ad7330b..67f85d8d8 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc @@ -1,868 +1,39 @@ /** * @file resolution_penalty.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 - * @date last modification: Wed Jun 09 2021 + * @date last modification: Fri Oct 21 2022 * * @brief Specialization of the resolution class for the penalty method * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ -#include "resolution_penalty.hh" +#include "resolution_penalty_tmpl.hh" #include "element_class_helper.hh" /* -------------------------------------------------------------------------- */ namespace akantu { -/* -------------------------------------------------------------------------- */ -ResolutionPenalty::ResolutionPenalty(ContactMechanicsModel & model, - const ID & id, PenaltyType penalty_type) - : Resolution(model, id) { - AKANTU_DEBUG_IN(); - this->initialize(); - AKANTU_DEBUG_OUT(); - this->penalty_type = penalty_type; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::initialize() { - this->registerParam("epsilon_n", epsilon_n, Real(0.), - _pat_parsable | _pat_modifiable, - "Normal penalty parameter"); - this->registerParam("epsilon_t", epsilon_t, Real(0.), - _pat_parsable | _pat_modifiable, - "Tangential penalty parameter"); -} - -/* -------------------------------------------------------------------------- */ -Real ResolutionPenalty::computeNormalTraction(Real & gap) const { - switch(this->penalty_type) { - case penalty_linear: - return epsilon_n * macaulay(gap); - case penalty_quadratic: - return epsilon_n * (macaulay(gap) * macaulay(gap) + macaulay(gap)); - default: AKANTU_ERROR("This case should not happen!!"); - } -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeNormalForce(const ContactElement & element, - Vector & force) { - - force.zero(); - - auto & gaps = model.getGaps(); - auto & projections = model.getProjections(); - auto & normals = model.getNormals(); - - auto surface_dimension = spatial_dimension - 1; - - Real gap(gaps.begin()[element.slave]); - Vector normal(normals.begin(spatial_dimension)[element.slave]); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & nodal_area = const_cast &>(model.getNodalArea()); - - // compute normal traction - Real p_n = computeNormalTraction(gap); - p_n *= nodal_area[element.slave]; - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix shape_matric(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - ResolutionUtils::computeShapeFunctionMatric(element, projection, - shape_matric); - - force.mul(shape_matric, normal, p_n); -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentialForce(const ContactElement & element, - Vector & force) { - - if (mu == 0) { - return; - } - - force.zero(); - - UInt surface_dimension = spatial_dimension - 1; - - // compute covariant basis - auto & projections = model.getProjections(); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & normals = model.getNormals(); - Vector normal(normals.begin(spatial_dimension)[element.slave]); - - auto & tangents = model.getTangents(); - Matrix covariant_basis( - tangents.begin(surface_dimension, spatial_dimension)[element.slave]); - - // check for no-contact to contact condition - // need a better way to check if new node added is not presnt in the - // previous master elemets - auto & previous_master_elements = model.getPreviousMasterElements(); - if (element.slave >= previous_master_elements.size()) { - return; - } - - auto & previous_element = previous_master_elements[element.slave]; - if (previous_element.type == _not_defined) { - return; - } - - // compute tangential traction using return map algorithm - auto & tangential_tractions = model.getTangentialTractions(); - Vector tangential_traction( - tangential_tractions.begin(surface_dimension)[element.slave]); - this->computeTangentialTraction(element, covariant_basis, - tangential_traction); - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix shape_matric(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - ResolutionUtils::computeShapeFunctionMatric(element, projection, - shape_matric); - - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - auto & nodal_area = const_cast &>(model.getNodalArea()); - - for (auto && values1 : enumerate(covariant_basis.transpose())) { - auto & alpha = std::get<0>(values1); - auto & tangent_alpha = std::get<1>(values1); - for (auto && values2 : enumerate(tangential_traction)) { - auto & beta = std::get<0>(values2); - auto & traction_beta = std::get<1>(values2); - Vector tmp(force.size()); - tmp.mul(shape_matric, tangent_alpha, traction_beta); - tmp *= - contravariant_metric_tensor(alpha, beta) * nodal_area[element.slave]; - force += tmp; - } - } -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentialTraction( - const ContactElement & element, const Matrix & covariant_basis, - Vector & traction_tangential) { - - UInt surface_dimension = spatial_dimension - 1; - - auto & gaps = model.getGaps(); - auto & gap = gaps.begin()[element.slave]; - - // Return map algorithm is employed - // compute trial traction - Vector traction_trial(surface_dimension); - this->computeTrialTangentialTraction(element, covariant_basis, - traction_trial); - - // compute norm of trial traction - Real traction_trial_norm = 0; - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - for (auto i : arange(surface_dimension)) { - for (auto j : arange(surface_dimension)) { - traction_trial_norm += traction_trial[i] * traction_trial[j] * - contravariant_metric_tensor(i, j); - } - } - traction_trial_norm = sqrt(traction_trial_norm); - - // check stick or slip condition - auto & contact_state = model.getContactState(); - auto & state = contact_state.begin()[element.slave]; - - Real p_n = computeNormalTraction(gap); - bool stick = (traction_trial_norm <= mu * p_n); - - if (stick) { - state = ContactState::_stick; - computeStickTangentialTraction(element, traction_trial, - traction_tangential); - } else { - state = ContactState::_slip; - computeSlipTangentialTraction(element, covariant_basis, traction_trial, - traction_tangential); - } -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTrialTangentialTraction( - const ContactElement & element, const Matrix & covariant_basis, - Vector & traction) { - - UInt surface_dimension = spatial_dimension - 1; - - auto & projections = model.getProjections(); - Vector current_projection( - projections.begin(surface_dimension)[element.slave]); - - auto & previous_projections = model.getPreviousProjections(); - Vector previous_projection( - previous_projections.begin(surface_dimension)[element.slave]); - - // method from Laursen et. al. - /*auto covariant_metric_tensor = - GeometryUtils::covariantMetricTensor(covariant_basis); auto - increment_projection = current_projection - previous_projection; - - traction.mul(covariant_metric_tensor, increment_projection, epsilon_t); - - auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); - Vector - previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]); - traction = previous_traction + traction;*/ - - // method from Schweizerhof - auto covariant_metric_tensor = - GeometryUtils::covariantMetricTensor(covariant_basis); - - auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); - Vector previous_traction( - previous_tangential_tractions.begin(surface_dimension)[element.slave]); - - auto & previous_tangents = model.getPreviousTangents(); - Matrix previous_covariant_basis(previous_tangents.begin( - surface_dimension, spatial_dimension)[element.slave]); - auto previous_contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(previous_covariant_basis); - - auto current_tangent = covariant_basis.transpose(); - auto previous_tangent = previous_covariant_basis.transpose(); - - for (auto alpha : arange(surface_dimension)) { - Vector tangent_alpha(current_tangent(alpha)); - for (auto gamma : arange(surface_dimension)) { - for (auto beta : arange(surface_dimension)) { - Vector tangent_beta(previous_tangent(beta)); - auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha); - traction[alpha] += previous_traction[gamma] * - previous_contravariant_metric_tensor(gamma, beta) * - t_alpha_t_beta; - } - } - } - - auto & previous_master_elements = model.getPreviousMasterElements(); - auto & previous_element = previous_master_elements[element.slave]; - - Vector previous_real_projection(spatial_dimension); - GeometryUtils::realProjection( - model.getMesh(), model.getContactDetector().getPositions(), - previous_element, previous_projection, previous_real_projection); - - Vector current_real_projection(spatial_dimension); - GeometryUtils::realProjection( - model.getMesh(), model.getContactDetector().getPositions(), - element.master, current_projection, current_real_projection); - - auto increment_real = current_real_projection - previous_real_projection; - Vector increment_xi(surface_dimension); - - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - // increment in natural coordinate - for (auto beta : arange(surface_dimension)) { - for (auto gamma : arange(surface_dimension)) { - auto temp = increment_real.dot(current_tangent(gamma)); - temp *= contravariant_metric_tensor(beta, gamma); - increment_xi[beta] += temp; - } - } - - Vector temp(surface_dimension); - temp.mul(covariant_metric_tensor, increment_xi, epsilon_t); - - traction -= temp; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeStickTangentialTraction( - const ContactElement & /*element*/, Vector & traction_trial, - Vector & traction_tangential) { - traction_tangential = traction_trial; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeSlipTangentialTraction( - const ContactElement & element, const Matrix & covariant_basis, - Vector & traction_trial, Vector & traction_tangential) { - UInt surface_dimension = spatial_dimension - 1; - - auto & gaps = model.getGaps(); - auto & gap = gaps.begin()[element.slave]; - - // compute norm of trial traction - Real traction_trial_norm = 0; - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - for (auto alpha : arange(surface_dimension)) { - for (auto beta : arange(surface_dimension)) { - traction_trial_norm += traction_trial[alpha] * traction_trial[beta] * - contravariant_metric_tensor(alpha, beta); - } - } - traction_trial_norm = sqrt(traction_trial_norm); - - auto slip_direction = traction_trial; - slip_direction /= traction_trial_norm; - - Real p_n = computeNormalTraction(gap); - traction_tangential = slip_direction; - traction_tangential *= mu * p_n; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeNormalModuli(const ContactElement & element, - Matrix & stiffness) { - - auto surface_dimension = spatial_dimension - 1; - - auto & gaps = model.getGaps(); - Real gap(gaps.begin()[element.slave]); - - auto & projections = model.getProjections(); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & nodal_areas = model.getNodalArea(); - auto & nodal_area = nodal_areas.begin()[element.slave]; - - auto & normals = model.getNormals(); - Vector normal(normals.begin(spatial_dimension)[element.slave]); - - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof - // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 - - // construct A matrix - const ElementType & type = element.master.type; - auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - A(j, i * spatial_dimension + j) = 1; - continue; - } - A(j, i * spatial_dimension + j) = -shapes[i - 1]; - } - } - - // construct the main part of normal matrix - Matrix k_main(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - Matrix n_outer_n(spatial_dimension, spatial_dimension); - Matrix mat_n(normal.storage(), normal.size(), 1.); - n_outer_n.mul(mat_n, mat_n); - - Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - tmp.mul(n_outer_n, A); - - k_main.mul(A, tmp); - switch(this->penalty_type) { - case penalty_linear: - k_main *= epsilon_n * heaviside(gap) * nodal_area; - break; - case penalty_quadratic: - k_main *= epsilon_n * heaviside(gap) * (2 * gap + 1) * nodal_area; - break; - default: AKANTU_ERROR("This case should not happen!!"); - } - - // construct the rotational part of the normal matrix - auto & tangents = model.getTangents(); - Matrix covariant_basis( - tangents.begin(surface_dimension, spatial_dimension)[element.slave]); - - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - // computing shape derivatives - auto && shape_derivatives = - ElementClassHelper<_ek_regular>::getDNDS(projection, type); - - // consists of 2 rotational parts - Matrix k_rot1(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix k_rot2(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - auto construct_Aj = [&](auto && dnds) { - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - Aj(j, i * spatial_dimension + j) = 0; - continue; - } - Aj(j, i * spatial_dimension + j) = dnds(i - 1); - } - } - }; - - for (auto && values1 : enumerate(covariant_basis.transpose())) { - auto & alpha = std::get<0>(values1); - auto & tangent = std::get<1>(values1); - - Matrix n_outer_t(spatial_dimension, spatial_dimension); - Matrix mat_t(tangent.storage(), tangent.size(), 1.); - n_outer_t.mul(mat_n, mat_t); - - Matrix t_outer_n(spatial_dimension, spatial_dimension); - t_outer_n.mul(mat_t, mat_n); - - for (auto && values2 : enumerate(shape_derivatives.transpose())) { - auto & beta = std::get<0>(values2); - auto & dnds = std::get<1>(values2); - // construct Aj from shape function wrt to jth natural - // coordinate - construct_Aj(dnds); - - Matrix tmp(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp1(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp.mul(n_outer_t, A); - tmp1.mul(Aj, tmp); - tmp1 *= contravariant_metric_tensor(alpha, beta); - k_rot1 += tmp1; - - tmp.mul(t_outer_n, Aj); - tmp1.mul(A, tmp); - tmp1 *= contravariant_metric_tensor(alpha, beta); - k_rot2 += tmp1; - } - } - - switch(this->penalty_type) { - case penalty_linear: - k_rot1 *= -epsilon_n * heaviside(gap) * gap * nodal_area; - k_rot2 *= -epsilon_n * heaviside(gap) * gap * nodal_area; - break; - case penalty_quadratic: - k_rot1 *= -epsilon_n * heaviside(gap) * (gap * gap + gap) * nodal_area; - k_rot2 *= -epsilon_n * heaviside(gap) * (gap * gap + gap) * nodal_area; - break; - default: AKANTU_ERROR("This case should not happen!!"); - } - - stiffness += k_main + k_rot1 + k_rot2; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentialModuli(const ContactElement & element, - Matrix & stiffness) { - if (mu == 0) { - return; - } - - stiffness.zero(); - - auto & contact_state = model.getContactState(); - auto state = contact_state.begin()[element.slave]; - - switch (state) { - case ContactState::_stick: { - computeStickModuli(element, stiffness); - break; - } - case ContactState::_slip: { - computeSlipModuli(element, stiffness); - break; - } - default: - break; - } -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeStickModuli(const ContactElement & element, - Matrix & stiffness) { - - auto surface_dimension = spatial_dimension - 1; - - auto & projections = model.getProjections(); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & nodal_areas = model.getNodalArea(); - auto & nodal_area = nodal_areas.begin()[element.slave]; - - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof - // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 - - // construct A matrix - const ElementType & type = element.master.type; - auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - A(j, i * spatial_dimension + j) = 1; - continue; - } - A(j, i * spatial_dimension + j) = -shapes[i - 1]; - } - } - - // computing shape derivatives - auto && shape_derivatives = - ElementClassHelper<_ek_regular>::getDNDS(projection, type); - - Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - auto construct_Aj = [&](auto && dnds) { - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - Aj(j, i * spatial_dimension + j) = 0; - continue; - } - Aj(j, i * spatial_dimension + j) = dnds(i - 1); - } - } - }; - - // tangents should have been calculated in normal modulii - auto & tangents = model.getTangents(); - Matrix covariant_basis( - tangents.begin(surface_dimension, spatial_dimension)[element.slave]); - - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - // construct 1st part of the stick modulii - Matrix k_main(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - for (auto && values1 : enumerate(covariant_basis.transpose())) { - auto & alpha = std::get<0>(values1); - auto & tangent_alpha = std::get<1>(values1); - - Matrix t_outer_t(spatial_dimension, spatial_dimension); - Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); - - for (auto && values2 : enumerate(covariant_basis.transpose())) { - auto & beta = std::get<0>(values2); - auto & tangent_beta = std::get<1>(values2); - - Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); - t_outer_t.mul(mat_t_alpha, mat_t_beta); - - Matrix tmp(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp1(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp.mul(t_outer_t, A); - tmp1.mul(A, tmp); - tmp1 *= contravariant_metric_tensor(alpha, beta); - k_main += tmp1; - } - } - - k_main *= -epsilon_t; - - // construct 2nd part of the stick modulii - auto & tangential_tractions = model.getTangentialTractions(); - Vector tangential_traction( - tangential_tractions.begin(surface_dimension)[element.slave]); - - Matrix k_second(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - for (auto alpha : arange(surface_dimension)) { - - Matrix k_sum(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - for (auto && values1 : enumerate(shape_derivatives.transpose())) { - auto & beta = std::get<0>(values1); - auto & dnds = std::get<1>(values1); - // construct Aj from shape function wrt to jth natural - // coordinate - construct_Aj(dnds); - for (auto && values2 : enumerate(covariant_basis.transpose())) { - auto & gamma = std::get<0>(values2); - auto & tangent_gamma = std::get<1>(values2); - - Matrix t_outer_t(spatial_dimension, spatial_dimension); - Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), - 1.); - - for (auto && values3 : enumerate(covariant_basis.transpose())) { - auto & theta = std::get<0>(values3); - auto & tangent_theta = std::get<1>(values3); - - Matrix mat_t_theta(tangent_theta.storage(), - tangent_theta.size(), 1.); - t_outer_t.mul(mat_t_gamma, mat_t_theta); - - Matrix tmp(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp1(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp.mul(t_outer_t, Aj); - tmp1.mul(A, tmp); - tmp1 *= contravariant_metric_tensor(alpha, theta) * - contravariant_metric_tensor(beta, gamma); - - Matrix tmp2(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp3(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp2.mul(t_outer_t, A); - tmp3.mul(Aj, tmp2); - tmp3 *= contravariant_metric_tensor(alpha, gamma) * - contravariant_metric_tensor(beta, theta); - - k_sum += tmp1 + tmp3; - } - } - } - - k_second += tangential_traction[alpha] * k_sum; - } - - stiffness += k_main * nodal_area - k_second * nodal_area; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeSlipModuli(const ContactElement & element, - Matrix & stiffness) { - - auto surface_dimension = spatial_dimension - 1; - - auto & gaps = model.getGaps(); - Real gap(gaps.begin()[element.slave]); - - auto & nodal_areas = model.getNodalArea(); - auto & nodal_area = nodal_areas.begin()[element.slave]; - - // compute normal traction - Real p_n = computeNormalTraction(gap); - - auto & projections = model.getProjections(); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & normals = model.getNormals(); - Vector normal(normals.begin(spatial_dimension)[element.slave]); - - // restructure normal as a matrix for an outer product - Matrix mat_n(normal.storage(), normal.size(), 1.); - - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof - // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 - - // construct A matrix - const ElementType & type = element.master.type; - auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - A(j, i * spatial_dimension + j) = 1; - continue; - } - A(j, i * spatial_dimension + j) = -shapes[i - 1]; - } - } - - // computing shape derivatives - auto && shape_derivatives = - ElementClassHelper<_ek_regular>::getDNDS(projection, type); - - Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - auto construct_Aj = [&](auto && dnds) { - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - Aj(j, i * spatial_dimension + j) = 0; - continue; - } - Aj(j, i * spatial_dimension + j) = dnds(i - 1); - } - } - }; - - // tangents should have been calculated in normal modulii - auto & tangents = model.getTangents(); - Matrix covariant_basis( - tangents.begin(surface_dimension, spatial_dimension)[element.slave]); - - auto & tangential_tractions = model.getTangentialTractions(); - Vector tangential_traction( - tangential_tractions.begin(surface_dimension)[element.slave]); - - // compute norm of trial traction - Real traction_norm = 0; - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - for (auto i : arange(surface_dimension)) { - for (auto j : arange(surface_dimension)) { - traction_norm += tangential_traction[i] * tangential_traction[j] * - contravariant_metric_tensor(i, j); - } - } - traction_norm = sqrt(traction_norm); - - // construct four parts of stick modulii (eq 107,107a-c) - Matrix k_first(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix k_second(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix k_third(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix k_fourth(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - for (auto && values1 : enumerate(covariant_basis.transpose())) { - auto & alpha = std::get<0>(values1); - auto & tangent_alpha = std::get<1>(values1); - - Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); - - Matrix t_outer_n(spatial_dimension, spatial_dimension); - Matrix t_outer_t(spatial_dimension, spatial_dimension); - - for (auto && values2 : - zip(arange(surface_dimension), covariant_basis.transpose(), - shape_derivatives.transpose())) { - auto & beta = std::get<0>(values2); - auto & tangent_beta = std::get<1>(values2); - auto & dnds = std::get<2>(values2); - // construct Aj from shape function wrt to jth natural - // coordinate - construct_Aj(dnds); - - // eq 107 - Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); - t_outer_n.mul(mat_t_beta, mat_n); - - Matrix tmp(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp1(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp.mul(t_outer_n, A); - tmp1.mul(A, tmp); - - tmp1 *= epsilon_n * mu * tangential_traction[alpha] * - contravariant_metric_tensor(alpha, beta); - tmp1 /= traction_norm; - - k_first += tmp1 * nodal_area; - - // eq 107a - t_outer_t.mul(mat_t_alpha, mat_t_beta); - - tmp.mul(t_outer_t, A); - tmp1.mul(A, tmp); - - tmp1 *= epsilon_t * mu * p_n * contravariant_metric_tensor(alpha, beta); - tmp1 /= traction_norm; - - k_second += tmp1 * nodal_area; - - for (auto && values3 : enumerate(covariant_basis.transpose())) { - auto & gamma = std::get<0>(values3); - auto & tangent_gamma = std::get<1>(values3); - - Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), - 1.); - - for (auto && values4 : enumerate(covariant_basis.transpose())) { - auto & theta = std::get<0>(values4); - auto & tangent_theta = std::get<1>(values4); - - Matrix mat_t_theta(tangent_theta.storage(), - tangent_theta.size(), 1.); - t_outer_t.mul(mat_t_gamma, mat_t_theta); - - // eq 107b - tmp.mul(t_outer_t, A); - tmp1.mul(A, tmp); - - tmp1 *= epsilon_t * mu * p_n * tangential_traction[alpha] * - tangential_traction[beta]; - tmp1 *= contravariant_metric_tensor(alpha, gamma) * - contravariant_metric_tensor(beta, theta); - tmp1 /= pow(traction_norm, 3); - - k_third += tmp1 * nodal_area; - - // eq 107c - tmp.mul(t_outer_t, Aj); - tmp1.mul(A, tmp); - tmp1 *= contravariant_metric_tensor(alpha, theta) * - contravariant_metric_tensor(beta, gamma); - tmp1 *= mu * p_n * tangential_traction[alpha]; - tmp1 /= traction_norm; - - Matrix tmp2(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp3(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp2.mul(t_outer_t, A); - tmp3.mul(Aj, tmp2); - tmp3 *= contravariant_metric_tensor(alpha, gamma) * - contravariant_metric_tensor(beta, theta); - tmp3 *= mu * p_n * tangential_traction[alpha]; - tmp3 /= traction_norm; - - k_fourth += (tmp1 + tmp3) * nodal_area; - } - } - } - } - - stiffness += k_third + k_fourth - k_first - k_second; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::beforeSolveStep() {} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::afterSolveStep(__attribute__((unused)) bool converged) { -} - -INSTANTIATE_RESOLUTION(penalty_linear, ResolutionPenalty); - } // namespace akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.hh b/src/model/contact_mechanics/resolutions/resolution_penalty.hh index f67d11ee6..47697dfa6 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.hh +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.hh @@ -1,142 +1,135 @@ /** * @file resolution_penalty.hh * * @author Mohit Pundir * * @date creation: Fri Jun 18 2010 - * @date last modification: Wed Jun 09 2021 + * @date last modification: Fri Oct 21 2022 * * @brief Linear Penalty Resolution for Contact Mechanics Model * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "resolution.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RESOLUTION_PENALTY_HH__ #define __AKANTU_RESOLUTION_PENALTY_HH__ namespace akantu { -enum PenaltyType { - penalty_linear, - penalty_quadratic -}; -class ResolutionPenalty : public Resolution { +template +class ResolutionPenaltyTemplate : public Resolution { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - ResolutionPenalty(ContactMechanicsModel & model, const ID & id = "", - PenaltyType penalty_type = penalty_linear); + ResolutionPenaltyTemplate(ContactMechanicsModel & model, const ID & id = ""); - ~ResolutionPenalty() override = default; + ~ResolutionPenaltyTemplate() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize the resolution void initialize(); /* ------------------------------------------------------------------------ */ /* Methods for stiffness computation */ /* ------------------------------------------------------------------------ */ protected: /// local computaion of stiffness matrix due to stick state void computeStickModuli(const ContactElement & /*element*/, Matrix & /*stiffness*/); /// local computation of stiffness matrix due to slip state void computeSlipModuli(const ContactElement & /*element*/, Matrix & /*stiffness*/); public: /// local computation of tangent moduli due to normal traction void computeNormalModuli(const ContactElement & /*element*/, Matrix & /*stiffness*/) override; /// local computation of tangent moduli due to tangential traction void computeTangentialModuli(const ContactElement & /*element*/, Matrix & /*stiffness*/) override; /* ------------------------------------------------------------------------ */ /* Methods for force computation */ /* ------------------------------------------------------------------------ */ public: /// local computation of normal force due to normal contact void computeNormalForce(const ContactElement & /*element*/, Vector & /*force*/) override; /// local computation of tangential force due to frictional traction void computeTangentialForce(const ContactElement & /*element*/, Vector & /*force*/) override; protected: /// local computation of normal traction due to penetration Real computeNormalTraction(Real & /*gap*/) const; /// local computation of trial tangential traction due to friction void computeTrialTangentialTraction(const ContactElement & /*element*/, const Matrix & /*covariant_basis*/, Vector & /*traction*/); /// local computation of tangential traction due to stick void computeStickTangentialTraction(const ContactElement & /*unused*/, Vector & /*traction_trial*/, Vector & /*traction_tangential*/); /// local computation of tangential traction due to slip void computeSlipTangentialTraction(const ContactElement & /*element*/, const Matrix & /*covariant_basis*/, Vector & /*traction_trial*/, Vector & /*traction_tangential*/); /// local computation of tangential traction due to friction void computeTangentialTraction(const ContactElement & /*element*/, const Matrix & /*covariant_basis*/, Vector & /*traction_tangential*/); public: void beforeSolveStep() override; void afterSolveStep(bool converged = true) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// penalty parameter for normal traction Real epsilon_n; /// penalty parameter for tangential traction Real epsilon_t; - - /// Penalty type to be used for penetration function - PenaltyType penalty_type; }; } // namespace akantu #endif /* __AKANTU_RESOLUTION_PENALTY_HH__ */ diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc index e7e9b2eeb..27856dcbd 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc @@ -1,51 +1,40 @@ /** * @file resolution_penalty_quadratic.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 - * @date last modification: Wed Jun 09 2021 + * @date last modification: Fri Oct 21 2022 * * @brief Specialization of the resolution class for the quadratic penalty * method * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "resolution_penalty_quadratic.hh" #include "element_class_helper.hh" /* -------------------------------------------------------------------------- */ namespace akantu { -/* -------------------------------------------------------------------------- */ -ResolutionPenaltyQuadratic::ResolutionPenaltyQuadratic( - ContactMechanicsModel & model, const ID & id) - : ResolutionPenalty(model, id, penalty_quadratic) { - AKANTU_DEBUG_IN(); - this->initialize(); - AKANTU_DEBUG_OUT(); -} - -INSTANTIATE_RESOLUTION(penalty_quadratic, ResolutionPenaltyQuadratic); - } // namespace akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh index 4f0930dc0..0e91a6f5a 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh +++ b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh @@ -1,55 +1,57 @@ /** * @file resolution_penalty_quadratic.hh * * @author Mohit Pundir * * @date creation: Fri Jun 18 2010 - * @date last modification: Mon Aug 10 2020 + * @date last modification: Fri Oct 21 2022 * * @brief Quadratic Penalty Resolution for Contact Mechanics Model * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "resolution_penalty.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RESOLUTION_PENALTY_QUADRATIC_HH__ #define __AKANTU_RESOLUTION_PENALTY_QUADRATIC_HH__ namespace akantu { -class ResolutionPenaltyQuadratic : public ResolutionPenalty { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: - ResolutionPenaltyQuadratic(ContactMechanicsModel & model, const ID & id = ""); +// Define function and derivative for quadratic penalty method +template T quadraticPenetrationFunction(T var) { + return var*var + var; +} - ~ResolutionPenaltyQuadratic() override = default; - -}; +template T quadraticPenetrationDerivative(T var) { + return 2.0*var + 1.0; +} + +// Instantiate quadratic penalty as a resolution +INSTANTIATE_RESOLUTION(penalty_quadratic, ResolutionPenaltyTemplate, Real, + quadraticPenetrationFunction, quadraticPenetrationDerivative); } // namespace akantu #endif /* __AKANTU_RESOLUTION_PENALTY_QUADRATIC_HH__ */ diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty_tmpl.hh similarity index 83% copy from src/model/contact_mechanics/resolutions/resolution_penalty.cc copy to src/model/contact_mechanics/resolutions/resolution_penalty_tmpl.hh index 35ad7330b..6ac8c1617 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty_tmpl.hh @@ -1,868 +1,872 @@ -/** - * @file resolution_penalty.cc - * - * @author Mohit Pundir - * - * @date creation: Thu Jan 17 2019 - * @date last modification: Wed Jun 09 2021 - * - * @brief Specialization of the resolution class for the penalty method - * - * - * @section LICENSE - * - * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * Akantu is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Akantu. If not, see . - * - */ - -/* -------------------------------------------------------------------------- */ -#include "resolution_penalty.hh" -#include "element_class_helper.hh" -/* -------------------------------------------------------------------------- */ - -namespace akantu { - -/* -------------------------------------------------------------------------- */ -ResolutionPenalty::ResolutionPenalty(ContactMechanicsModel & model, - const ID & id, PenaltyType penalty_type) - : Resolution(model, id) { - AKANTU_DEBUG_IN(); - this->initialize(); - AKANTU_DEBUG_OUT(); - this->penalty_type = penalty_type; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::initialize() { - this->registerParam("epsilon_n", epsilon_n, Real(0.), - _pat_parsable | _pat_modifiable, - "Normal penalty parameter"); - this->registerParam("epsilon_t", epsilon_t, Real(0.), - _pat_parsable | _pat_modifiable, - "Tangential penalty parameter"); -} - -/* -------------------------------------------------------------------------- */ -Real ResolutionPenalty::computeNormalTraction(Real & gap) const { - switch(this->penalty_type) { - case penalty_linear: - return epsilon_n * macaulay(gap); - case penalty_quadratic: - return epsilon_n * (macaulay(gap) * macaulay(gap) + macaulay(gap)); - default: AKANTU_ERROR("This case should not happen!!"); - } -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeNormalForce(const ContactElement & element, - Vector & force) { - - force.zero(); - - auto & gaps = model.getGaps(); - auto & projections = model.getProjections(); - auto & normals = model.getNormals(); - - auto surface_dimension = spatial_dimension - 1; - - Real gap(gaps.begin()[element.slave]); - Vector normal(normals.begin(spatial_dimension)[element.slave]); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & nodal_area = const_cast &>(model.getNodalArea()); - - // compute normal traction - Real p_n = computeNormalTraction(gap); - p_n *= nodal_area[element.slave]; - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix shape_matric(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - ResolutionUtils::computeShapeFunctionMatric(element, projection, - shape_matric); - - force.mul(shape_matric, normal, p_n); -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentialForce(const ContactElement & element, - Vector & force) { - - if (mu == 0) { - return; - } - - force.zero(); - - UInt surface_dimension = spatial_dimension - 1; - - // compute covariant basis - auto & projections = model.getProjections(); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & normals = model.getNormals(); - Vector normal(normals.begin(spatial_dimension)[element.slave]); - - auto & tangents = model.getTangents(); - Matrix covariant_basis( - tangents.begin(surface_dimension, spatial_dimension)[element.slave]); - - // check for no-contact to contact condition - // need a better way to check if new node added is not presnt in the - // previous master elemets - auto & previous_master_elements = model.getPreviousMasterElements(); - if (element.slave >= previous_master_elements.size()) { - return; - } - - auto & previous_element = previous_master_elements[element.slave]; - if (previous_element.type == _not_defined) { - return; - } - - // compute tangential traction using return map algorithm - auto & tangential_tractions = model.getTangentialTractions(); - Vector tangential_traction( - tangential_tractions.begin(surface_dimension)[element.slave]); - this->computeTangentialTraction(element, covariant_basis, - tangential_traction); - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix shape_matric(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - ResolutionUtils::computeShapeFunctionMatric(element, projection, - shape_matric); - - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - auto & nodal_area = const_cast &>(model.getNodalArea()); - - for (auto && values1 : enumerate(covariant_basis.transpose())) { - auto & alpha = std::get<0>(values1); - auto & tangent_alpha = std::get<1>(values1); - for (auto && values2 : enumerate(tangential_traction)) { - auto & beta = std::get<0>(values2); - auto & traction_beta = std::get<1>(values2); - Vector tmp(force.size()); - tmp.mul(shape_matric, tangent_alpha, traction_beta); - tmp *= - contravariant_metric_tensor(alpha, beta) * nodal_area[element.slave]; - force += tmp; - } - } -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentialTraction( - const ContactElement & element, const Matrix & covariant_basis, - Vector & traction_tangential) { - - UInt surface_dimension = spatial_dimension - 1; - - auto & gaps = model.getGaps(); - auto & gap = gaps.begin()[element.slave]; - - // Return map algorithm is employed - // compute trial traction - Vector traction_trial(surface_dimension); - this->computeTrialTangentialTraction(element, covariant_basis, - traction_trial); - - // compute norm of trial traction - Real traction_trial_norm = 0; - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - for (auto i : arange(surface_dimension)) { - for (auto j : arange(surface_dimension)) { - traction_trial_norm += traction_trial[i] * traction_trial[j] * - contravariant_metric_tensor(i, j); - } - } - traction_trial_norm = sqrt(traction_trial_norm); - - // check stick or slip condition - auto & contact_state = model.getContactState(); - auto & state = contact_state.begin()[element.slave]; - - Real p_n = computeNormalTraction(gap); - bool stick = (traction_trial_norm <= mu * p_n); - - if (stick) { - state = ContactState::_stick; - computeStickTangentialTraction(element, traction_trial, - traction_tangential); - } else { - state = ContactState::_slip; - computeSlipTangentialTraction(element, covariant_basis, traction_trial, - traction_tangential); - } -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTrialTangentialTraction( - const ContactElement & element, const Matrix & covariant_basis, - Vector & traction) { - - UInt surface_dimension = spatial_dimension - 1; - - auto & projections = model.getProjections(); - Vector current_projection( - projections.begin(surface_dimension)[element.slave]); - - auto & previous_projections = model.getPreviousProjections(); - Vector previous_projection( - previous_projections.begin(surface_dimension)[element.slave]); - - // method from Laursen et. al. - /*auto covariant_metric_tensor = - GeometryUtils::covariantMetricTensor(covariant_basis); auto - increment_projection = current_projection - previous_projection; - - traction.mul(covariant_metric_tensor, increment_projection, epsilon_t); - - auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); - Vector - previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]); - traction = previous_traction + traction;*/ - - // method from Schweizerhof - auto covariant_metric_tensor = - GeometryUtils::covariantMetricTensor(covariant_basis); - - auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); - Vector previous_traction( - previous_tangential_tractions.begin(surface_dimension)[element.slave]); - - auto & previous_tangents = model.getPreviousTangents(); - Matrix previous_covariant_basis(previous_tangents.begin( - surface_dimension, spatial_dimension)[element.slave]); - auto previous_contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(previous_covariant_basis); - - auto current_tangent = covariant_basis.transpose(); - auto previous_tangent = previous_covariant_basis.transpose(); - - for (auto alpha : arange(surface_dimension)) { - Vector tangent_alpha(current_tangent(alpha)); - for (auto gamma : arange(surface_dimension)) { - for (auto beta : arange(surface_dimension)) { - Vector tangent_beta(previous_tangent(beta)); - auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha); - traction[alpha] += previous_traction[gamma] * - previous_contravariant_metric_tensor(gamma, beta) * - t_alpha_t_beta; - } - } - } - - auto & previous_master_elements = model.getPreviousMasterElements(); - auto & previous_element = previous_master_elements[element.slave]; - - Vector previous_real_projection(spatial_dimension); - GeometryUtils::realProjection( - model.getMesh(), model.getContactDetector().getPositions(), - previous_element, previous_projection, previous_real_projection); - - Vector current_real_projection(spatial_dimension); - GeometryUtils::realProjection( - model.getMesh(), model.getContactDetector().getPositions(), - element.master, current_projection, current_real_projection); - - auto increment_real = current_real_projection - previous_real_projection; - Vector increment_xi(surface_dimension); - - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - // increment in natural coordinate - for (auto beta : arange(surface_dimension)) { - for (auto gamma : arange(surface_dimension)) { - auto temp = increment_real.dot(current_tangent(gamma)); - temp *= contravariant_metric_tensor(beta, gamma); - increment_xi[beta] += temp; - } - } - - Vector temp(surface_dimension); - temp.mul(covariant_metric_tensor, increment_xi, epsilon_t); - - traction -= temp; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeStickTangentialTraction( - const ContactElement & /*element*/, Vector & traction_trial, - Vector & traction_tangential) { - traction_tangential = traction_trial; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeSlipTangentialTraction( - const ContactElement & element, const Matrix & covariant_basis, - Vector & traction_trial, Vector & traction_tangential) { - UInt surface_dimension = spatial_dimension - 1; - - auto & gaps = model.getGaps(); - auto & gap = gaps.begin()[element.slave]; - - // compute norm of trial traction - Real traction_trial_norm = 0; - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - for (auto alpha : arange(surface_dimension)) { - for (auto beta : arange(surface_dimension)) { - traction_trial_norm += traction_trial[alpha] * traction_trial[beta] * - contravariant_metric_tensor(alpha, beta); - } - } - traction_trial_norm = sqrt(traction_trial_norm); - - auto slip_direction = traction_trial; - slip_direction /= traction_trial_norm; - - Real p_n = computeNormalTraction(gap); - traction_tangential = slip_direction; - traction_tangential *= mu * p_n; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeNormalModuli(const ContactElement & element, - Matrix & stiffness) { - - auto surface_dimension = spatial_dimension - 1; - - auto & gaps = model.getGaps(); - Real gap(gaps.begin()[element.slave]); - - auto & projections = model.getProjections(); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & nodal_areas = model.getNodalArea(); - auto & nodal_area = nodal_areas.begin()[element.slave]; - - auto & normals = model.getNormals(); - Vector normal(normals.begin(spatial_dimension)[element.slave]); - - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof - // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 - - // construct A matrix - const ElementType & type = element.master.type; - auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - A(j, i * spatial_dimension + j) = 1; - continue; - } - A(j, i * spatial_dimension + j) = -shapes[i - 1]; - } - } - - // construct the main part of normal matrix - Matrix k_main(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - Matrix n_outer_n(spatial_dimension, spatial_dimension); - Matrix mat_n(normal.storage(), normal.size(), 1.); - n_outer_n.mul(mat_n, mat_n); - - Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - tmp.mul(n_outer_n, A); - - k_main.mul(A, tmp); - switch(this->penalty_type) { - case penalty_linear: - k_main *= epsilon_n * heaviside(gap) * nodal_area; - break; - case penalty_quadratic: - k_main *= epsilon_n * heaviside(gap) * (2 * gap + 1) * nodal_area; - break; - default: AKANTU_ERROR("This case should not happen!!"); - } - - // construct the rotational part of the normal matrix - auto & tangents = model.getTangents(); - Matrix covariant_basis( - tangents.begin(surface_dimension, spatial_dimension)[element.slave]); - - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - // computing shape derivatives - auto && shape_derivatives = - ElementClassHelper<_ek_regular>::getDNDS(projection, type); - - // consists of 2 rotational parts - Matrix k_rot1(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix k_rot2(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - auto construct_Aj = [&](auto && dnds) { - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - Aj(j, i * spatial_dimension + j) = 0; - continue; - } - Aj(j, i * spatial_dimension + j) = dnds(i - 1); - } - } - }; - - for (auto && values1 : enumerate(covariant_basis.transpose())) { - auto & alpha = std::get<0>(values1); - auto & tangent = std::get<1>(values1); - - Matrix n_outer_t(spatial_dimension, spatial_dimension); - Matrix mat_t(tangent.storage(), tangent.size(), 1.); - n_outer_t.mul(mat_n, mat_t); - - Matrix t_outer_n(spatial_dimension, spatial_dimension); - t_outer_n.mul(mat_t, mat_n); - - for (auto && values2 : enumerate(shape_derivatives.transpose())) { - auto & beta = std::get<0>(values2); - auto & dnds = std::get<1>(values2); - // construct Aj from shape function wrt to jth natural - // coordinate - construct_Aj(dnds); - - Matrix tmp(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp1(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp.mul(n_outer_t, A); - tmp1.mul(Aj, tmp); - tmp1 *= contravariant_metric_tensor(alpha, beta); - k_rot1 += tmp1; - - tmp.mul(t_outer_n, Aj); - tmp1.mul(A, tmp); - tmp1 *= contravariant_metric_tensor(alpha, beta); - k_rot2 += tmp1; - } - } - - switch(this->penalty_type) { - case penalty_linear: - k_rot1 *= -epsilon_n * heaviside(gap) * gap * nodal_area; - k_rot2 *= -epsilon_n * heaviside(gap) * gap * nodal_area; - break; - case penalty_quadratic: - k_rot1 *= -epsilon_n * heaviside(gap) * (gap * gap + gap) * nodal_area; - k_rot2 *= -epsilon_n * heaviside(gap) * (gap * gap + gap) * nodal_area; - break; - default: AKANTU_ERROR("This case should not happen!!"); - } - - stiffness += k_main + k_rot1 + k_rot2; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentialModuli(const ContactElement & element, - Matrix & stiffness) { - if (mu == 0) { - return; - } - - stiffness.zero(); - - auto & contact_state = model.getContactState(); - auto state = contact_state.begin()[element.slave]; - - switch (state) { - case ContactState::_stick: { - computeStickModuli(element, stiffness); - break; - } - case ContactState::_slip: { - computeSlipModuli(element, stiffness); - break; - } - default: - break; - } -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeStickModuli(const ContactElement & element, - Matrix & stiffness) { - - auto surface_dimension = spatial_dimension - 1; - - auto & projections = model.getProjections(); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & nodal_areas = model.getNodalArea(); - auto & nodal_area = nodal_areas.begin()[element.slave]; - - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof - // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 - - // construct A matrix - const ElementType & type = element.master.type; - auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - A(j, i * spatial_dimension + j) = 1; - continue; - } - A(j, i * spatial_dimension + j) = -shapes[i - 1]; - } - } - - // computing shape derivatives - auto && shape_derivatives = - ElementClassHelper<_ek_regular>::getDNDS(projection, type); - - Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - auto construct_Aj = [&](auto && dnds) { - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - Aj(j, i * spatial_dimension + j) = 0; - continue; - } - Aj(j, i * spatial_dimension + j) = dnds(i - 1); - } - } - }; - - // tangents should have been calculated in normal modulii - auto & tangents = model.getTangents(); - Matrix covariant_basis( - tangents.begin(surface_dimension, spatial_dimension)[element.slave]); - - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - // construct 1st part of the stick modulii - Matrix k_main(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - for (auto && values1 : enumerate(covariant_basis.transpose())) { - auto & alpha = std::get<0>(values1); - auto & tangent_alpha = std::get<1>(values1); - - Matrix t_outer_t(spatial_dimension, spatial_dimension); - Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); - - for (auto && values2 : enumerate(covariant_basis.transpose())) { - auto & beta = std::get<0>(values2); - auto & tangent_beta = std::get<1>(values2); - - Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); - t_outer_t.mul(mat_t_alpha, mat_t_beta); - - Matrix tmp(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp1(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp.mul(t_outer_t, A); - tmp1.mul(A, tmp); - tmp1 *= contravariant_metric_tensor(alpha, beta); - k_main += tmp1; - } - } - - k_main *= -epsilon_t; - - // construct 2nd part of the stick modulii - auto & tangential_tractions = model.getTangentialTractions(); - Vector tangential_traction( - tangential_tractions.begin(surface_dimension)[element.slave]); - - Matrix k_second(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - for (auto alpha : arange(surface_dimension)) { - - Matrix k_sum(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - for (auto && values1 : enumerate(shape_derivatives.transpose())) { - auto & beta = std::get<0>(values1); - auto & dnds = std::get<1>(values1); - // construct Aj from shape function wrt to jth natural - // coordinate - construct_Aj(dnds); - for (auto && values2 : enumerate(covariant_basis.transpose())) { - auto & gamma = std::get<0>(values2); - auto & tangent_gamma = std::get<1>(values2); - - Matrix t_outer_t(spatial_dimension, spatial_dimension); - Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), - 1.); - - for (auto && values3 : enumerate(covariant_basis.transpose())) { - auto & theta = std::get<0>(values3); - auto & tangent_theta = std::get<1>(values3); - - Matrix mat_t_theta(tangent_theta.storage(), - tangent_theta.size(), 1.); - t_outer_t.mul(mat_t_gamma, mat_t_theta); - - Matrix tmp(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp1(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp.mul(t_outer_t, Aj); - tmp1.mul(A, tmp); - tmp1 *= contravariant_metric_tensor(alpha, theta) * - contravariant_metric_tensor(beta, gamma); - - Matrix tmp2(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp3(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp2.mul(t_outer_t, A); - tmp3.mul(Aj, tmp2); - tmp3 *= contravariant_metric_tensor(alpha, gamma) * - contravariant_metric_tensor(beta, theta); - - k_sum += tmp1 + tmp3; - } - } - } - - k_second += tangential_traction[alpha] * k_sum; - } - - stiffness += k_main * nodal_area - k_second * nodal_area; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeSlipModuli(const ContactElement & element, - Matrix & stiffness) { - - auto surface_dimension = spatial_dimension - 1; - - auto & gaps = model.getGaps(); - Real gap(gaps.begin()[element.slave]); - - auto & nodal_areas = model.getNodalArea(); - auto & nodal_area = nodal_areas.begin()[element.slave]; - - // compute normal traction - Real p_n = computeNormalTraction(gap); - - auto & projections = model.getProjections(); - Vector projection(projections.begin(surface_dimension)[element.slave]); - - auto & normals = model.getNormals(); - Vector normal(normals.begin(spatial_dimension)[element.slave]); - - // restructure normal as a matrix for an outer product - Matrix mat_n(normal.storage(), normal.size(), 1.); - - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof - // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 - - // construct A matrix - const ElementType & type = element.master.type; - auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); - - UInt nb_nodes_per_contact = element.getNbNodes(); - Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - A(j, i * spatial_dimension + j) = 1; - continue; - } - A(j, i * spatial_dimension + j) = -shapes[i - 1]; - } - } - - // computing shape derivatives - auto && shape_derivatives = - ElementClassHelper<_ek_regular>::getDNDS(projection, type); - - Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); - - auto construct_Aj = [&](auto && dnds) { - for (auto i : arange(nb_nodes_per_contact)) { - for (auto j : arange(spatial_dimension)) { - if (i == 0) { - Aj(j, i * spatial_dimension + j) = 0; - continue; - } - Aj(j, i * spatial_dimension + j) = dnds(i - 1); - } - } - }; - - // tangents should have been calculated in normal modulii - auto & tangents = model.getTangents(); - Matrix covariant_basis( - tangents.begin(surface_dimension, spatial_dimension)[element.slave]); - - auto & tangential_tractions = model.getTangentialTractions(); - Vector tangential_traction( - tangential_tractions.begin(surface_dimension)[element.slave]); - - // compute norm of trial traction - Real traction_norm = 0; - auto contravariant_metric_tensor = - GeometryUtils::contravariantMetricTensor(covariant_basis); - - for (auto i : arange(surface_dimension)) { - for (auto j : arange(surface_dimension)) { - traction_norm += tangential_traction[i] * tangential_traction[j] * - contravariant_metric_tensor(i, j); - } - } - traction_norm = sqrt(traction_norm); - - // construct four parts of stick modulii (eq 107,107a-c) - Matrix k_first(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix k_second(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix k_third(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - Matrix k_fourth(nb_nodes_per_contact * spatial_dimension, - nb_nodes_per_contact * spatial_dimension); - - for (auto && values1 : enumerate(covariant_basis.transpose())) { - auto & alpha = std::get<0>(values1); - auto & tangent_alpha = std::get<1>(values1); - - Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); - - Matrix t_outer_n(spatial_dimension, spatial_dimension); - Matrix t_outer_t(spatial_dimension, spatial_dimension); - - for (auto && values2 : - zip(arange(surface_dimension), covariant_basis.transpose(), - shape_derivatives.transpose())) { - auto & beta = std::get<0>(values2); - auto & tangent_beta = std::get<1>(values2); - auto & dnds = std::get<2>(values2); - // construct Aj from shape function wrt to jth natural - // coordinate - construct_Aj(dnds); - - // eq 107 - Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); - t_outer_n.mul(mat_t_beta, mat_n); - - Matrix tmp(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp1(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp.mul(t_outer_n, A); - tmp1.mul(A, tmp); - - tmp1 *= epsilon_n * mu * tangential_traction[alpha] * - contravariant_metric_tensor(alpha, beta); - tmp1 /= traction_norm; - - k_first += tmp1 * nodal_area; - - // eq 107a - t_outer_t.mul(mat_t_alpha, mat_t_beta); - - tmp.mul(t_outer_t, A); - tmp1.mul(A, tmp); - - tmp1 *= epsilon_t * mu * p_n * contravariant_metric_tensor(alpha, beta); - tmp1 /= traction_norm; - - k_second += tmp1 * nodal_area; - - for (auto && values3 : enumerate(covariant_basis.transpose())) { - auto & gamma = std::get<0>(values3); - auto & tangent_gamma = std::get<1>(values3); - - Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), - 1.); - - for (auto && values4 : enumerate(covariant_basis.transpose())) { - auto & theta = std::get<0>(values4); - auto & tangent_theta = std::get<1>(values4); - - Matrix mat_t_theta(tangent_theta.storage(), - tangent_theta.size(), 1.); - t_outer_t.mul(mat_t_gamma, mat_t_theta); - - // eq 107b - tmp.mul(t_outer_t, A); - tmp1.mul(A, tmp); - - tmp1 *= epsilon_t * mu * p_n * tangential_traction[alpha] * - tangential_traction[beta]; - tmp1 *= contravariant_metric_tensor(alpha, gamma) * - contravariant_metric_tensor(beta, theta); - tmp1 /= pow(traction_norm, 3); - - k_third += tmp1 * nodal_area; - - // eq 107c - tmp.mul(t_outer_t, Aj); - tmp1.mul(A, tmp); - tmp1 *= contravariant_metric_tensor(alpha, theta) * - contravariant_metric_tensor(beta, gamma); - tmp1 *= mu * p_n * tangential_traction[alpha]; - tmp1 /= traction_norm; - - Matrix tmp2(spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - Matrix tmp3(nb_nodes_per_contact * spatial_dimension, - spatial_dimension * nb_nodes_per_contact); - tmp2.mul(t_outer_t, A); - tmp3.mul(Aj, tmp2); - tmp3 *= contravariant_metric_tensor(alpha, gamma) * - contravariant_metric_tensor(beta, theta); - tmp3 *= mu * p_n * tangential_traction[alpha]; - tmp3 /= traction_norm; - - k_fourth += (tmp1 + tmp3) * nodal_area; - } - } - } - } - - stiffness += k_third + k_fourth - k_first - k_second; -} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::beforeSolveStep() {} - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::afterSolveStep(__attribute__((unused)) bool converged) { -} - -INSTANTIATE_RESOLUTION(penalty_linear, ResolutionPenalty); - -} // namespace akantu +/** +* @file resolution_penalty_tmpl.hh +* +* @author Fabio Matti +* +* @date creation: Fri Oct 21 2022 +* @date last modification: Fri Oct 21 2022 +* +* @brief Implementation of the templated penalty method +* +* +* @section LICENSE +* +* Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) +* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +* +* Akantu is free software: you can redistribute it and/or modify it under the +* terms of the GNU Lesser General Public License as published by the Free +* Software Foundation, either version 3 of the License, or (at your option) any +* later version. +* +* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +* details. +* +* You should have received a copy of the GNU Lesser General Public License +* along with Akantu. If not, see . +* +*/ + +#include "resolution_penalty.hh" +#include "element_class_helper.hh" + +#ifndef AKANTU_RESOLUTION_PENALTY_TMPL_HH +#define AKANTU_RESOLUTION_PENALTY_TMPL_HH + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +template +ResolutionPenaltyTemplate::ResolutionPenaltyTemplate(ContactMechanicsModel & model, + const ID & id) + : Resolution(model, id) { + AKANTU_DEBUG_IN(); + this->initialize(); + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::initialize() { + this->registerParam("epsilon_n", epsilon_n, Real(0.), + _pat_parsable | _pat_modifiable, + "Normal penalty parameter"); + this->registerParam("epsilon_t", epsilon_t, Real(0.), + _pat_parsable | _pat_modifiable, + "Tangential penalty parameter"); +} + +/* -------------------------------------------------------------------------- */ +template +Real ResolutionPenaltyTemplate::computeNormalTraction(Real & gap) const { + return epsilon_n * penetrationFunction(macaulay(gap)); +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeNormalForce(const ContactElement & element, + Vector & force) { + + force.zero(); + + auto & gaps = model.getGaps(); + auto & projections = model.getProjections(); + auto & normals = model.getNormals(); + + auto surface_dimension = spatial_dimension - 1; + + Real gap(gaps.begin()[element.slave]); + Vector normal(normals.begin(spatial_dimension)[element.slave]); + Vector projection(projections.begin(surface_dimension)[element.slave]); + + auto & nodal_area = const_cast &>(model.getNodalArea()); + + // compute normal traction + Real p_n = computeNormalTraction(gap); + p_n *= nodal_area[element.slave]; + + UInt nb_nodes_per_contact = element.getNbNodes(); + Matrix shape_matric(spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + ResolutionUtils::computeShapeFunctionMatric(element, projection, + shape_matric); + + force.mul(shape_matric, normal, p_n); +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeTangentialForce(const ContactElement & element, + Vector & force) { + + if (mu == 0) { + return; + } + + force.zero(); + + UInt surface_dimension = spatial_dimension - 1; + + // compute covariant basis + auto & projections = model.getProjections(); + Vector projection(projections.begin(surface_dimension)[element.slave]); + + auto & normals = model.getNormals(); + Vector normal(normals.begin(spatial_dimension)[element.slave]); + + auto & tangents = model.getTangents(); + Matrix covariant_basis( + tangents.begin(surface_dimension, spatial_dimension)[element.slave]); + + // check for no-contact to contact condition + // need a better way to check if new node added is not presnt in the + // previous master elemets + auto & previous_master_elements = model.getPreviousMasterElements(); + if (element.slave >= previous_master_elements.size()) { + return; + } + + auto & previous_element = previous_master_elements[element.slave]; + if (previous_element.type == _not_defined) { + return; + } + + // compute tangential traction using return map algorithm + auto & tangential_tractions = model.getTangentialTractions(); + Vector tangential_traction( + tangential_tractions.begin(surface_dimension)[element.slave]); + this->computeTangentialTraction(element, covariant_basis, + tangential_traction); + + UInt nb_nodes_per_contact = element.getNbNodes(); + Matrix shape_matric(spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + ResolutionUtils::computeShapeFunctionMatric(element, projection, + shape_matric); + + auto contravariant_metric_tensor = + GeometryUtils::contravariantMetricTensor(covariant_basis); + + auto & nodal_area = const_cast &>(model.getNodalArea()); + + for (auto && values1 : enumerate(covariant_basis.transpose())) { + auto & alpha = std::get<0>(values1); + auto & tangent_alpha = std::get<1>(values1); + for (auto && values2 : enumerate(tangential_traction)) { + auto & beta = std::get<0>(values2); + auto & traction_beta = std::get<1>(values2); + Vector tmp(force.size()); + tmp.mul(shape_matric, tangent_alpha, traction_beta); + tmp *= + contravariant_metric_tensor(alpha, beta) * nodal_area[element.slave]; + force += tmp; + } + } +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeTangentialTraction( + const ContactElement & element, const Matrix & covariant_basis, + Vector & traction_tangential) { + + UInt surface_dimension = spatial_dimension - 1; + + auto & gaps = model.getGaps(); + auto & gap = gaps.begin()[element.slave]; + + // Return map algorithm is employed + // compute trial traction + Vector traction_trial(surface_dimension); + this->computeTrialTangentialTraction(element, covariant_basis, + traction_trial); + + // compute norm of trial traction + Real traction_trial_norm = 0; + auto contravariant_metric_tensor = + GeometryUtils::contravariantMetricTensor(covariant_basis); + for (auto i : arange(surface_dimension)) { + for (auto j : arange(surface_dimension)) { + traction_trial_norm += traction_trial[i] * traction_trial[j] * + contravariant_metric_tensor(i, j); + } + } + traction_trial_norm = sqrt(traction_trial_norm); + + // check stick or slip condition + auto & contact_state = model.getContactState(); + auto & state = contact_state.begin()[element.slave]; + + Real p_n = computeNormalTraction(gap); + bool stick = (traction_trial_norm <= mu * p_n); + + if (stick) { + state = ContactState::_stick; + computeStickTangentialTraction(element, traction_trial, + traction_tangential); + } else { + state = ContactState::_slip; + computeSlipTangentialTraction(element, covariant_basis, traction_trial, + traction_tangential); + } +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeTrialTangentialTraction( + const ContactElement & element, const Matrix & covariant_basis, + Vector & traction) { + + UInt surface_dimension = spatial_dimension - 1; + + auto & projections = model.getProjections(); + Vector current_projection( + projections.begin(surface_dimension)[element.slave]); + + auto & previous_projections = model.getPreviousProjections(); + Vector previous_projection( + previous_projections.begin(surface_dimension)[element.slave]); + + // method from Laursen et. al. + /*auto covariant_metric_tensor = + GeometryUtils::covariantMetricTensor(covariant_basis); auto + increment_projection = current_projection - previous_projection; + + traction.mul(covariant_metric_tensor, increment_projection, epsilon_t); + + auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); + Vector + previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]); + traction = previous_traction + traction;*/ + + // method from Schweizerhof + auto covariant_metric_tensor = + GeometryUtils::covariantMetricTensor(covariant_basis); + + auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); + Vector previous_traction( + previous_tangential_tractions.begin(surface_dimension)[element.slave]); + + auto & previous_tangents = model.getPreviousTangents(); + Matrix previous_covariant_basis(previous_tangents.begin( + surface_dimension, spatial_dimension)[element.slave]); + auto previous_contravariant_metric_tensor = + GeometryUtils::contravariantMetricTensor(previous_covariant_basis); + + auto current_tangent = covariant_basis.transpose(); + auto previous_tangent = previous_covariant_basis.transpose(); + + for (auto alpha : arange(surface_dimension)) { + Vector tangent_alpha(current_tangent(alpha)); + for (auto gamma : arange(surface_dimension)) { + for (auto beta : arange(surface_dimension)) { + Vector tangent_beta(previous_tangent(beta)); + auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha); + traction[alpha] += previous_traction[gamma] * + previous_contravariant_metric_tensor(gamma, beta) * + t_alpha_t_beta; + } + } + } + + auto & previous_master_elements = model.getPreviousMasterElements(); + auto & previous_element = previous_master_elements[element.slave]; + + Vector previous_real_projection(spatial_dimension); + GeometryUtils::realProjection( + model.getMesh(), model.getContactDetector().getPositions(), + previous_element, previous_projection, previous_real_projection); + + Vector current_real_projection(spatial_dimension); + GeometryUtils::realProjection( + model.getMesh(), model.getContactDetector().getPositions(), + element.master, current_projection, current_real_projection); + + auto increment_real = current_real_projection - previous_real_projection; + Vector increment_xi(surface_dimension); + + auto contravariant_metric_tensor = + GeometryUtils::contravariantMetricTensor(covariant_basis); + + // increment in natural coordinate + for (auto beta : arange(surface_dimension)) { + for (auto gamma : arange(surface_dimension)) { + auto temp = increment_real.dot(current_tangent(gamma)); + temp *= contravariant_metric_tensor(beta, gamma); + increment_xi[beta] += temp; + } + } + + Vector temp(surface_dimension); + temp.mul(covariant_metric_tensor, increment_xi, epsilon_t); + + traction -= temp; +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeStickTangentialTraction( + const ContactElement & /*element*/, Vector & traction_trial, + Vector & traction_tangential) { + traction_tangential = traction_trial; +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeSlipTangentialTraction( + const ContactElement & element, const Matrix & covariant_basis, + Vector & traction_trial, Vector & traction_tangential) { + UInt surface_dimension = spatial_dimension - 1; + + auto & gaps = model.getGaps(); + auto & gap = gaps.begin()[element.slave]; + + // compute norm of trial traction + Real traction_trial_norm = 0; + auto contravariant_metric_tensor = + GeometryUtils::contravariantMetricTensor(covariant_basis); + + for (auto alpha : arange(surface_dimension)) { + for (auto beta : arange(surface_dimension)) { + traction_trial_norm += traction_trial[alpha] * traction_trial[beta] * + contravariant_metric_tensor(alpha, beta); + } + } + traction_trial_norm = sqrt(traction_trial_norm); + + auto slip_direction = traction_trial; + slip_direction /= traction_trial_norm; + + Real p_n = computeNormalTraction(gap); + traction_tangential = slip_direction; + traction_tangential *= mu * p_n; +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeNormalModuli(const ContactElement & element, + Matrix & stiffness) { + + auto surface_dimension = spatial_dimension - 1; + + auto & gaps = model.getGaps(); + Real gap(gaps.begin()[element.slave]); + + auto & projections = model.getProjections(); + Vector projection(projections.begin(surface_dimension)[element.slave]); + + auto & nodal_areas = model.getNodalArea(); + auto & nodal_area = nodal_areas.begin()[element.slave]; + + auto & normals = model.getNormals(); + Vector normal(normals.begin(spatial_dimension)[element.slave]); + + // method from Schweizerhof and A. Konyukhov, K. Schweizerhof + // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 + + // construct A matrix + const ElementType & type = element.master.type; + auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); + + UInt nb_nodes_per_contact = element.getNbNodes(); + Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); + + for (auto i : arange(nb_nodes_per_contact)) { + for (auto j : arange(spatial_dimension)) { + if (i == 0) { + A(j, i * spatial_dimension + j) = 1; + continue; + } + A(j, i * spatial_dimension + j) = -shapes[i - 1]; + } + } + + // construct the main part of normal matrix + Matrix k_main(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + + Matrix n_outer_n(spatial_dimension, spatial_dimension); + Matrix mat_n(normal.storage(), normal.size(), 1.); + n_outer_n.mul(mat_n, mat_n); + + Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); + tmp.mul(n_outer_n, A); + + k_main.mul(A, tmp); + k_main *= epsilon_n * heaviside(gap) * penetrationDerivative(gap) * nodal_area; + + // construct the rotational part of the normal matrix + auto & tangents = model.getTangents(); + Matrix covariant_basis( + tangents.begin(surface_dimension, spatial_dimension)[element.slave]); + + auto contravariant_metric_tensor = + GeometryUtils::contravariantMetricTensor(covariant_basis); + + // computing shape derivatives + auto && shape_derivatives = + ElementClassHelper<_ek_regular>::getDNDS(projection, type); + + // consists of 2 rotational parts + Matrix k_rot1(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + Matrix k_rot2(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); + + auto construct_Aj = [&](auto && dnds) { + for (auto i : arange(nb_nodes_per_contact)) { + for (auto j : arange(spatial_dimension)) { + if (i == 0) { + Aj(j, i * spatial_dimension + j) = 0; + continue; + } + Aj(j, i * spatial_dimension + j) = dnds(i - 1); + } + } + }; + + for (auto && values1 : enumerate(covariant_basis.transpose())) { + auto & alpha = std::get<0>(values1); + auto & tangent = std::get<1>(values1); + + Matrix n_outer_t(spatial_dimension, spatial_dimension); + Matrix mat_t(tangent.storage(), tangent.size(), 1.); + n_outer_t.mul(mat_n, mat_t); + + Matrix t_outer_n(spatial_dimension, spatial_dimension); + t_outer_n.mul(mat_t, mat_n); + + for (auto && values2 : enumerate(shape_derivatives.transpose())) { + auto & beta = std::get<0>(values2); + auto & dnds = std::get<1>(values2); + // construct Aj from shape function wrt to jth natural + // coordinate + construct_Aj(dnds); + + Matrix tmp(spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + Matrix tmp1(nb_nodes_per_contact * spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + tmp.mul(n_outer_t, A); + tmp1.mul(Aj, tmp); + tmp1 *= contravariant_metric_tensor(alpha, beta); + k_rot1 += tmp1; + + tmp.mul(t_outer_n, Aj); + tmp1.mul(A, tmp); + tmp1 *= contravariant_metric_tensor(alpha, beta); + k_rot2 += tmp1; + } + } + + k_rot1 *= -epsilon_n * heaviside(gap) * penetrationFunction(gap) * nodal_area; + k_rot2 *= -epsilon_n * heaviside(gap) * penetrationFunction(gap) * nodal_area; + + stiffness += k_main + k_rot1 + k_rot2; +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeTangentialModuli(const ContactElement & element, + Matrix & stiffness) { + if (mu == 0) { + return; + } + + stiffness.zero(); + + auto & contact_state = model.getContactState(); + auto state = contact_state.begin()[element.slave]; + + switch (state) { + case ContactState::_stick: { + computeStickModuli(element, stiffness); + break; + } + case ContactState::_slip: { + computeSlipModuli(element, stiffness); + break; + } + default: + break; + } +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeStickModuli(const ContactElement & element, + Matrix & stiffness) { + + auto surface_dimension = spatial_dimension - 1; + + auto & projections = model.getProjections(); + Vector projection(projections.begin(surface_dimension)[element.slave]); + + auto & nodal_areas = model.getNodalArea(); + auto & nodal_area = nodal_areas.begin()[element.slave]; + + // method from Schweizerhof and A. Konyukhov, K. Schweizerhof + // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 + + // construct A matrix + const ElementType & type = element.master.type; + auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); + + UInt nb_nodes_per_contact = element.getNbNodes(); + Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); + + for (auto i : arange(nb_nodes_per_contact)) { + for (auto j : arange(spatial_dimension)) { + if (i == 0) { + A(j, i * spatial_dimension + j) = 1; + continue; + } + A(j, i * spatial_dimension + j) = -shapes[i - 1]; + } + } + + // computing shape derivatives + auto && shape_derivatives = + ElementClassHelper<_ek_regular>::getDNDS(projection, type); + + Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); + + auto construct_Aj = [&](auto && dnds) { + for (auto i : arange(nb_nodes_per_contact)) { + for (auto j : arange(spatial_dimension)) { + if (i == 0) { + Aj(j, i * spatial_dimension + j) = 0; + continue; + } + Aj(j, i * spatial_dimension + j) = dnds(i - 1); + } + } + }; + + // tangents should have been calculated in normal modulii + auto & tangents = model.getTangents(); + Matrix covariant_basis( + tangents.begin(surface_dimension, spatial_dimension)[element.slave]); + + auto contravariant_metric_tensor = + GeometryUtils::contravariantMetricTensor(covariant_basis); + + // construct 1st part of the stick modulii + Matrix k_main(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + + for (auto && values1 : enumerate(covariant_basis.transpose())) { + auto & alpha = std::get<0>(values1); + auto & tangent_alpha = std::get<1>(values1); + + Matrix t_outer_t(spatial_dimension, spatial_dimension); + Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); + + for (auto && values2 : enumerate(covariant_basis.transpose())) { + auto & beta = std::get<0>(values2); + auto & tangent_beta = std::get<1>(values2); + + Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); + t_outer_t.mul(mat_t_alpha, mat_t_beta); + + Matrix tmp(spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + Matrix tmp1(nb_nodes_per_contact * spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + tmp.mul(t_outer_t, A); + tmp1.mul(A, tmp); + tmp1 *= contravariant_metric_tensor(alpha, beta); + k_main += tmp1; + } + } + + k_main *= -epsilon_t; + + // construct 2nd part of the stick modulii + auto & tangential_tractions = model.getTangentialTractions(); + Vector tangential_traction( + tangential_tractions.begin(surface_dimension)[element.slave]); + + Matrix k_second(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + + for (auto alpha : arange(surface_dimension)) { + + Matrix k_sum(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + + for (auto && values1 : enumerate(shape_derivatives.transpose())) { + auto & beta = std::get<0>(values1); + auto & dnds = std::get<1>(values1); + // construct Aj from shape function wrt to jth natural + // coordinate + construct_Aj(dnds); + for (auto && values2 : enumerate(covariant_basis.transpose())) { + auto & gamma = std::get<0>(values2); + auto & tangent_gamma = std::get<1>(values2); + + Matrix t_outer_t(spatial_dimension, spatial_dimension); + Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), + 1.); + + for (auto && values3 : enumerate(covariant_basis.transpose())) { + auto & theta = std::get<0>(values3); + auto & tangent_theta = std::get<1>(values3); + + Matrix mat_t_theta(tangent_theta.storage(), + tangent_theta.size(), 1.); + t_outer_t.mul(mat_t_gamma, mat_t_theta); + + Matrix tmp(spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + Matrix tmp1(nb_nodes_per_contact * spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + tmp.mul(t_outer_t, Aj); + tmp1.mul(A, tmp); + tmp1 *= contravariant_metric_tensor(alpha, theta) * + contravariant_metric_tensor(beta, gamma); + + Matrix tmp2(spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + Matrix tmp3(nb_nodes_per_contact * spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + tmp2.mul(t_outer_t, A); + tmp3.mul(Aj, tmp2); + tmp3 *= contravariant_metric_tensor(alpha, gamma) * + contravariant_metric_tensor(beta, theta); + + k_sum += tmp1 + tmp3; + } + } + } + + k_second += tangential_traction[alpha] * k_sum; + } + + stiffness += k_main * nodal_area - k_second * nodal_area; +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::computeSlipModuli(const ContactElement & element, + Matrix & stiffness) { + + auto surface_dimension = spatial_dimension - 1; + + auto & gaps = model.getGaps(); + Real gap(gaps.begin()[element.slave]); + + auto & nodal_areas = model.getNodalArea(); + auto & nodal_area = nodal_areas.begin()[element.slave]; + + // compute normal traction + Real p_n = computeNormalTraction(gap); + + auto & projections = model.getProjections(); + Vector projection(projections.begin(surface_dimension)[element.slave]); + + auto & normals = model.getNormals(); + Vector normal(normals.begin(spatial_dimension)[element.slave]); + + // restructure normal as a matrix for an outer product + Matrix mat_n(normal.storage(), normal.size(), 1.); + + // method from Schweizerhof and A. Konyukhov, K. Schweizerhof + // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 + + // construct A matrix + const ElementType & type = element.master.type; + auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); + + UInt nb_nodes_per_contact = element.getNbNodes(); + Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); + + for (auto i : arange(nb_nodes_per_contact)) { + for (auto j : arange(spatial_dimension)) { + if (i == 0) { + A(j, i * spatial_dimension + j) = 1; + continue; + } + A(j, i * spatial_dimension + j) = -shapes[i - 1]; + } + } + + // computing shape derivatives + auto && shape_derivatives = + ElementClassHelper<_ek_regular>::getDNDS(projection, type); + + Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); + + auto construct_Aj = [&](auto && dnds) { + for (auto i : arange(nb_nodes_per_contact)) { + for (auto j : arange(spatial_dimension)) { + if (i == 0) { + Aj(j, i * spatial_dimension + j) = 0; + continue; + } + Aj(j, i * spatial_dimension + j) = dnds(i - 1); + } + } + }; + + // tangents should have been calculated in normal modulii + auto & tangents = model.getTangents(); + Matrix covariant_basis( + tangents.begin(surface_dimension, spatial_dimension)[element.slave]); + + auto & tangential_tractions = model.getTangentialTractions(); + Vector tangential_traction( + tangential_tractions.begin(surface_dimension)[element.slave]); + + // compute norm of trial traction + Real traction_norm = 0; + auto contravariant_metric_tensor = + GeometryUtils::contravariantMetricTensor(covariant_basis); + + for (auto i : arange(surface_dimension)) { + for (auto j : arange(surface_dimension)) { + traction_norm += tangential_traction[i] * tangential_traction[j] * + contravariant_metric_tensor(i, j); + } + } + traction_norm = sqrt(traction_norm); + + // construct four parts of stick modulii (eq 107,107a-c) + Matrix k_first(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + Matrix k_second(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + Matrix k_third(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + Matrix k_fourth(nb_nodes_per_contact * spatial_dimension, + nb_nodes_per_contact * spatial_dimension); + + for (auto && values1 : enumerate(covariant_basis.transpose())) { + auto & alpha = std::get<0>(values1); + auto & tangent_alpha = std::get<1>(values1); + + Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); + + Matrix t_outer_n(spatial_dimension, spatial_dimension); + Matrix t_outer_t(spatial_dimension, spatial_dimension); + + for (auto && values2 : + zip(arange(surface_dimension), covariant_basis.transpose(), + shape_derivatives.transpose())) { + auto & beta = std::get<0>(values2); + auto & tangent_beta = std::get<1>(values2); + auto & dnds = std::get<2>(values2); + // construct Aj from shape function wrt to jth natural + // coordinate + construct_Aj(dnds); + + // eq 107 + Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); + t_outer_n.mul(mat_t_beta, mat_n); + + Matrix tmp(spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + Matrix tmp1(nb_nodes_per_contact * spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + tmp.mul(t_outer_n, A); + tmp1.mul(A, tmp); + + tmp1 *= epsilon_n * mu * tangential_traction[alpha] * + contravariant_metric_tensor(alpha, beta); + tmp1 /= traction_norm; + + k_first += tmp1 * nodal_area; + + // eq 107a + t_outer_t.mul(mat_t_alpha, mat_t_beta); + + tmp.mul(t_outer_t, A); + tmp1.mul(A, tmp); + + tmp1 *= epsilon_t * mu * p_n * contravariant_metric_tensor(alpha, beta); + tmp1 /= traction_norm; + + k_second += tmp1 * nodal_area; + + for (auto && values3 : enumerate(covariant_basis.transpose())) { + auto & gamma = std::get<0>(values3); + auto & tangent_gamma = std::get<1>(values3); + + Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), + 1.); + + for (auto && values4 : enumerate(covariant_basis.transpose())) { + auto & theta = std::get<0>(values4); + auto & tangent_theta = std::get<1>(values4); + + Matrix mat_t_theta(tangent_theta.storage(), + tangent_theta.size(), 1.); + t_outer_t.mul(mat_t_gamma, mat_t_theta); + + // eq 107b + tmp.mul(t_outer_t, A); + tmp1.mul(A, tmp); + + tmp1 *= epsilon_t * mu * p_n * tangential_traction[alpha] * + tangential_traction[beta]; + tmp1 *= contravariant_metric_tensor(alpha, gamma) * + contravariant_metric_tensor(beta, theta); + tmp1 /= pow(traction_norm, 3); + + k_third += tmp1 * nodal_area; + + // eq 107c + tmp.mul(t_outer_t, Aj); + tmp1.mul(A, tmp); + tmp1 *= contravariant_metric_tensor(alpha, theta) * + contravariant_metric_tensor(beta, gamma); + tmp1 *= mu * p_n * tangential_traction[alpha]; + tmp1 /= traction_norm; + + Matrix tmp2(spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + Matrix tmp3(nb_nodes_per_contact * spatial_dimension, + spatial_dimension * nb_nodes_per_contact); + tmp2.mul(t_outer_t, A); + tmp3.mul(Aj, tmp2); + tmp3 *= contravariant_metric_tensor(alpha, gamma) * + contravariant_metric_tensor(beta, theta); + tmp3 *= mu * p_n * tangential_traction[alpha]; + tmp3 /= traction_norm; + + k_fourth += (tmp1 + tmp3) * nodal_area; + } + } + } + } + + stiffness += k_third + k_fourth - k_first - k_second; +} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::beforeSolveStep() {} + +/* -------------------------------------------------------------------------- */ +template +void ResolutionPenaltyTemplate::afterSolveStep(__attribute__((unused)) bool converged) {} + +} // namespace akantu + +/* -------------------------------------------------------------------------- */ +namespace akantu { + +// Define function and derivative for linear penalty method +template T linearPenetrationFunction(T var) { return var; } +template T linearPenetrationDerivative(T var) { return 1.0; } + +// Instantiate linear penalty as a resolution +INSTANTIATE_RESOLUTION(penalty_linear, ResolutionPenaltyTemplate, Real, + linearPenetrationFunction, linearPenetrationDerivative); + +} // namespace akantu + +#endif // AKANTU_RESOLUTION_PENALTY_TMPL_HH