diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc index cbe3d8e2f..4d72d20a9 100644 --- a/src/model/contact_mechanics/resolution.cc +++ b/src/model/contact_mechanics/resolution.cc @@ -1,228 +1,228 @@ /** * @file resolution.cc * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Implementation of common part of the contact resolution class * * @section LICENSE * * Copyright (©) 2010-2018 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) : Memory(id, model.getMemoryID()), Parsable(ParserType::_contact_resolution, id), fem(model.getFEEngine()), name(""), 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; for (Int i = 0; i < indent; i++, space += 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 (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; 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); Real alpha{1.0}; - if(slave_result != std::end(slave_list) and master_result != std::end(master_list)) + if(is_master_deformable) alpha = 0.5; 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 (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