diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index dca325db0..fc197e906 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,623 +1,630 @@ /** * @file coontact_mechanics_model.cc * * @author Mohit Pundir * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Contact mechanics model * * @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 "contact_mechanics_model.hh" #include "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" +#include "group_manager_inline_impl.cc" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id, memory_id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, id + ":contact_detector"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::~ContactMechanicsModel() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // initalize the resolutions if (this->parser.getLastParsedFile() != "") { this->instantiateResolutions(); } this->initResolutions(); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::instantiateResolutions() { ParserSection model_section; bool is_empty; std::tie(model_section, is_empty) = this->getParserSection(); if (not is_empty) { auto model_resolutions = model_section.getSubSections(ParserType::_contact_resolution); for (const auto & section : model_resolutions) { this->registerNewResolution(section); } } auto sub_sections = this->parser.getSubSections(ParserType::_contact_resolution); for (const auto & section : sub_sections) { this->registerNewResolution(section); } if (resolutions.empty()) AKANTU_EXCEPTION("No contact resolutions where instantiated for the model" << getID()); are_resolutions_instantiated = true; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution(const ParserSection & section) { std::string res_name; std::string res_type = section.getName(); std::string opt_param = section.getOption(); try { std::string tmp = section.getParameter("name"); res_name = tmp; /** this can seem weird, but there is an ambiguous operator * overload that i couldn't solve. @todo remove the * weirdness of this code */ } catch (debug::Exception &) { AKANTU_ERROR("A contact resolution of type \'" << res_type << "\' in the input file has been defined without a name!"); } Resolution & res = this->registerNewResolution(res_name, res_type, opt_param); res.parseSection(section); return res; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution( const ID & res_name, const ID & res_type, const ID & opt_param) { AKANTU_DEBUG_ASSERT(resolutions_names_to_id.find(res_name) == resolutions_names_to_id.end(), "A resolution with this name '" << res_name << "' has already been registered. " << "Please use unique names for resolutions"); UInt res_count = resolutions.size(); resolutions_names_to_id[res_name] = res_count; std::stringstream sstr_res; sstr_res << this->id << ":" << res_count << ":" << res_type; ID res_id = sstr_res.str(); std::unique_ptr resolution = ResolutionFactory::getInstance().allocate(res_type, spatial_dimension, opt_param, *this, res_id); resolutions.push_back(std::move(resolution)); return *(resolutions.back()); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initResolutions() { AKANTU_DEBUG_ASSERT(resolutions.size() != 0, "No resolutions to initialize !"); if (!are_resolutions_instantiated) instantiateResolutions(); // \TODO check if each resolution needs a initResolution() method } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initSolver( TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType) { // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->normal_force, spatial_dimension, "normal_force"); this->allocNodalField(this->tangential_force, spatial_dimension, "tangential_force"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->nodal_area, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); this->allocNodalField(this->stick_or_slip, 1, "stick_or_slip"); this->allocNodalField(this->normals, spatial_dimension, "normals"); this->allocNodalField(this->tangents, spatial_dimension, "tangents"); this->allocNodalField(this->projections, spatial_dimension - 1, "projections"); // todo register multipliers as dofs for lagrange multipliers } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", TimeStepSolverType::_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", TimeStepSolverType::_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", TimeStepSolverType::_dynamic_lumped); break; } default: return std::make_tuple("unkown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the contact forces"); UInt nb_nodes = mesh.getNbNodes(); this->internal_force->clear(); this->normal_force->clear(); this->tangential_force->clear(); internal_force->resize(nb_nodes, 0.); normal_force->resize(nb_nodes, 0.); tangential_force->resize(nb_nodes, 0.); // assemble the forces due to contact auto assemble = [&](auto && ghost_type) { for (auto & resolution : resolutions) { resolution->assembleInternalForces(ghost_type); } }; AKANTU_DEBUG_INFO("Assemble residual for local elements"); assemble(_not_ghost); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); assemble(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { UInt nb_nodes = mesh.getNbNodes(); contact_elements.clear(); contact_elements.resize(0); // this to resize if cohesive elements are added gaps->clear(); gaps->resize(nb_nodes, 0.); normals->clear(); normals->resize(nb_nodes, 0.); projections->clear(); projections->resize(nb_nodes, 0.); this->detector->search(contact_elements, *gaps, *normals, *projections); for (auto & gap : *gaps) { if (gap < 0) gap = std::abs(gap); else gap = -gap; } this->computeNodalAreas(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleFieldsFromContactMap() { UInt nb_nodes = mesh.getNbNodes(); this->gaps->clear(); gaps->resize(nb_nodes, 0.); normals->resize(nb_nodes, 0.); tangents->resize(nb_nodes, 0.); if (this->contact_map.empty()) return; for (auto & entry : contact_map) { const auto & element = entry.second; auto connectivity = element.connectivity; auto node = connectivity(0); (*gaps)[node] = element.gap; for (UInt i = 0; i < spatial_dimension; ++i) { (*normals)(node, i) = element.normal[i]; (*tangents)(node, i) = element.tangents(0, i); } } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::computeNodalAreas() { UInt nb_nodes = mesh.getNbNodes(); nodal_area->clear(); external_force->clear(); nodal_area->resize(nb_nodes, 0.); external_force->resize(nb_nodes, 0.); auto & fem_boundary = this->getFEEngineBoundary(); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); this->applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), mesh.getElementGroup("contact_surface")); for (auto && tuple : zip(*nodal_area, make_view(*external_force, spatial_dimension))) { auto & area = std::get<0>(tuple); auto & force = std::get<1>(tuple); for (auto & f : force) area += pow(f, 2); area = sqrt(area); } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep() {} /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Contact Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + resolutions [" << std::endl; for (auto & resolution : resolutions) { resolution->printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ MatrixType ContactMechanicsModel::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } for (auto & resolution : resolutions) { resolution->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER +/* -------------------------------------------------------------------------- */ std::shared_ptr -ContactMechanicsModel::createNodalFieldBool(const std::string & /*field_name*/, - const std::string & /*group_name*/, - bool /*padding_flag*/) { +ContactMechanicsModel::createNodalFieldBool(const std::string &, + const std::string &, bool) { - std::shared_ptr field; - return field; + return nullptr; } + /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["contact_force"] = this->internal_force; real_nodal_fields["normal_force"] = this->normal_force; real_nodal_fields["tangential_force"] = this->tangential_force; real_nodal_fields["blocked_dofs"] = this->blocked_dofs; real_nodal_fields["normals"] = this->normals; real_nodal_fields["tangents"] = this->tangents; real_nodal_fields["gaps"] = this->gaps; real_nodal_fields["areas"] = this->nodal_area; real_nodal_fields["stick_or_slip"] = this->stick_or_slip; - if (padding_flag) { - return this->mesh.createNodalField(real_nodal_fields[field_name], - group_name, 3); - } else { - return this->mesh.createNodalField(real_nodal_fields[field_name], - group_name); - } - std::shared_ptr field; + if (padding_flag) + field = this->mesh.createNodalField(real_nodal_fields[field_name], + group_name, 3); + else + field = this->mesh.createNodalField(real_nodal_fields[field_name], + group_name); return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr -ContactMechanicsModel::createNodalFieldReal(const std::string & /*field_name*/, - const std::string & /*group_name*/, - bool /*padding_flag*/) { +ContactMechanicsModel::createNodalFieldBool(const std::string &, + const std::string &, bool) { + + return nullptr; +} + + +/* -------------------------------------------------------------------------- */ +std::shared_ptr +ContactMechanicsModel::createNodalFieldReal(const std::string & , + const std::string & , bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData( const Array & elements, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData( const Array & dofs, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); UInt size = 0; AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/, const Array & /*dofs*/, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/, const Array & /*dofs*/, const SynchronizationTag & /*tag*/) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_cohesive_contact.cc b/src/model/model_couplers/coupler_solid_cohesive_contact.cc index da2e01f6b..ed632e883 100644 --- a/src/model/model_couplers/coupler_solid_cohesive_contact.cc +++ b/src/model/model_couplers/coupler_solid_cohesive_contact.cc @@ -1,529 +1,527 @@ /** * @file coupler_solid_cohesive_contact.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 * @date last modification: Thu May 22 2019 * * @brief class for coupling of solid mechanics and conatct mechanics * model * * @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 "coupler_solid_cohesive_contact.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidCohesiveContact::CouplerSolidCohesiveContact( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject( "CohesiveFEEngine", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_cohesive_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_cohesive_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModelCohesive(mesh, Model::spatial_dimension, "solid_mechanics_model_cohesive", 0, this->dof_manager); contact = new ContactMechanicsModel(mesh, Model::spatial_dimension, "contact_mechanics_model", 0, this->dof_manager); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidCohesiveContact::~CouplerSolidCohesiveContact() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initModel() { getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidCohesiveContact::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initSolver(TimeStepSolverType, NonLinearSolverType) {} /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidCohesiveContact::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", TimeStepSolverType::_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", TimeStepSolverType::_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", TimeStepSolverType::_dynamic_lumped); break; } default: return std::make_tuple("unkown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidCohesiveContact::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidCohesiveContact::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleResidual() { // computes the internal forces this->assembleInternalForces(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); auto & contact_map = contact->getContactMap(); switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { for (auto & pair : contact_map) { auto & connectivity = pair.second.connectivity; for (auto node : connectivity) { for (auto s : arange(spatial_dimension)) external_force(node, s) = contact_force(node, s); } } break; } default: break; } /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); auto & contact_map = contact->getContactMap(); switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { for (auto & pair : contact_map) { auto & connectivity = pair.second.connectivity; for (auto node : connectivity) { for (auto s : arange(spatial_dimension)) external_force(node, s) = contact_force(node, s); } } break; } default: break; } if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::afterSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::predictor() { switch (method) { case _explicit_dynamic_contact: { Array displacement(0, Model::spatial_dimension); auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(mesh.getNodes()); auto us = this->getDOFManager().getDOFs("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); if (not bld) cp += u; } contact->search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::corrector() { switch (method) { case _implicit_contact: case _explicit_contact: { Array displacement(0, Model::spatial_dimension); auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(mesh.getNodes()); auto us = this->getDOFManager().getDOFs("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); if (not bld) cp += u; } contact->search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidCohesiveContact::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); contact->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); switch (method) { case _implicit_contact: { contact->assembleStiffnessMatrix(); break; } default: break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMassLumped(GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMass(GhostType ghost_type) { solid->assembleMass(ghost_type); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { - return solid->createElementalField(field_name, group_name, padding_flag, - spatial_dimension, kind); - std::shared_ptr field; + field = solid->createElementalField(field_name, group_name, padding_flag, + spatial_dimension, kind); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldReal( const std::string & field_name, const std::string & group_name, bool padding_flag) { + std::shared_ptr field; if (field_name == "contact_force" or field_name == "normals" or field_name == "gaps" or field_name == "previous_gaps" or field_name == "areas" or field_name == "tangents") { - return contact->createNodalFieldReal(field_name, group_name, padding_flag); + field = contact->createNodalFieldReal(field_name, group_name, padding_flag); } else { - return solid->createNodalFieldReal(field_name, group_name, padding_flag); + field = solid->createNodalFieldReal(field_name, group_name, padding_flag); } - std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldBool( const std::string & field_name, const std::string & group_name, bool padding_flag) { - return solid->createNodalFieldBool(field_name, group_name, padding_flag); - std::shared_ptr field; + field = solid->createNodalFieldBool(field_name, group_name, padding_flag); return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createElementalField(const std::string &, const std::string &, bool, const UInt &, const ElementKind &) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_contact.cc b/src/model/model_couplers/coupler_solid_contact.cc index 38fd9fb00..c11e7996c 100644 --- a/src/model/model_couplers/coupler_solid_contact.cc +++ b/src/model/model_couplers/coupler_solid_contact.cc @@ -1,546 +1,546 @@ /** * @file coupler_solid_contact.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 * @date last modification: Thu May 22 2019 * * @brief class for coupling of solid mechanics and conatct mechanics * model * * @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 "coupler_solid_contact.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidContact::CouplerSolidContact( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id, memory_id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("CouplerSolidContact", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModel(mesh, Model::spatial_dimension, "solid_mechanics_model", 0, this->dof_manager); contact = new ContactMechanicsModel(mesh, Model::spatial_dimension, "contact_mechanics_model", 0, this->dof_manager); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidContact::~CouplerSolidContact() {} /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidContact::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initSolver(TimeStepSolverType, NonLinearSolverType) {} /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidContact::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", TimeStepSolverType::_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", TimeStepSolverType::_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", TimeStepSolverType::_dynamic_lumped); break; } default: return std::make_tuple("unkown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidContact::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidContact::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleResidual() { // computes the internal forces this->assembleInternalForces(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn(const_cast(mesh).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; }; switch(method) { case _explicit_contact: case _implicit_contact: case _explicit_dynamic_contact: { for (auto & element : contact->getContactElements()) { for (auto & conn : get_connectivity(element.slave, element.master)) { for (auto dim : arange(spatial_dimension)) { external_force(conn, dim) = contact_force(conn, dim); } } } } default: break; } /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _explicit_contact: case _implicit_contact: { //this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); //contact->assembleInternalForces(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); /*auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn(const_cast(mesh).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; }; switch(method) { case _explicit_dynamic_contact: { for (auto & element : contact->getContactElements()) { for (auto & conn : get_connectivity(element.slave, element.master)) { for (auto dim : arange(spatial_dimension)) { external_force(conn, dim) = contact_force(conn, dim); } } } } default: break; }*/ if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", contact_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _explicit_contact: case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidContact::afterSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidContact::predictor() { switch (method) { case _implicit_contact: case _explicit_dynamic_contact: { Array displacement(0, Model::spatial_dimension); auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(mesh.getNodes()); auto us = this->getDOFManager().getDOFs("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); //const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); //if (not bld) cp += u; } contact->search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::corrector() { switch (method) { case _explicit_contact: { Array displacement(0, Model::spatial_dimension); auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(mesh.getNodes()); auto us = this->getDOFManager().getDOFs("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); //const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); //if (not bld) cp += u; } contact->search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidContact::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); contact->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); switch (method) { case _explicit_contact: case _implicit_contact: { //contact->assembleStiffnessMatrix(); break; } default: break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMassLumped(GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMass(GhostType ghost_type) { solid->assembleMass(ghost_type); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { - return solid->createElementalField(field_name, group_name, padding_flag, - spatial_dimension, kind); - std::shared_ptr field; + field = solid->createElementalField(field_name, group_name, padding_flag, + spatial_dimension, kind); + return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { + std::shared_ptr field; if (field_name == "contact_force" or field_name == "normals" or field_name == "normal_force" or field_name == "tangential_force" or field_name == "stick_or_slip" or field_name == "gaps" or field_name == "previous_gaps" or field_name == "areas" or field_name == "tangents") { - return contact->createNodalFieldReal(field_name, group_name, padding_flag); + field = contact->createNodalFieldReal(field_name, group_name, padding_flag); } else { - return solid->createNodalFieldReal(field_name, group_name, padding_flag); + field = solid->createNodalFieldReal(field_name, group_name, padding_flag); } - std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { - return solid->createNodalFieldBool(field_name, group_name, padding_flag); - + std::shared_ptr field; + field = solid->createNodalFieldBool(field_name, group_name, padding_flag); return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createElementalField(const std::string &, const std::string &, bool, const UInt &, const ElementKind &) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidContact::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu