diff --git a/src/model/common/model_solver.cc b/src/model/common/model_solver.cc index 3c6de5ad3..3f3af09e9 100644 --- a/src/model/common/model_solver.cc +++ b/src/model/common/model_solver.cc @@ -1,404 +1,399 @@ /** * @file model_solver.cc * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Tue Mar 30 2021 * * @brief Implementation of ModelSolver * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "model_solver.hh" #include "dof_manager.hh" #include "dof_manager_default.hh" #include "mesh.hh" #include "non_linear_solver.hh" #include "time_step_solver.hh" #if defined(AKANTU_USE_PETSC) #include "dof_manager_petsc.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template static T getOptionToType(const std::string & opt_str) { std::stringstream sstr(opt_str); T opt; sstr >> opt; return opt; } /* -------------------------------------------------------------------------- */ ModelSolver::ModelSolver(Mesh & mesh, const ModelType & type, const ID & id) : Parsable(ParserType::_model, id), model_type(type), parent_id(id), mesh(mesh) {} -/* -------------------------------------------------------------------------- */ -ModelSolver::ModelSolver(Mesh & mesh, const ModelType & type, const ID & id, - std::shared_ptr dof_manager) - : ModelSolver(mesh, type, id) { - if (not dof_manager) { - this->initDOFManager(); - } else { - this->dof_manager = dof_manager; - this->setDOFManager(*this->dof_manager); - } -} - /* -------------------------------------------------------------------------- */ ModelSolver::~ModelSolver() = default; /* -------------------------------------------------------------------------- */ std::tuple ModelSolver::getParserSection() { auto sub_sections = getStaticParser().getSubSections(ParserType::_model); auto it = std::find_if( sub_sections.begin(), sub_sections.end(), [&](auto && section) { auto type = getOptionToType(section.getName()); // default id should be the model type if not defined std::string name = section.getParameter("name", this->parent_id); return type == model_type and name == this->parent_id; }); if (it == sub_sections.end()) { return std::make_tuple(ParserSection(), true); } return std::make_tuple(*it, false); } /* -------------------------------------------------------------------------- */ -std::shared_ptr ModelSolver::initDOFManager() { +std::shared_ptr +ModelSolver::initDOFManager(std::shared_ptr dof_manager) { + if (dof_manager) { + this->dof_manager = dof_manager; + this->setDOFManager(*this->dof_manager); + return this->dof_manager; + } + // default without external solver activated at compilation same as mumps that // is the historical solver but with only the lumped solver ID solver_type = "default"; #if defined(AKANTU_USE_MUMPS) solver_type = "default"; #elif defined(AKANTU_USE_PETSC) solver_type = "petsc"; #endif ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { solver_type = section.getOption(solver_type); return this->initDOFManager(section, solver_type); } return this->initDOFManager(solver_type); } /* -------------------------------------------------------------------------- */ std::shared_ptr ModelSolver::initDOFManager(const ID & solver_type) { if (dof_manager) { AKANTU_EXCEPTION("The DOF manager for this model is already initialized !"); } try { this->dof_manager = DOFManagerFactory::getInstance().allocate( solver_type, mesh, this->parent_id + ":dof_manager_" + solver_type); } catch (...) { AKANTU_EXCEPTION( "To use the solver " << solver_type << " you will have to code it. This is an unknown solver type."); } this->setDOFManager(*this->dof_manager); return this->dof_manager; } /* -------------------------------------------------------------------------- */ std::shared_ptr ModelSolver::initDOFManager(const ParserSection & section, const ID & solver_type) { this->initDOFManager(solver_type); auto sub_sections = section.getSubSections(ParserType::_time_step_solver); // parsing the time step solvers for (auto && section : sub_sections) { ID type = section.getName(); ID solver_id = section.getParameter("name", type); auto tss_type = getOptionToType(type); auto tss_options = this->getDefaultSolverOptions(tss_type); auto sub_solvers_sect = section.getSubSections(ParserType::_non_linear_solver); auto nb_non_linear_solver_section = section.getNbSubSections(ParserType::_non_linear_solver); auto nls_type = tss_options.non_linear_solver_type; if (nb_non_linear_solver_section == 1) { auto && nls_section = *(sub_solvers_sect.first); nls_type = getOptionToType(nls_section.getName()); } else if (nb_non_linear_solver_section > 0) { AKANTU_EXCEPTION("More than one non linear solver are provided for the " "time step solver " << solver_id); } this->getNewSolver(solver_id, tss_type, nls_type); if (nb_non_linear_solver_section == 1) { const auto & nls_section = *(sub_solvers_sect.first); this->dof_manager->getNonLinearSolver(solver_id).parseSection( nls_section); } auto sub_integrator_sections = section.getSubSections(ParserType::_integration_scheme); for (auto && is_section : sub_integrator_sections) { const auto & dof_type_str = is_section.getName(); ID dof_id; try { ID tmp = is_section.getParameter("name"); dof_id = tmp; } catch (...) { AKANTU_EXCEPTION("No degree of freedom name specified for the " "integration scheme of type " << dof_type_str); } auto it_type = getOptionToType(dof_type_str); IntegrationScheme::SolutionType s_type = is_section.getParameter( "solution_type", tss_options.solution_type[dof_id]); this->setIntegrationScheme(solver_id, dof_id, it_type, s_type); } for (auto & is_type : tss_options.integration_scheme_type) { if (!this->hasIntegrationScheme(solver_id, is_type.first)) { this->setIntegrationScheme(solver_id, is_type.first, is_type.second, tss_options.solution_type[is_type.first]); } } } if (section.hasParameter("default_solver")) { ID default_solver = section.getParameter("default_solver"); if (this->hasSolver(default_solver)) { this->setDefaultSolver(default_solver); } else { AKANTU_EXCEPTION( "The solver \"" << default_solver << "\" was not created, it cannot be set as default solver"); } } return this->dof_manager; } /* -------------------------------------------------------------------------- */ TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) { ID tmp_solver_id = solver_id; if (tmp_solver_id.empty()) { tmp_solver_id = this->default_solver_id; } TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id); return tss; } /* -------------------------------------------------------------------------- */ const TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) const { ID tmp_solver_id = solver_id; if (solver_id.empty()) { tmp_solver_id = this->default_solver_id; } const TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id); return tss; } /* -------------------------------------------------------------------------- */ TimeStepSolver & ModelSolver::getTimeStepSolver(const ID & solver_id) { return this->getSolver(solver_id); } /* -------------------------------------------------------------------------- */ const TimeStepSolver & ModelSolver::getTimeStepSolver(const ID & solver_id) const { return this->getSolver(solver_id); } /* -------------------------------------------------------------------------- */ NonLinearSolver & ModelSolver::getNonLinearSolver(const ID & solver_id) { return this->getSolver(solver_id).getNonLinearSolver(); } /* -------------------------------------------------------------------------- */ const NonLinearSolver & ModelSolver::getNonLinearSolver(const ID & solver_id) const { return this->getSolver(solver_id).getNonLinearSolver(); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasSolver(const ID & solver_id) const { ID tmp_solver_id = solver_id; if (solver_id.empty()) { tmp_solver_id = this->default_solver_id; } if (not this->dof_manager) { AKANTU_EXCEPTION("No DOF manager was initialized"); } return this->dof_manager->hasTimeStepSolver(tmp_solver_id); } /* -------------------------------------------------------------------------- */ void ModelSolver::setDefaultSolver(const ID & solver_id) { AKANTU_DEBUG_ASSERT( this->hasSolver(solver_id), "Cannot set the default solver to a solver that does not exists"); this->default_solver_id = solver_id; } /* -------------------------------------------------------------------------- */ void ModelSolver::solveStep(SolverCallback & callback, const ID & solver_id) { AKANTU_DEBUG_IN(); TimeStepSolver & tss = this->getSolver(solver_id); // make one non linear solve tss.solveStep(callback); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ModelSolver::solveStep(const ID & solver_id) { solveStep(*this, solver_id); } /* -------------------------------------------------------------------------- */ void ModelSolver::getNewSolver(const ID & solver_id, TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { if (this->default_solver_id.empty()) { this->default_solver_id = solver_id; } if (non_linear_solver_type == NonLinearSolverType::_auto) { switch (time_step_solver_type) { case TimeStepSolverType::_dynamic: case TimeStepSolverType::_static: non_linear_solver_type = NonLinearSolverType::_newton_raphson; break; case TimeStepSolverType::_dynamic_lumped: non_linear_solver_type = NonLinearSolverType::_lumped; break; case TimeStepSolverType::_not_defined: AKANTU_EXCEPTION(time_step_solver_type << " is not a valid time step solver type"); break; } } this->initSolver(time_step_solver_type, non_linear_solver_type); NonLinearSolver & nls = this->dof_manager->getNewNonLinearSolver( solver_id, non_linear_solver_type); this->dof_manager->getNewTimeStepSolver(solver_id, time_step_solver_type, nls, *this); } /* -------------------------------------------------------------------------- */ Real ModelSolver::getTimeStep(const ID & solver_id) const { const TimeStepSolver & tss = this->getSolver(solver_id); return tss.getTimeStep(); } /* -------------------------------------------------------------------------- */ void ModelSolver::setTimeStep(Real time_step, const ID & solver_id) { TimeStepSolver & tss = this->getSolver(solver_id); return tss.setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void ModelSolver::setIntegrationScheme( const ID & solver_id, const ID & dof_id, const IntegrationSchemeType & integration_scheme_type, IntegrationScheme::SolutionType solution_type) { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); tss.setIntegrationScheme(dof_id, integration_scheme_type, solution_type); } /* -------------------------------------------------------------------------- */ void ModelSolver::setIntegrationScheme( const ID & solver_id, const ID & dof_id, std::unique_ptr & integration_scheme, IntegrationScheme::SolutionType solution_type) { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); tss.setIntegrationScheme(dof_id, integration_scheme, solution_type); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasDefaultSolver() const { return (not this->default_solver_id.empty()); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasIntegrationScheme(const ID & solver_id, const ID & dof_id) const { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); return tss.hasIntegrationScheme(dof_id); } /* -------------------------------------------------------------------------- */ void ModelSolver::predictor() {} /* -------------------------------------------------------------------------- */ void ModelSolver::corrector() {} /* -------------------------------------------------------------------------- */ TimeStepSolverType ModelSolver::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions ModelSolver::getDefaultSolverOptions(__attribute__((unused)) const TimeStepSolverType & type) const { ModelSolverOptions options; options.non_linear_solver_type = NonLinearSolverType::_auto; return options; } } // namespace akantu diff --git a/src/model/common/model_solver.hh b/src/model/common/model_solver.hh index 3cb511baf..d2c41a0da 100644 --- a/src/model/common/model_solver.hh +++ b/src/model/common/model_solver.hh @@ -1,204 +1,203 @@ /** * @file model_solver.hh * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Fri May 14 2021 * * @brief Class regrouping the common solve interface to the different models * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "integration_scheme.hh" #include "parsable.hh" #include "solver_callback.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MODEL_SOLVER_HH_ #define AKANTU_MODEL_SOLVER_HH_ namespace akantu { class Mesh; class DOFManager; class TimeStepSolver; class NonLinearSolver; struct ModelSolverOptions; } // namespace akantu namespace akantu { class ModelSolver : public Parsable, public SolverCallback, public SynchronizerRegistry { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ModelSolver(Mesh & mesh, const ModelType & type, const ID & id); - ModelSolver(Mesh & mesh, const ModelType & type, const ID & id, - std::shared_ptr dof_manager); ~ModelSolver() override; /// initialize the dof manager based on solver type passed in the input file - std::shared_ptr initDOFManager(); + std::shared_ptr + initDOFManager(std::shared_ptr dof_manager = nullptr); /// initialize the dof manager based on the used chosen solver type std::shared_ptr initDOFManager(const ID & solver_type); protected: /// initialize the dof manager based on the used chosen solver type std::shared_ptr initDOFManager(const ParserSection & section, const ID & solver_type); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Callback for the model to instantiate the matricees when needed virtual void initSolver(TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType /*non_linear_solver_type*/) {} /// get the section in the input file (if it exsits) corresponding to this /// model std::tuple getParserSection(); /// solve a step using a given pre instantiated time step solver and /// non linear solver virtual void solveStep(const ID & solver_id = ""); /// solve a step using a given pre instantiated time step solver and /// non linear solver with a user defined callback instead of the /// model itself /!\ This can mess up everything virtual void solveStep(SolverCallback & callback, const ID & solver_id = ""); /// Initialize a time solver that can be used afterwards with its id void getNewSolver( const ID & solver_id, TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type = NonLinearSolverType::_auto); /// set an integration scheme for a given dof and a given solver void setIntegrationScheme(const ID & solver_id, const ID & dof_id, const IntegrationSchemeType & integration_scheme_type, IntegrationScheme::SolutionType solution_type = IntegrationScheme::_not_defined); /// set an externally instantiated integration scheme void setIntegrationScheme(const ID & solver_id, const ID & dof_id, std::unique_ptr & integration_scheme, IntegrationScheme::SolutionType solution_type = IntegrationScheme::_not_defined); /* ------------------------------------------------------------------------ */ /* SolverCallback interface */ /* ------------------------------------------------------------------------ */ public: /// Predictor interface for the callback void predictor() override; /// Corrector interface for the callback void corrector() override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Default time step solver to instantiate for this model virtual TimeStepSolverType getDefaultSolverType() const; /// Default configurations for a given time step solver virtual ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /// get access to the internal dof manager DOFManager & getDOFManager() { return *this->dof_manager; } /// get the time step of a given solver Real getTimeStep(const ID & solver_id = "") const; /// set the time step of a given solver virtual void setTimeStep(Real time_step, const ID & solver_id = ""); /// set the parameter 'param' of the solver 'solver_id' // template // void set(const ID & param, const T & value, const ID & solver_id = ""); /// get the parameter 'param' of the solver 'solver_id' // const Parameter & get(const ID & param, const ID & solver_id = "") const; /// answer to the question "does the solver exists ?" bool hasSolver(const ID & solver_id) const; /// changes the current default solver void setDefaultSolver(const ID & solver_id); /// is a default solver defined bool hasDefaultSolver() const; /// is an integration scheme set for a given solver and a given dof bool hasIntegrationScheme(const ID & solver_id, const ID & dof_id) const; TimeStepSolver & getTimeStepSolver(const ID & solver_id = ""); NonLinearSolver & getNonLinearSolver(const ID & solver_id = ""); const TimeStepSolver & getTimeStepSolver(const ID & solver_id = "") const; const NonLinearSolver & getNonLinearSolver(const ID & solver_id = "") const; private: TimeStepSolver & getSolver(const ID & solver_id); const TimeStepSolver & getSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: ModelType model_type; /// Underlying dof_manager (the brain...) std::shared_ptr dof_manager; private: ID parent_id; /// Underlying mesh Mesh & mesh; /// Default time step solver to use ID default_solver_id; }; struct ModelSolverOptions { NonLinearSolverType non_linear_solver_type; std::map integration_scheme_type; std::map solution_type; }; } // namespace akantu #endif /* AKANTU_MODEL_SOLVER_HH_ */ diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 817c09f72..2f4ffa575 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,704 +1,706 @@ /** * @file contact_mechanics_model.cc * * @author Mohit Pundir * * @date creation: Thu Feb 21 2013 * @date last modification: Wed Jul 28 2021 * * @brief Contact mechanics model * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "contact_mechanics_model.hh" #include "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #include "group_manager_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) - : Model(mesh, model_type, dof_manager, dim, id) { + : Model(mesh, model_type, dim, id) { AKANTU_DEBUG_IN(); + this->initDOFManager(dof_manager); + this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension - 1, _not_ghost, _ek_regular); 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() = default; /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // initalize the resolutions if (not this->parser.getLastParsedFile().empty()) { 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(); } } /* -------------------------------------------------------------------------- */ 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 /*unused*/) { // 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); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { // save the previous state this->savePreviousState(); contact_elements.clear(); // this needs to be resized if cohesive elements are added UInt nb_nodes = mesh.getNbNodes(); auto resize_arrays = [&](auto & internal_array) { internal_array->resize(nb_nodes); internal_array->zero(); }; 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, *tangents, *projections); // interpenetration 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.empty()) { this->computeNodalAreas(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::savePreviousState() { AKANTU_DEBUG_IN(); // saving previous natural projections (*previous_projections).copy(*projections); // saving previous tangents (*previous_tangents).copy(*tangents); // saving previous tangential traction (*previous_tangential_tractions).copy(*tangential_tractions); previous_master_elements->clear(); previous_master_elements->resize(projections->size()); previous_master_elements->set(ElementNull); for (auto & element : contact_elements) { (*previous_master_elements)[element.slave] = element.master; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::computeNodalAreas(GhostType ghost_type) { UInt nb_nodes = mesh.getNbNodes(); nodal_area->resize(nb_nodes); nodal_area->zero(); external_force->resize(nb_nodes); external_force->zero(); 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); IntegrationPoint quad_point; quad_point.ghost_type = ghost_type; auto & group = mesh.getElementGroup("contact_surface"); UInt nb_degree_of_freedom = external_force->getNbComponent(); for (auto && type : group.elementTypes(spatial_dimension - 1, ghost_type)) { const auto & element_ids = group.getElements(type, ghost_type); UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(type, ghost_type); UInt nb_elements = element_ids.size(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array dual_before_integ(nb_elements * nb_quad_points, nb_degree_of_freedom, 0.); Array quad_coords(nb_elements * nb_quad_points, spatial_dimension); const auto & normals_on_quad = fem_boundary.getNormalsOnIntegrationPoints(type, ghost_type); auto normals_begin = normals_on_quad.begin(spatial_dimension); decltype(normals_begin) normals_iter; auto quad_coords_iter = quad_coords.begin(spatial_dimension); auto dual_iter = dual_before_integ.begin(nb_degree_of_freedom); quad_point.type = type; Element subelement; subelement.type = type; subelement.ghost_type = ghost_type; for (auto el : element_ids) { subelement.element = el; const auto & element_to_subelement = mesh.getElementToSubelement(type)(el); Vector outside(spatial_dimension); mesh.getBarycenter(subelement, outside); Vector inside(spatial_dimension); if (mesh.isMeshFacets()) { mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); } else { mesh.getBarycenter(element_to_subelement[0], inside); } Vector inside_to_outside(spatial_dimension); inside_to_outside = outside - inside; normals_iter = normals_begin + el * nb_quad_points; quad_point.element = el; for (auto q : arange(nb_quad_points)) { quad_point.num_point = q; auto ddot = inside_to_outside.dot(*normals_iter); Vector normal(*normals_iter); if (ddot < 0) { normal *= -1.0; } (*dual_iter) .mul(Matrix::eye(spatial_dimension, 1), normal); ++dual_iter; ++quad_coords_iter; ++normals_iter; } } Array dual_by_shapes(nb_elements * nb_quad_points, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.computeNtb(dual_before_integ, dual_by_shapes, type, ghost_type, element_ids); Array dual_by_shapes_integ(nb_elements, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.integrate(dual_by_shapes, dual_by_shapes_integ, nb_degree_of_freedom * nb_nodes_per_element, type, ghost_type, element_ids); this->getDOFManager().assembleElementalArrayLocalArray( dual_by_shapes_integ, *external_force, type, ghost_type, 1., element_ids); } 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(); } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space(indent, 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 (const 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) const { 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_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); } } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string & /*unused*/, const std::string & /*unused*/, bool /*unused*/) { 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.get(); real_nodal_fields["normal_force"] = this->normal_force.get(); real_nodal_fields["tangential_force"] = this->tangential_force.get(); real_nodal_fields["blocked_dofs"] = this->blocked_dofs.get(); real_nodal_fields["normals"] = this->normals.get(); real_nodal_fields["tangents"] = this->tangents.get(); real_nodal_fields["gaps"] = this->gaps.get(); real_nodal_fields["areas"] = this->nodal_area.get(); real_nodal_fields["tangential_traction"] = this->tangential_tractions.get(); 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; } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldUInt(const std::string & field_name, const std::string & group_name, bool /*padding_flag*/) { std::shared_ptr field; if (field_name == "contact_state") { auto && func = std::make_unique>(); field = mesh.createNodalField(this->contact_state.get(), group_name); field = dumpers::FieldComputeProxy::createFieldCompute(field, std::move(func)); } return field; } /* -------------------------------------------------------------------------- */ 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/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index d00ef47ff..fcf8b2032 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,884 +1,886 @@ /** * @file heat_transfer_model.cc * * @author Guillaume Anciaux * @author Lucas Frerot * @author Emil Gallyamov * @author David Simon Kammer * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Fri Apr 09 2021 * * @brief Implementation of HeatTransferModel class * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_model.hh" #include "dumpable_inline_impl.hh" #include "element_synchronizer.hh" #include "fe_engine_template.hh" #include "generalized_trapezoidal.hh" #include "group_manager_inline_impl.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "parser.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ namespace akantu { namespace heat_transfer { namespace details { class ComputeRhoFunctor { public: ComputeRhoFunctor(const HeatTransferModel & model) : model(model){}; void operator()(Matrix & rho, const Element & /*unused*/) { rho.set(model.getCapacity() * model.getDensity()); } private: const HeatTransferModel & model; }; } // namespace details } // namespace heat_transfer /* -------------------------------------------------------------------------- */ HeatTransferModel::HeatTransferModel(Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager) - : Model(mesh, ModelType::_heat_transfer_model, dof_manager, dim, id), + : Model(mesh, ModelType::_heat_transfer_model, dim, id), temperature_gradient("temperature_gradient", id), temperature_on_qpoints("temperature_on_qpoints", id), conductivity_on_qpoints("conductivity_on_qpoints", id), k_gradt_on_qpoints("k_gradt_on_qpoints", id) { AKANTU_DEBUG_IN(); + this->initDOFManager(dof_manager); + conductivity = Matrix(this->spatial_dimension, this->spatial_dimension); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_htm_temperature); this->registerSynchronizer(synchronizer, SynchronizationTag::_htm_gradient_temperature); } registerFEEngineObject(id + ":fem", mesh, spatial_dimension); this->mesh.registerDumper("heat_transfer", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); this->registerParam("conductivity", conductivity, _pat_parsmod); this->registerParam("conductivity_variation", conductivity_variation, 0., _pat_parsmod); this->registerParam("temperature_reference", T_ref, 0., _pat_parsmod); this->registerParam("capacity", capacity, _pat_parsmod); this->registerParam("density", density, _pat_parsmod); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initModel() { auto & fem = this->getFEEngine(); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); temperature_on_qpoints.initialize(fem, _nb_component = 1); temperature_gradient.initialize(fem, _nb_component = spatial_dimension); conductivity_on_qpoints.initialize(fem, _nb_component = spatial_dimension * spatial_dimension); k_gradt_on_qpoints.initialize(fem, _nb_component = spatial_dimension); } /* -------------------------------------------------------------------------- */ FEEngine & HeatTransferModel::getFEEngineBoundary(const ID & name) { return aka::as_type(getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ HeatTransferModel::~HeatTransferModel() = default; /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor compute_rho(*this); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldLumped(compute_rho, "M", "temperature", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType HeatTransferModel::getMatrixType(const ID & matrix_id) const { if (matrix_id == "K" or matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleConductivityMatrix(); } else if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacity(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacityLumped(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleResidual() { AKANTU_DEBUG_IN(); this->assembleInternalHeatRate(); this->getDOFManager().assembleToResidual("temperature", *this->external_heat_rate, 1); this->getDOFManager().assembleToResidual("temperature", *this->internal_heat_rate, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::predictor() { ++temperature_release; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped() { AKANTU_DEBUG_IN(); if (!this->getDOFManager().hasLumpedMatrix("M")) { this->getDOFManager().getNewLumpedMatrix("M"); } this->getDOFManager().zeroLumpedMatrix("M"); assembleCapacityLumped(_not_ghost); assembleCapacityLumped(_ghost); need_to_reassemble_capacity_lumped = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->temperature, 1, "temperature"); this->allocNodalField(this->external_heat_rate, 1, "external_heat_rate"); this->allocNodalField(this->internal_heat_rate, 1, "internal_heat_rate"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); if (!dof_manager.hasDOFs("temperature")) { dof_manager.registerDOFs("temperature", *this->temperature, _dst_nodal); dof_manager.registerBlockedDOFs("temperature", *this->blocked_dofs); } if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(this->temperature_rate, 1, "temperature_rate"); if (!dof_manager.hasDOFsDerivatives("temperature", 1)) { dof_manager.registerDOFsDerivative("temperature", 1, *this->temperature_rate); } } } /* -------------------------------------------------------------------------- */ std::tuple HeatTransferModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } 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 HeatTransferModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_forward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_pseudo_time; options.solution_type["temperature"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_forward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; } else { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_backward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleConductivityMatrix() { AKANTU_DEBUG_IN(); this->computeConductivityOnQuadPoints(_not_ghost); if (conductivity_release[_not_ghost] == conductivity_matrix_release) { return; } AKANTU_DEBUG_ASSERT(this->getDOFManager().hasMatrix("K"), "The K matrix has not been initialized yet."); this->getDOFManager().zeroMatrix("K"); auto & fem = this->getFEEngine(); for (auto && type : mesh.elementTypes(spatial_dimension)) { auto nb_element = mesh.getNbElement(type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); fem.computeBtDB(conductivity_on_qpoints(type), *bt_d_b, 2, type); /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto K_e = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_e"); fem.integrate(*bt_d_b, *K_e, nb_nodes_per_element * nb_nodes_per_element, type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "temperature", *K_e, type, _not_ghost, _symmetric); } conductivity_matrix_release = conductivity_release[_not_ghost]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeConductivityOnQuadPoints(GhostType ghost_type) { // if already computed once check if need to compute if (not initial_conductivity[ghost_type]) { // if temperature did not change, conductivity will not vary if (temperature_release == conductivity_release[ghost_type]) { return; } // if conductivity_variation is 0 no need to recompute if (conductivity_variation == 0.) { return; } } for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & temperature_interpolated = temperature_on_qpoints(type, ghost_type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type, ghost_type); auto & cond = conductivity_on_qpoints(type, ghost_type); for (auto && tuple : zip(make_view(cond, spatial_dimension, spatial_dimension), temperature_interpolated)) { auto & C = std::get<0>(tuple); auto & T = std::get<1>(tuple); C = conductivity; Matrix variation(spatial_dimension, spatial_dimension, conductivity_variation * (T - T_ref)); // @TODO: Guillaume are you sure ? why due you compute variation then ? C += conductivity_variation; } } conductivity_release[ghost_type] = temperature_release; initial_conductivity[ghost_type] = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeKgradT(GhostType ghost_type) { computeConductivityOnQuadPoints(ghost_type); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & gradient = temperature_gradient(type, ghost_type); this->getFEEngine().gradientOnIntegrationPoints(*temperature, gradient, 1, type, ghost_type); for (auto && values : zip(make_view(conductivity_on_qpoints(type, ghost_type), spatial_dimension, spatial_dimension), make_view(gradient, spatial_dimension), make_view(k_gradt_on_qpoints(type, ghost_type), spatial_dimension))) { const auto & C = std::get<0>(values); const auto & BT = std::get<1>(values); auto & k_BT = std::get<2>(values); k_BT.mul(C, BT); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleInternalHeatRate() { AKANTU_DEBUG_IN(); this->internal_heat_rate->zero(); this->synchronize(SynchronizationTag::_htm_temperature); auto & fem = this->getFEEngine(); for (auto ghost_type : ghost_types) { // compute k \grad T computeKgradT(ghost_type); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto & k_gradt_on_qpoints_vect = k_gradt_on_qpoints(type, ghost_type); UInt nb_quad_points = k_gradt_on_qpoints_vect.size(); Array bt_k_gT(nb_quad_points, nb_nodes_per_element); fem.computeBtD(k_gradt_on_qpoints_vect, bt_k_gT, type, ghost_type); UInt nb_elements = mesh.getNbElement(type, ghost_type); Array int_bt_k_gT(nb_elements, nb_nodes_per_element); fem.integrate(bt_k_gT, int_bt_k_gT, nb_nodes_per_element, type, ghost_type); this->getDOFManager().assembleElementalArrayLocalArray( int_bt_k_gT, *this->internal_heat_rate, type, ghost_type, -1); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real el_size; Real min_el_size = std::numeric_limits::max(); Real conductivitymax = conductivity(0, 0); // get the biggest parameter from k11 until k33// for (UInt i = 0; i < spatial_dimension; i++) { for (UInt j = 0; j < spatial_dimension; j++) { conductivitymax = std::max(conductivity(i, j), conductivitymax); } } for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array coord(0, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, _not_ghost); auto el_coord = coord.begin(spatial_dimension, nb_nodes_per_element); UInt nb_element = mesh.getNbElement(type); for (UInt el = 0; el < nb_element; ++el, ++el_coord) { el_size = getFEEngine().getElementInradius(*el_coord, type); min_el_size = std::min(min_el_size, el_size); } AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size << " and the max conductivity is : " << conductivitymax); } Real min_dt = 2. * min_el_size * min_el_size / 4. * density * capacity / conductivitymax; mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); this->mesh.getDumper("heat_transfer").setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::readMaterials() { auto sect = this->getParserSection(); if (not std::get<1>(sect)) { const auto & section = std::get<0>(sect); this->parseSection(section); } conductivity_on_qpoints.set(conductivity); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); readMaterials(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacity() { AKANTU_DEBUG_IN(); auto ghost_type = _not_ghost; this->getDOFManager().zeroMatrix("M"); auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor rho_functor(*this); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldMatrix(rho_functor, "M", "temperature", this->getDOFManager(), type, ghost_type); } need_to_reassemble_capacity = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeRho(Array & rho, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); FEEngine & fem = this->getFEEngine(); UInt nb_element = mesh.getNbElement(type, ghost_type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); rho.resize(nb_element * nb_quadrature_points); rho.set(this->capacity); // Real * rho_1_val = rho.storage(); // /// compute @f$ rho @f$ for each nodes of each element // for (UInt el = 0; el < nb_element; ++el) { // for (UInt n = 0; n < nb_quadrature_points; ++n) { // *rho_1_val++ = this->capacity; // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::computeThermalEnergyByNode() { AKANTU_DEBUG_IN(); Real ethermal = 0.; for (auto && pair : enumerate(make_view( *internal_heat_rate, internal_heat_rate->getNbComponent()))) { auto n = std::get<0>(pair); auto & heat_rate = std::get<1>(pair); Real heat = 0.; bool is_local_node = mesh.isLocalOrMasterNode(n); bool count_node = is_local_node; for (UInt i = 0; i < heat_rate.size(); ++i) { if (count_node) { heat += heat_rate[i] * time_step; } } ethermal += heat; } mesh.getCommunicator().allReduce(ethermal, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ethermal; } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::getThermalEnergy( iterator Eth, Array::const_iterator T_it, const Array::const_iterator & T_end) const { for (; T_it != T_end; ++T_it, ++Eth) { *Eth = capacity * density * *T_it; } } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy(ElementType type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Vector Eth_on_quarature_points(nb_quadrature_points); auto T_it = this->temperature_on_qpoints(type).begin(); T_it += index * nb_quadrature_points; auto T_end = T_it + nb_quadrature_points; getThermalEnergy(Eth_on_quarature_points.storage(), T_it, T_end); return getFEEngine().integrate(Eth_on_quarature_points, type, index); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy() { Real Eth = 0; auto & fem = getFEEngine(); for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_element = mesh.getNbElement(type, _not_ghost); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, _not_ghost); Array Eth_per_quad(nb_element * nb_quadrature_points, 1); auto & temperature_interpolated = temperature_on_qpoints(type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type); auto T_it = temperature_interpolated.begin(); auto T_end = temperature_interpolated.end(); getThermalEnergy(Eth_per_quad.begin(), T_it, T_end); Eth += fem.integrate(Eth_per_quad, type); } return Eth; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id) { AKANTU_DEBUG_IN(); Real energy = 0; if (id == "thermal") { energy = getThermalEnergy(); } // reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); Real energy = 0.; if (id == "thermal") { energy = getThermalEnergy(type, index); } AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ std::shared_ptr HeatTransferModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs.get(); auto field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr HeatTransferModel::createNodalFieldReal( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { if (field_name == "capacity_lumped") { AKANTU_EXCEPTION( "Capacity lumped is a nodal field now stored in the DOF manager." "Therefore it cannot be used by a dumper anymore"); } std::map *> real_nodal_fields; real_nodal_fields["temperature"] = temperature.get(); real_nodal_fields["temperature_rate"] = temperature_rate.get(); real_nodal_fields["external_heat_rate"] = external_heat_rate.get(); real_nodal_fields["internal_heat_rate"] = internal_heat_rate.get(); real_nodal_fields["increment"] = increment.get(); std::shared_ptr field = mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr HeatTransferModel::createElementalField( const std::string & field_name, const std::string & group_name, bool /*padding_flag*/, UInt /*spatial_dimension*/, ElementKind element_kind) { std::shared_ptr field; if (field_name == "partitions") { field = mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); } else if (field_name == "temperature_gradient") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(temperature_gradient); field = mesh.createElementalField( temperature_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "conductivity") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(conductivity_on_qpoints); field = mesh.createElementalField( conductivity_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbData(const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = indexes.size(); switch (tag) { case SynchronizationTag::_htm_temperature: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case SynchronizationTag::_htm_temperature: { buffer << (*temperature)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case SynchronizationTag::_htm_temperature: { buffer >> (*temperature)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; Array::const_iterator it = elements.begin(); Array::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case SynchronizationTag::_htm_temperature: { size += nb_nodes_per_element * sizeof(Real); // temperature break; } case SynchronizationTag::_htm_gradient_temperature: { // temperature gradient size += getNbIntegrationPoints(elements) * spatial_dimension * sizeof(Real); size += nb_nodes_per_element * sizeof(Real); // nodal temperatures break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { switch (tag) { case SynchronizationTag::_htm_temperature: { packNodalDataHelper(*temperature, buffer, elements, mesh); break; } case SynchronizationTag::_htm_gradient_temperature: { packElementalDataHelper(temperature_gradient, buffer, elements, true, getFEEngine()); packNodalDataHelper(*temperature, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { switch (tag) { case SynchronizationTag::_htm_temperature: { unpackNodalDataHelper(*temperature, buffer, elements, mesh); break; } case SynchronizationTag::_htm_gradient_temperature: { unpackElementalDataHelper(temperature_gradient, buffer, elements, true, getFEEngine()); unpackNodalDataHelper(*temperature, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/model.cc b/src/model/model.cc index 71f166251..9f84d962b 100644 --- a/src/model/model.cc +++ b/src/model/model.cc @@ -1,362 +1,352 @@ /** * @file model.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Mon Oct 03 2011 * @date last modification: Fri Apr 09 2021 * * @brief implementation of model common parts * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "model.hh" #include "communicator.hh" #include "data_accessor.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ namespace akantu { -/* -------------------------------------------------------------------------- */ -Model::Model(Mesh & mesh, const ModelType & type, - std::shared_ptr dof_manager, UInt dim, const ID & id) - : ModelSolver(mesh, type, id, std::move(dof_manager)), mesh(mesh), - spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension() - : dim), - parser(getStaticParser()) { - this->mesh.registerEventHandler(*this, _ehp_model); -} - /* -------------------------------------------------------------------------- */ Model::Model(Mesh & mesh, const ModelType & type, UInt dim, const ID & id) : ModelSolver(mesh, type, id), mesh(mesh), spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension() : dim), parser(getStaticParser()) { this->mesh.registerEventHandler(*this, _ehp_model); } /* -------------------------------------------------------------------------- */ Model::~Model() = default; /* -------------------------------------------------------------------------- */ void Model::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); method = options.analysis_method; if (!this->hasDefaultSolver()) { this->initNewSolver(this->method); } initModel(); initFEEngineBoundary(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Model::initNewSolver(const AnalysisMethod & method) { ID solver_name; TimeStepSolverType tss_type; std::tie(solver_name, tss_type) = this->getDefaultSolverID(method); if (not this->hasSolver(solver_name)) { ModelSolverOptions options = this->getDefaultSolverOptions(tss_type); this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type); for (auto && is_type : options.integration_scheme_type) { if (!this->hasIntegrationScheme(solver_name, is_type.first)) { this->setIntegrationScheme(solver_name, is_type.first, is_type.second, options.solution_type[is_type.first]); } } } this->method = method; this->setDefaultSolver(solver_name); } /* -------------------------------------------------------------------------- */ void Model::initFEEngineBoundary() { try { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } catch (debug::Exception & /*e*/) { } } /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name, const std::string & dumper_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup() { for (auto & group : mesh.iterateElementGroups()) { group.dump(); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory) { for (auto & group : mesh.iterateElementGroups()) { group.setDirectory(directory); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Model::setGroupBaseName(const std::string & basename, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setBaseName(basename); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Model::getGroupDumper(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); return group.getDumper(); } /* -------------------------------------------------------------------------- */ // DUMPER stuff /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & field_id, std::shared_ptr field, DumperIOHelper & dumper) { dumper.registerField(field_id, std::move(field)); } /* -------------------------------------------------------------------------- */ void Model::addDumpField(const std::string & field_id) { this->addDumpFieldToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldVector(const std::string & field_id) { this->addDumpFieldVectorToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldTensor(const std::string & field_id) { this->addDumpFieldTensorToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::setBaseName(const std::string & field_id) { mesh.setBaseName(field_id); } /* -------------------------------------------------------------------------- */ void Model::setBaseNameToDumper(const std::string & dumper_name, const std::string & basename) { mesh.setBaseNameToDumper(dumper_name, basename); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", dumper_default_element_kind, false); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id, group_name, dumper_default_element_kind, false); } /* -------------------------------------------------------------------------- */ void Model::removeDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void Model::removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.removeDumpFieldFromDumper(dumper_name, field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", dumper_default_element_kind, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name, dumper_default_element_kind, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", dumper_default_element_kind, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag) { this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name, this->spatial_dimension, element_kind, padding_flag); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, ElementKind element_kind, bool padding_flag) { std::shared_ptr field; if (!field) { field = this->createNodalFieldReal(field_id, group_name, padding_flag); } if (!field) { field = this->createNodalFieldUInt(field_id, group_name, padding_flag); } if (!field) { field = this->createNodalFieldBool(field_id, group_name, padding_flag); } if (!field) { field = this->createElementalField(field_id, group_name, padding_flag, spatial_dimension, element_kind); } if (!field) { field = this->mesh.createFieldFromAttachedData(field_id, group_name, element_kind); } if (!field) { field = this->mesh.createFieldFromAttachedData(field_id, group_name, element_kind); } #ifndef AKANTU_NDEBUG if (!field) { AKANTU_DEBUG_WARNING("No field could be found based on name: " << field_id); } #endif if (field) { DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name, group_name); this->addDumpGroupFieldToDumper(field_id, field, dumper); } } /* -------------------------------------------------------------------------- */ void Model::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void Model::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void Model::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void Model::dump() { auto default_dumper = mesh.getDefaultDumperName(); this->dump(default_dumper); } /* -------------------------------------------------------------------------- */ void Model::dump(UInt step) { auto default_dumper = mesh.getDefaultDumperName(); this->dump(default_dumper, step); } /* -------------------------------------------------------------------------- */ void Model::dump(Real time, UInt step) { auto default_dumper = mesh.getDefaultDumperName(); this->dump(default_dumper, time, step); } /* -------------------------------------------------------------------------- */ void Model::setDirectory(const std::string & directory) { mesh.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Model::setDirectoryToDumper(const std::string & dumper_name, const std::string & directory) { mesh.setDirectoryToDumper(dumper_name, directory); } /* -------------------------------------------------------------------------- */ void Model::setTextModeToDumper() { mesh.setTextModeToDumper(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/model.hh b/src/model/model.hh index d0040fe98..b19da1d7a 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,390 +1,384 @@ /** * @file model.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Apr 09 2021 * * @brief Interface of a model * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_named_argument.hh" #include "fe_engine.hh" #include "mesh.hh" #include "model_options.hh" #include "model_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MODEL_HH_ #define AKANTU_MODEL_HH_ namespace akantu { class SynchronizerRegistry; class Parser; class DumperIOHelper; class DOFManager; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { class Model : public ModelSolver, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Normal constructor where the DOFManager is created internally Model(Mesh & mesh, const ModelType & type, UInt dim = _all_dimensions, const ID & id = "model"); - /// Model constructor the the dof manager is created externally, for example - /// in a ModelCoupler - Model(Mesh & mesh, const ModelType & type, - std::shared_ptr dof_manager, UInt dim = _all_dimensions, - const ID & id = "model"); - ~Model() override; using FEEngineMap = std::map>; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: virtual void initFullImpl(const ModelOptions & options); public: template std::enable_if_t::value> initFull(pack &&... _pack) { switch (this->model_type) { #ifdef AKANTU_SOLID_MECHANICS case ModelType::_solid_mechanics_model: this->initFullImpl(SolidMechanicsModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_COHESIVE_ELEMENT case ModelType::_solid_mechanics_model_cohesive: this->initFullImpl(SolidMechanicsModelCohesiveOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_HEAT_TRANSFER case ModelType::_heat_transfer_model: this->initFullImpl(HeatTransferModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_PHASE_FIELD case ModelType::_phase_field_model: this->initFullImpl(PhaseFieldModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_EMBEDDED case ModelType::_embedded_model: this->initFullImpl(EmbeddedInterfaceModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_CONTACT_MECHANICS case ModelType::_contact_mechanics_model: this->initFullImpl(ContactMechanicsModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_MODEL_COUPLERS case ModelType::_coupler_solid_contact: this->initFullImpl(CouplerSolidContactOptions{ use_named_args, std::forward(_pack)...}); break; case ModelType::_coupler_solid_cohesive_contact: this->initFullImpl(CouplerSolidCohesiveContactOptions{ use_named_args, std::forward(_pack)...}); break; #endif default: this->initFullImpl(ModelOptions{use_named_args, std::forward(_pack)...}); } } template std::enable_if_t::value> initFull(pack &&... _pack) { this->initFullImpl(std::forward(_pack)...); } /// initialize a new solver if needed void initNewSolver(const AnalysisMethod & method); protected: /// get some default values for derived classes virtual std::tuple getDefaultSolverID(const AnalysisMethod & method) = 0; virtual void initModel(){}; virtual void initFEEngineBoundary(); /// function to print the containt of the class void printself(std::ostream & /*stream*/, int /*indent*/ = 0) const override{}; public: /* ------------------------------------------------------------------------ */ /* Access to the dumpable interface of the boundaries */ /* ------------------------------------------------------------------------ */ /// Dump the data for a given group void dumpGroup(const std::string & group_name); void dumpGroup(const std::string & group_name, const std::string & dumper_name); /// Dump the data for all boundaries void dumpGroup(); /// Set the directory for a given group void setGroupDirectory(const std::string & directory, const std::string & group_name); /// Set the directory for all boundaries void setGroupDirectory(const std::string & directory); /// Set the base name for a given group void setGroupBaseName(const std::string & basename, const std::string & group_name); /// Get the internal dumper of a given group DumperIOHelper & getGroupDumper(const std::string & group_name); /* ------------------------------------------------------------------------ */ /* Function for non local capabilities */ /* ------------------------------------------------------------------------ */ virtual void updateDataForNonLocalCriterion(__attribute__((unused)) ElementTypeMapReal & criterion) { AKANTU_TO_IMPLEMENT(); } protected: template void allocNodalField(std::unique_ptr> & array, UInt nb_component, const ID & name) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get id of model AKANTU_GET_MACRO(ID, id, const ID &) /// get the number of surfaces AKANTU_GET_MACRO(Mesh, mesh, Mesh &) /// synchronize the boundary in case of parallel run virtual void synchronizeBoundaries(){}; /// return the fem object associated with a provided name inline FEEngine & getFEEngine(const ID & name = "") const; /// return the fem boundary object associated with a provided name virtual FEEngine & getFEEngineBoundary(const ID & name = ""); inline bool hasFEEngineBoundary(const ID & name = ""); /// register a fem object associated with name template inline void registerFEEngineObject(const std::string & name, Mesh & mesh, UInt spatial_dimension); /// unregister a fem object associated with name inline void unRegisterFEEngineObject(const std::string & name); /// return the synchronizer registry SynchronizerRegistry & getSynchronizerRegistry(); /// return the fem object associated with a provided name template inline FEEngineClass & getFEEngineClass(std::string name = "") const; /// return the fem boundary object associated with a provided name template inline FEEngineClass & getFEEngineClassBoundary(std::string name = ""); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /* ------------------------------------------------------------------------ */ /* Pack and unpack hexlper functions */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbIntegrationPoints(const Array & elements, const ID & fem_id = ID()) const; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ void setTextModeToDumper(); virtual void addDumpGroupFieldToDumper(const std::string & field_id, std::shared_ptr field, DumperIOHelper & dumper); virtual void addDumpField(const std::string & field_id); virtual void addDumpFieldVector(const std::string & field_id); virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensor(const std::string & field_id); virtual void setBaseName(const std::string & field_id); virtual void setBaseNameToDumper(const std::string & dumper_name, const std::string & basename); virtual void addDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, ElementKind element_kind, bool padding_flag); virtual void removeDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual std::shared_ptr createNodalFieldReal(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createNodalFieldUInt(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createNodalFieldBool(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createElementalField( const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/, UInt /*spatial_dimension*/, ElementKind /*kind*/) { return nullptr; } void setDirectory(const std::string & directory); void setDirectoryToDumper(const std::string & dumper_name, const std::string & directory); /* ------------------------------------------------------------------------ */ virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); /* ------------------------------------------------------------------------ */ virtual void dump(); virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: friend std::ostream & operator<<(std::ostream & /*stream*/, const Model & /*_this*/); ID id; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// Mesh Mesh & mesh; /// Spatial dimension of the problem UInt spatial_dimension; /// the main fem object present in all models FEEngineMap fems; /// the fem object present in all models for boundaries FEEngineMap fems_boundary; /// default fem object std::string default_fem; /// parser to the pointer to use Parser & parser; /// default ElementKind for dumper ElementKind dumper_default_element_kind{_ek_regular}; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Model & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "model_inline_impl.hh" #endif /* AKANTU_MODEL_HH_ */ diff --git a/src/model/model_couplers/coupler_solid_cohesive_contact.cc b/src/model/model_couplers/coupler_solid_cohesive_contact.cc index 3a2099ad0..c5e6ced73 100644 --- a/src/model/model_couplers/coupler_solid_cohesive_contact.cc +++ b/src/model/model_couplers/coupler_solid_cohesive_contact.cc @@ -1,76 +1,77 @@ /** * @file coupler_solid_cohesive_contact.cc * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Mon Jan 21 2019 * @date last modification: Wed Jun 23 2021 * * @brief class for coupling of solid mechanics cohesive and conatct mechanics * model * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "coupler_solid_cohesive_contact.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template <> CouplerSolidContactTemplate:: CouplerSolidContactTemplate(Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager) - : Model(mesh, ModelType::_coupler_solid_cohesive_contact, dof_manager, dim, - id) { + : Model(mesh, ModelType::_coupler_solid_cohesive_contact, dim, id) { + this->initDOFManager(dof_manager); + this->mesh.registerDumper("coupler_solid_cohesive_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_cohesive_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); this->registerDataAccessor(*this); solid = std::make_unique( mesh, Model::spatial_dimension, "solid_mechanics_model_cohesive", this->dof_manager); contact = std::make_unique(mesh.getMeshFacets(), Model::spatial_dimension, "contact_mechanics_model"); } /* -------------------------------------------------------------------------- */ template <> void CouplerSolidContactTemplate::initFullImpl( const ModelOptions & options) { Model::initFullImpl(options); const auto & cscc_options = aka::as_type(options); solid->initFull(_analysis_method = cscc_options.analysis_method, _is_extrinsic = cscc_options.is_extrinsic); contact->initFull(_analysis_method = cscc_options.analysis_method); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_contact.cc b/src/model/model_couplers/coupler_solid_contact.cc index 5f111e9e2..9f6b0961f 100644 --- a/src/model/model_couplers/coupler_solid_contact.cc +++ b/src/model/model_couplers/coupler_solid_contact.cc @@ -1,69 +1,72 @@ /** * @file coupler_solid_contact.cc * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Mon Jan 21 2019 * @date last modification: Wed Jun 23 2021 * * @brief class for coupling of solid mechanics and conatct mechanics * model * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "coupler_solid_contact.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template <> CouplerSolidContactTemplate::CouplerSolidContactTemplate( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager) - : Model(mesh, ModelType::_coupler_solid_contact, dof_manager, dim, id) { + : Model(mesh, ModelType::_coupler_solid_contact, dim, id) { + + this->initDOFManager(dof_manager); + this->mesh.registerDumper("coupler_solid_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); this->registerDataAccessor(*this); solid = std::make_unique(mesh, Model::spatial_dimension, "solid_mechanics_model", this->dof_manager); contact = std::make_unique( mesh, Model::spatial_dimension, "contact_mechanics_model", this->dof_manager); } /* -------------------------------------------------------------------------- */ template <> void CouplerSolidContactTemplate::initFullImpl( const ModelOptions & options) { Model::initFullImpl(options); solid->initFull(_analysis_method = this->method); contact->initFull(_analysis_method = this->method); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 228a7ab7d..7136890bb 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1258 +1,1259 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Mauro Corrado * @author Aurelia Isabel Cuba Ramos * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Clement Roux * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Apr 09 2021 * * @brief Implementation of the SolidMechanicsModel class * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "solid_mechanics_model_tmpl.hh" #include "communicator.hh" #include "element_synchronizer.hh" #include "sparse_matrix.hh" #include "synchronizer_registry.hh" #include "dumpable_inline_impl.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model * @param model_type this is an internal parameter for inheritance purposes */ SolidMechanicsModel::SolidMechanicsModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) - : Model(mesh, model_type, std::move(dof_manager), dim, id), - material_index("material index", id), + : Model(mesh, model_type, dim, id), material_index("material index", id), material_local_numbering("material local numbering", id) { AKANTU_DEBUG_IN(); + this->initDOFManager(dof_manager); + this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); this->mesh.registerDumper("solid_mechanics_model", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); material_selector = std::make_shared(material_index); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_material_id); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_mass); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_stress); this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); this->mesh.getDumper().setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ /* Initialization */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param options * \parblock * contains the different options to initialize the model * \li \c analysis_method specify the type of solver to use * \endparblock */ void SolidMechanicsModel::initFullImpl(const ModelOptions & options) { material_index.initialize(mesh, _element_kind = _ek_not_defined, _default_value = UInt(-1), _with_nb_element = true); material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true); Model::initFullImpl(options); // initialize the materials if (not this->parser.getLastParsedFile().empty()) { this->instantiateMaterials(); this->initMaterials(); } this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModel::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::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ std::tuple SolidMechanicsModel::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); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { auto & dof_manager = this->getDOFManager(); /* ------------------------------------------------------------------------ */ // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->previous_displacement, spatial_dimension, "previous_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->blocked_dofs, spatial_dimension, "blocked_dofs"); this->allocNodalField(this->current_position, spatial_dimension, "current_position"); // initialize the current positions this->current_position->copy(this->mesh.getNodes()); /* ------------------------------------------------------------------------ */ if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); dof_manager.registerDOFsIncrement("displacement", *this->displacement_increment); dof_manager.registerDOFsPrevious("displacement", *this->previous_displacement); } /* ------------------------------------------------------------------------ */ // for dynamic if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(this->velocity, spatial_dimension, "velocity"); this->allocNodalField(this->acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->assembleInternalForces(); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) const { // \TODO check the materials to know what is the correct answer if (matrix_id == "C") { return _mt_not_defined; } if (matrix_id == "K") { auto matrix_type = _unsymmetric; for (auto & material : materials) { matrix_type = std::max(matrix_type, material->getMatrixType(matrix_id)); } } return _symmetric; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { this->assembleMass(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { this->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::beforeSolveStep() { for (auto & material : materials) { material->beforeSolveStep(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::afterSolveStep(bool converged) { for (auto & material : materials) { material->afterSolveStep(converged); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::predictor() { ++displacement_release; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::corrector() { ++displacement_release; } /* -------------------------------------------------------------------------- */ /** * This function computes the internal forces as \f$F_{int} = \int_{\Omega} N * \sigma d\Omega@\f$ */ void SolidMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->zero(); // compute the stresses of local elements AKANTU_DEBUG_INFO("Compute local stresses"); for (auto & material : materials) { material->computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ if (this->non_local_manager) { this->non_local_manager->computeAllNonLocalStresses(); } // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_smm_stress); // assemble the forces due to local stresses AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & material : materials) { material->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(SynchronizationTag::_smm_stress); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for (auto & material : materials) { material->assembleInternalForces(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix(bool need_to_reassemble) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); if (not this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", this->getMatrixType("K")); } // Check if materials need to recompute the matrix for (auto & material : materials) { need_to_reassemble |= material->hasMatrixChanged("K"); } if (need_to_reassemble) { this->getDOFManager().getMatrix("K").zero(); // call compute stiffness matrix on each local elements for (auto & material : materials) { material->assembleStiffnessMatrix(_not_ghost); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { if (this->current_position_release == this->displacement_release) { return; } this->current_position->copy(this->mesh.getNodes()); auto cpos_it = this->current_position->begin(Model::spatial_dimension); auto cpos_end = this->current_position->end(Model::spatial_dimension); auto disp_it = this->displacement->begin(Model::spatial_dimension); for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) { *cpos_it += *disp_it; } this->current_position_release = this->displacement_release; } /* -------------------------------------------------------------------------- */ const Array & SolidMechanicsModel::getCurrentPosition() { this->updateCurrentPosition(); return *this->current_position; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateDataForNonLocalCriterion( ElementTypeMapReal & criterion) { const ID field_name = criterion.getName(); for (auto & material : materials) { if (!material->isInternal(field_name, _ek_regular)) { continue; } for (auto ghost_type : ghost_types) { material->flattenInternal(field_name, criterion, ghost_type, _ek_regular); } } } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(GhostType ghost_type) { AKANTU_DEBUG_IN(); Real min_dt = std::numeric_limits::max(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { elem.type = type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); auto mat_indexes = material_index(type, ghost_type).begin(); auto mat_loc_num = material_local_numbering(type, ghost_type).begin(); Array X(0, nb_nodes_per_element * Model::spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, type, _not_ghost); auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, type); Real el_c = this->materials[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); if (this->getDOFManager().hasLumpedMatrix("M")) { this->assembleLumpedMatrix("M"); auto m_it = this->mass->begin(Model::spatial_dimension); auto m_end = this->mass->end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) { const auto & v = *v_it; const auto & m = *m_it; Real mv2 = 0.; auto is_local_node = mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); auto count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (m(i) > std::numeric_limits::epsilon()) { mv2 += v(i) * v(i) * m(i); } } } ekin += mv2; } } else if (this->getDOFManager().hasMatrix("M")) { this->assembleMatrix("M"); Array Mv(nb_nodes, Model::spatial_dimension); this->getDOFManager().assembleMatMulVectToArray("displacement", "M", *this->velocity, Mv); for (auto && data : zip(arange(nb_nodes), make_view(Mv, spatial_dimension), make_view(*this->velocity, spatial_dimension))) { ekin += std::get<2>(data).dot(std::get<1>(data)) * static_cast(mesh.isLocalOrMasterNode(std::get<0>(data))); } } else { AKANTU_ERROR("No function called to assemble the mass matrix."); } mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(ElementType type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array vel_on_quad(nb_quadrature_points, Model::spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, Model::spatial_dimension, type, _not_ghost, filter_element); auto vit = vel_on_quad.begin(Model::spatial_dimension); auto vend = vel_on_quad.end(Model::spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5 * getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); auto ext_force_it = external_force->begin(Model::spatial_dimension); auto int_force_it = internal_force->begin(Model::spatial_dimension); auto boun_it = blocked_dofs->begin(Model::spatial_dimension); decltype(ext_force_it) incr_or_velo_it; if (this->method == _static) { incr_or_velo_it = this->displacement_increment->begin(Model::spatial_dimension); } else { incr_or_velo_it = this->velocity->begin(Model::spatial_dimension); } Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) { const auto & int_force = *int_force_it; const auto & ext_force = *ext_force_it; const auto & boun = *boun_it; const auto & incr_or_velo = *incr_or_velo_it; bool is_local_node = this->mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (boun(i)) { work -= int_force(i) * incr_or_velo(i); } else { work += ext_force(i) * incr_or_velo(i); } } } } mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum); if (this->method != _static) { work *= this->getTimeStep(); } AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } if (energy_id == "external work") { return getExternalWork(); } Real energy = 0.; for (auto & material : materials) { energy += material->getEnergy(energy_id); } /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const ID & energy_id, const ID & group_id) { auto && group = mesh.getElementGroup(group_id); auto energy = 0.; for (auto && type : group.elementTypes()) { for (auto el : group.getElementsIterable(type)) { energy += getEnergy(energy_id, el); } } /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->material_index.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true, _default_value = UInt(-1)); this->material_local_numbering.initialize( mesh, _element_kind = _ek_not_defined, _with_nb_element = true, _default_value = UInt(-1)); ElementTypeMapArray filter("new_element_filter", this->getID()); for (const auto & elem : element_list) { if (mesh.getSpatialDimension(elem.type) != spatial_dimension) { continue; } if (!filter.exists(elem.type, elem.ghost_type)) { filter.alloc(0, 1, elem.type, elem.ghost_type); } filter(elem.type, elem.ghost_type).push_back(elem.element); } // this fails in parallel if the event is sent on facet between constructor // and initFull \todo: to debug... this->assignMaterialToElements(&filter); for (auto & material : materials) { material->onElementsAdded(element_list, event); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { for (auto & material : materials) { material->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if (displacement) { displacement->resize(nb_nodes, 0.); ++displacement_release; } if (mass) { mass->resize(nb_nodes, 0.); } if (velocity) { velocity->resize(nb_nodes, 0.); } if (acceleration) { acceleration->resize(nb_nodes, 0.); } if (external_force) { external_force->resize(nb_nodes, 0.); } if (internal_force) { internal_force->resize(nb_nodes, 0.); } if (blocked_dofs) { blocked_dofs->resize(nb_nodes, false); } if (current_position) { current_position->resize(nb_nodes, 0.); } if (previous_displacement) { previous_displacement->resize(nb_nodes, 0.); } if (displacement_increment) { displacement_increment->resize(nb_nodes, 0.); } for (auto & material : materials) { material->onNodesAdded(nodes_list, event); } need_to_reassemble_lumped_mass = true; need_to_reassemble_mass = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(const Array & /*element_list*/, const Array & new_numbering, const RemovedNodesEvent & /*event*/) { if (displacement) { mesh.removeNodesFromArray(*displacement, new_numbering); ++displacement_release; } if (mass) { mesh.removeNodesFromArray(*mass, new_numbering); } if (velocity) { mesh.removeNodesFromArray(*velocity, new_numbering); } if (acceleration) { mesh.removeNodesFromArray(*acceleration, new_numbering); } if (internal_force) { mesh.removeNodesFromArray(*internal_force, new_numbering); } if (external_force) { mesh.removeNodesFromArray(*external_force, new_numbering); } if (blocked_dofs) { mesh.removeNodesFromArray(*blocked_dofs, new_numbering); } // if (increment_acceleration) // mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if (displacement_increment) { mesh.removeNodesFromArray(*displacement_increment, new_numbering); } if (previous_displacement) { mesh.removeNodesFromArray(*previous_displacement, new_numbering); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Solid 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 << " ]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); if (velocity) { velocity->printself(stream, indent + 2); } if (acceleration) { acceleration->printself(stream, indent + 2); } if (mass) { mass->printself(stream, indent + 2); } external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + materials [" << std::endl; for (const auto & material : materials) { material->printself(stream, indent + 2); } stream << space << " ]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeNonLocal() { this->non_local_manager->synchronize(*this, SynchronizationTag::_material_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods( GhostType ghost_type) { for (auto & mat : materials) { MaterialNonLocalInterface * mat_non_local; if ((mat_non_local = dynamic_cast(mat.get())) == nullptr) { continue; } ElementTypeMapArray quadrature_points_coordinates( "quadrature_points_coordinates_tmp_nl", this->id); quadrature_points_coordinates.initialize(this->getFEEngine(), _nb_component = spatial_dimension, _ghost_type = ghost_type); for (const auto & type : quadrature_points_coordinates.elementTypes( Model::spatial_dimension, ghost_type)) { this->getFEEngine().computeIntegrationPointsCoordinates( quadrature_points_coordinates(type, ghost_type), type, ghost_type); } mat_non_local->initMaterialNonLocal(); mat_non_local->insertIntegrationPointsInNeighborhoods( ghost_type, quadrature_points_coordinates); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeNonLocalStresses(GhostType ghost_type) { for (auto & mat : materials) { if (not aka::is_of_type(*mat)) { continue; } auto & mat_non_local = dynamic_cast(*mat); mat_non_local.computeNonLocalStresses(ghost_type); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateLocalInternal( ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) { const ID field_name = internal_flat.getName(); for (auto & material : materials) { if (material->isInternal(field_name, kind)) { material->flattenInternal(field_name, internal_flat, ghost_type, kind); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateNonLocalInternal( ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) { const ID field_name = internal_flat.getName(); for (auto & mat : materials) { if (not aka::is_of_type(*mat)) { continue; } auto & mat_non_local = dynamic_cast(*mat); mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) { return getFEEngineClassBoundary(name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::splitElementByMaterial( const Array & elements, std::vector> & elements_per_mat) const { for (const auto & el : elements) { Element mat_el = el; mat_el.element = this->material_local_numbering(el); elements_per_mat[this->material_index(el)].push_back(mat_el); } } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::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); } switch (tag) { case SynchronizationTag::_material_id: { size += elements.size() * sizeof(UInt); break; } case SynchronizationTag::_smm_mass: { size += nb_nodes_per_element * sizeof(Real) * Model::spatial_dimension; // mass vector break; } case SynchronizationTag::_smm_for_gradu: { size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real); // displacement break; } case SynchronizationTag::_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * Model::spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } case SynchronizationTag::_for_dump: { // displacement, velocity, acceleration, residual, force size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5; break; } default: { } } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { size += mat.getNbData(elements, tag); }); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_material_id: { packElementalDataHelper(material_index, buffer, elements, false, getFEEngine()); break; } case SynchronizationTag::_smm_mass: { packNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { packNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_for_dump: { packNodalDataHelper(*displacement, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*acceleration, buffer, elements, mesh); packNodalDataHelper(*internal_force, buffer, elements, mesh); packNodalDataHelper(*external_force, buffer, elements, mesh); break; } case SynchronizationTag::_smm_boundary: { packNodalDataHelper(*external_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.packData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_material_id: { for (auto && element : elements) { UInt recv_mat_index; buffer >> recv_mat_index; UInt & mat_index = material_index(element); if (mat_index != UInt(-1)) { continue; } // add ghosts element to the correct material mat_index = recv_mat_index; UInt index = materials[mat_index]->addElement(element); material_local_numbering(element) = index; } break; } case SynchronizationTag::_smm_mass: { unpackNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_for_dump: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*acceleration, buffer, elements, mesh); unpackNodalDataHelper(*internal_force, buffer, elements, mesh); unpackNodalDataHelper(*external_force, buffer, elements, mesh); break; } case SynchronizationTag::_smm_boundary: { unpackNodalDataHelper(*external_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.unpackData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::getNbData(const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch (tag) { case SynchronizationTag::_smm_uv: { size += sizeof(Real) * Model::spatial_dimension * 2; break; } case SynchronizationTag::_smm_res: /* FALLTHRU */ case SynchronizationTag::_smm_mass: { size += sizeof(Real) * Model::spatial_dimension; break; } case SynchronizationTag::_for_dump: { size += sizeof(Real) * Model::spatial_dimension * 5; break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { packDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { packDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); packDOFDataHelper(*acceleration, buffer, dofs); packDOFDataHelper(*internal_force, buffer, dofs); packDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { unpackDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { unpackDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); unpackDOFDataHelper(*acceleration, buffer, dofs); unpackDOFDataHelper(*internal_force, buffer, dofs); unpackDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } } // namespace akantu