diff --git a/src/model/common/integration_scheme/newmark-beta.hh b/src/model/common/integration_scheme/newmark-beta.hh index 237310f75..39f460001 100644 --- a/src/model/common/integration_scheme/newmark-beta.hh +++ b/src/model/common/integration_scheme/newmark-beta.hh @@ -1,194 +1,194 @@ /** * @file newmark-beta.hh * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Tue Oct 05 2010 * @date last modification: Wed Jan 31 2018 * * @brief implementation of the newmark-@f$\beta@f$ integration scheme. This * implementation is taken from Méthodes numériques en mécanique des solides by * Alain Curnier \note{ISBN: 2-88074-247-1} * * * 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 "integration_scheme_2nd_order.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NEWMARK_BETA_HH__ #define __AKANTU_NEWMARK_BETA_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { /** * The three differentiate equations (dynamic and cinematic) are : * \f{eqnarray*}{ * M \ddot{u}_{n+1} + C \dot{u}_{n+1} + K u_{n+1} &=& q_{n+1} \\ * u_{n+1} &=& u_{n} + (1 - \alpha) \Delta t \dot{u}_{n} + \alpha \Delta t *\dot{u}_{n+1} + (1/2 - \alpha) \Delta t^2 \ddot{u}_n \\ * \dot{u}_{n+1} &=& \dot{u}_{n} + (1 - \beta) \Delta t \ddot{u}_{n} + \beta *\Delta t \ddot{u}_{n+1} * \f} * * Predictor: * \f{eqnarray*}{ * u^{0}_{n+1} &=& u_{n} + \Delta t \dot{u}_n + \frac{\Delta t^2}{2} *\ddot{u}_n \\ * \dot{u}^{0}_{n+1} &=& \dot{u}_{n} + \Delta t \ddot{u}_{n} \\ * \ddot{u}^{0}_{n+1} &=& \ddot{u}_{n} * \f} * * Solve : * \f[ (c M + d C + e K^i_{n+1}) w = = q_{n+1} - f^i_{n+1} - C \dot{u}^i_{n+1} *- M \ddot{u}^i_{n+1} \f] * * Corrector : * \f{eqnarray*}{ * \ddot{u}^{i+1}_{n+1} &=& \ddot{u}^{i}_{n+1} + c w \\ * \dot{u}^{i+1}_{n+1} &=& \dot{u}^{i}_{n+1} + d w \\ * u^{i+1}_{n+1} &=& u^{i}_{n+1} + e w * \f} * * c, d and e are parameters depending on the method used to solve the equations *\n * For acceleration : \f$ w = \delta \ddot{u}, e = \alpha \beta \Delta t^2, d = *\beta \Delta t, c = 1 \f$ \n * For velocity : \f$ w = \delta \dot{u}, e = 1/\beta \Delta t, d = *1, c = \alpha \Delta t \f$ \n * For displacement : \f$ w = \delta u, e = 1, d = *1/\alpha \Delta t, c = 1/\alpha \beta \Delta t^2 \f$ */ class NewmarkBeta : public IntegrationScheme2ndOrder { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NewmarkBeta(DOFManager & dof_manager, const ID & dof_id, Real alpha = 0., Real beta = 0.); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void predictor(Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs) const override; void corrector(const SolutionType & type, Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs, const Array & delta) const override; void assembleJacobian(const SolutionType & type, Real delta_t) override; public: Real getAccelerationCoefficient(const SolutionType & type, Real delta_t) const override; Real getVelocityCoefficient(const SolutionType & type, Real delta_t) const override; Real getDisplacementCoefficient(const SolutionType & type, Real delta_t) const override; private: template void allCorrector(Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs, const Array & delta) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Beta, beta, Real); AKANTU_GET_MACRO(Alpha, alpha, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the \f$\beta\f$ parameter Real beta; /// the \f$\alpha\f$ parameter Real alpha; Real k; Real h; /// last release of M matrix UInt m_release; /// last release of K matrix UInt k_release; /// last release of C matrix UInt c_release; }; /** * central difference method (explicit) * undamped stability condition : * \f$ \Delta t = \alpha \Delta t_{crit} = \frac{2}{\omega_{max}} \leq \min_{e} *\frac{l_e}{c_e}\f$ * */ class CentralDifference : public NewmarkBeta { public: CentralDifference(DOFManager & dof_manager, const ID & dof_id) - : NewmarkBeta(dof_manager, dof_id, 0., 1. / 2.){}; + : NewmarkBeta(dof_manager, dof_id, 0., 1./2.){}; std::vector getNeededMatrixList() override { return {"M", "C"}; } }; //#include "integration_scheme/central_difference.hh" /// undamped trapezoidal rule (implicit) class TrapezoidalRule2 : public NewmarkBeta { public: TrapezoidalRule2(DOFManager & dof_manager, const ID & dof_id) : NewmarkBeta(dof_manager, dof_id, 1. / 2., 1. / 2.){}; }; /// Fox-Goodwin rule (implicit) class FoxGoodwin : public NewmarkBeta { public: FoxGoodwin(DOFManager & dof_manager, const ID & dof_id) : NewmarkBeta(dof_manager, dof_id, 1. / 6., 1. / 2.){}; }; /// Linear acceleration (implicit) class LinearAceleration : public NewmarkBeta { public: LinearAceleration(DOFManager & dof_manager, const ID & dof_id) : NewmarkBeta(dof_manager, dof_id, 1. / 3., 1. / 2.){}; }; /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_NEWMARK_BETA_HH__ */ diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index b3bdd3753..c89767d35 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,693 +1,713 @@ /** * @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.hh" #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"); registerFEEngineObject("ContactFacetsFEEngine", mesh, Model::spatial_dimension - 1); 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() { AKANTU_DEBUG_IN(); getFEEngine("ContactMechanicsModel").initShapeFunctions(_not_ghost); getFEEngine("ContactMechanicsModel").initShapeFunctions(_ghost); getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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->contact_state, 1, "contact_state"); this->allocNodalField(this->previous_master_elements, 1, "previous_master_elements"); this->allocNodalField(this->normals, spatial_dimension, "normals"); auto surface_dimension = spatial_dimension - 1; this->allocNodalField(this->tangents, surface_dimension*spatial_dimension, "tangents"); this->allocNodalField(this->projections, surface_dimension, "projections"); this->allocNodalField(this->previous_projections, surface_dimension, "previous_projections"); this->allocNodalField(this->previous_tangents, surface_dimension*spatial_dimension, "previous_tangents"); this->allocNodalField(this->tangential_tractions, surface_dimension, "tangential_tractions"); this->allocNodalField(this->previous_tangential_tractions, surface_dimension, "previous_tangential_tractions"); // todo register multipliers as dofs for lagrange multipliers } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", 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); // TODO : uncomment when developing code for parallelization, // currently it addes the force twice for not ghost elements // hence source of error AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { // save the previous state this->savePreviousState(); contact_elements.clear(); contact_elements.resize(0); // this needs to be resized if cohesive elements are added UInt nb_nodes = mesh.getNbNodes(); auto resize_arrays = [&](auto & internal_array) { internal_array->clear(); internal_array->resize(nb_nodes, 0.); }; resize_arrays(gaps); resize_arrays(normals); resize_arrays(tangents); resize_arrays(projections); resize_arrays(tangential_tractions); resize_arrays(contact_state); resize_arrays(nodal_area); resize_arrays(external_force); this->detector->search(contact_elements, *gaps, *normals, *projections); // intepenetration value must be positive for contact mechanics // model to work by default the gap value from detector is negative std::for_each((*gaps).begin(), (*gaps).end(), [](Real & gap){ gap *= -1.; }); if (contact_elements.size() != 0) { this->computeNodalAreas(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::savePreviousState() { AKANTU_DEBUG_IN(); // saving previous natural projections previous_projections->clear(); previous_projections->resize(projections->size(), 0.); (*previous_projections).copy(*projections); // saving previous tangents previous_tangents->clear(); previous_tangents->resize(tangents->size(), 0.); (*previous_tangents).copy(*tangents); // saving previous tangential traction previous_tangential_tractions->clear(); previous_tangential_tractions->resize(tangential_tractions->size(), 0.); (*previous_tangential_tractions).copy(*tangential_tractions); previous_master_elements->clear(); previous_master_elements->resize(projections->size()); for (auto & element : contact_elements) { (*previous_master_elements)[element.slave] = element.master; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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 = getFEEngineClassBoundary("ContactMechanicsModel"); fem_boundary.initShapeFunctions(getPositions(), _not_ghost); fem_boundary.initShapeFunctions(getPositions(), _ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); switch (spatial_dimension) { case 1: { std::fill((*nodal_area).begin(), (*nodal_area).end(), 1.); break; } case 2: case 3: { this->applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), mesh.getElementGroup("contact_surface")); - for (auto && tuple : + /*for (auto && tuple : zip(*nodal_area, make_view(*external_force, spatial_dimension))) { auto & area = std::get<0>(tuple); Vector force(std::get<1>(tuple)); area = force.norm(); + }*/ + + + for (auto && tuple : + zip(*nodal_area, + make_view(*external_force, spatial_dimension), + make_view(*normals, 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);*/ + auto & area = std::get<0>(tuple); + Vector force(std::get<1>(tuple)); + Vector normal(std::get<2>(tuple)); + area = abs(force.dot(normal)); } + break; } default: break; } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ 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(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() { for (auto & resolution : resolutions) resolution->beforeSolveStep(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep(bool converged) { for (auto & resolution : resolutions) resolution->afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string &, const std::string &, bool) { 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["contact_state"] = this->contact_state; real_nodal_fields["tangential_traction"] = this->tangential_tractions; 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::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_contact.cc b/src/model/model_couplers/coupler_solid_contact.cc index 869fc3125..a17ef1ff5 100644 --- a/src/model/model_couplers/coupler_solid_contact.cc +++ b/src/model/model_couplers/coupler_solid_contact.cc @@ -1,538 +1,482 @@ /** * @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); solid->initFull( _analysis_method = this->method); contact->initFull(_analysis_method = this->method); } /* -------------------------------------------------------------------------- */ 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 time_step_solver_type, NonLinearSolverType non_linear_solver_type) { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); auto & contact_model_solver = aka::as_type(*contact); contact_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); } /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidContact::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", 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 + + switch (method) { + case _explicit_lumped_mass: { + auto & current_positions = contact->getContactDetector().getPositions(); + current_positions.copy(solid->getCurrentPosition()); + contact->search(); + break; + } + default: + break; + } + 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); this->getDOFManager().assembleToResidual("displacement", contact_force, 1); } /* -------------------------------------------------------------------------- */ 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); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::predictor() { - - switch (method) { - case _explicit_lumped_mass: { - auto & current_positions = contact->getContactDetector().getPositions(); - current_positions.copy(solid->getCurrentPosition()); - contact->search(); - break; - } - default: - break; - } auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.predictor(); - - - - + } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::corrector() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.corrector(); switch (method) { case _static: case _implicit_dynamic: { auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(solid->getCurrentPosition()); 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::beforeSolveStep() { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.beforeSolveStep(); auto & contact_solver_callback = aka::as_type(*contact); contact_solver_callback.beforeSolveStep(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::afterSolveStep(bool converged) { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.afterSolveStep(converged); auto & contact_solver_callback = aka::as_type(*contact); contact_solver_callback.afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ 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(true); switch (method) { case _static: case _implicit_dynamic: { 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) { 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 == "contact_state" or field_name == "gaps" or field_name == "previous_gaps" or field_name == "areas" or field_name == "tangents") { field = contact->createNodalFieldReal(field_name, group_name, padding_flag); } else { field = solid->createNodalFieldReal(field_name, group_name, padding_flag); } return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool 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