diff --git a/examples/python/solver_callback/solver_callback.py b/examples/python/solver_callback/solver_callback.py index 347c069c8..3517ece2a 100644 --- a/examples/python/solver_callback/solver_callback.py +++ b/examples/python/solver_callback/solver_callback.py @@ -1,178 +1,172 @@ #!/usr/bin/env python3 """ solver_callback.py: solver_callback overload example""" __author__ = "Nicolas Richart" __credits__ = [ "Guillaume Anciaux ", "Nicolas Richart ", ] __copyright__ = "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __license__ = "LGPLv3" import numpy as np import akantu as aka -class SolverCallback(aka.SolverCallback): +class SolverCallback(aka.InterceptSolverCallback): def __init__(self, model): super().__init__(model.getDOFManager()) self.model = model mesh = model.getMesh() left = mesh.getElementGroup("Left").getNodeGroup().getNodes() right = mesh.getElementGroup("Right").getNodeGroup().getNodes() position = mesh.getNodes() self.pair = [] for node_l in left: node_l = int(node_l) for node_r in right: node_r = int(node_r) if abs(position[node_r, 1] - position[node_l, 1]) < 1e-6: self.pair.append([node_l, node_r]) blocked_dofs = model.getBlockedDOFs() self.periodic_K_modif = aka.TermsToAssemble("displacement", "displacement") matrix_type = self.model.getMatrixType("K") for p in self.pair: #blocked_dofs[p[1]] = True self.periodic_K_modif(p[0]*2, p[0]*2, 1) self.periodic_K_modif(p[0]*2, p[1]*2, -1) if matrix_type == aka._unsymmetric: self.periodic_K_modif(p[1]*2, p[0]*2, -1) self.periodic_K_modif(p[1]*2, p[1]*2, 1) self.first = True self.k_release = -1 - def getMatrixType(self, matrix_id): - return self.model.getMatrixType(matrix_id) - def assembleMatrix(self, matrix_id): self.model.assembleMatrix(matrix_id) if matrix_id == "K": release = self.model.getDOFManager().getMatrix("K").getRelease() if release == self.k_release: return if self.first: self.model.getDOFManager().getMatrix("K").saveMatrix("K0.mtx") self.model.getDOFManager().assemblePreassembledMatrix( "K", self.periodic_K_modif) if self.first: self.model.getDOFManager().getMatrix("K").saveMatrix("K1.mtx") self.k_release = self.model.getDOFManager().getMatrix("K").getRelease() self.first = False - def assembleLumpedMatrix(self, matrix_id): - self.model.assembleLumpedMatrix(matrix_id) - def assembleResidual(self): displacement = self.model.getDisplacement() force = np.zeros(displacement.shape) for p in self.pair: force[p[0], 0] += displacement[p[0], 0] - displacement[p[1], 0] force[p[1], 0] += displacement[p[1], 0] - displacement[p[0], 0] self.model.getDOFManager().assembleToResidual('displacement', force, -1.); self.model.assembleResidual(); # ----------------------------------------------------------------------------- def main(): spatial_dimension = 2 mesh_file = 'bar.msh' max_steps = 250 time_step = 1e-3 aka.parseInput('material.dat') # ------------------------------------------------------------------------- # Initialization # ------------------------------------------------------------------------- mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) model = aka.SolidMechanicsModel(mesh) model.initFull(_analysis_method=aka._implicit_dynamic) model.setBaseName("solver_callback") model.addDumpFieldVector("displacement") model.addDumpFieldVector("acceleration") model.addDumpFieldVector("velocity") model.addDumpFieldVector("internal_force") model.addDumpFieldVector("external_force") model.addDumpField("strain") model.addDumpField("stress") model.addDumpField("blocked_dofs") # ------------------------------------------------------------------------- # boundary conditions # ------------------------------------------------------------------------- model.applyBC(aka.FixedValue(0, aka._y), "YBlocked") # ------------------------------------------------------------------------- # initial conditions # ------------------------------------------------------------------------- displacement = model.getDisplacement() velocity = model.getVelocity() nb_nodes = mesh.getNbNodes() position = mesh.getNodes() L = 1 # pulse_width A = 0.01 v = np.sqrt(model.getMaterial(0).getReal('E') / model.getMaterial(0).getReal('rho')) k = 0.1 * 2 * np.pi * 3 / L t = 0. velocity[:, 0] = k * v * A * np.sin(k * ((position[:, 0] - 5.) - v * t)) displacement[:, 0] = A * np.cos(k * ((position[:, 0] - 5.) - v * t)) # ------------------------------------------------------------------------- # timestep value computation # ------------------------------------------------------------------------- time_factor = 0.8 stable_time_step = model.getStableTimeStep() * time_factor print("Stable Time Step = {0}".format(stable_time_step)) print("Required Time Step = {0}".format(time_step)) time_step = stable_time_step * time_factor model.setTimeStep(time_step) solver_callback = SolverCallback(model) solver = model.getNonLinearSolver() solver.set("max_iterations", 100) solver.set("threshold", 1e-7) # ------------------------------------------------------------------------- # loop for evolution of motion dynamics # ------------------------------------------------------------------------- print("step,step * time_step,epot,ekin,epot + ekin") for step in range(0, max_steps + 1): model.solveStep(solver_callback) #model.solveStep() if step % 10 == 0: model.dump() epot = model.getEnergy('potential') ekin = model.getEnergy('kinetic') # output energy calculation to screen print("{0},{1},{2},{3},{4}".format(step, step * time_step, epot, ekin, (epot + ekin))) return # ----------------------------------------------------------------------------- if __name__ == "__main__": main() diff --git a/python/py_dof_manager.cc b/python/py_dof_manager.cc index 627e2a742..5c8263521 100644 --- a/python/py_dof_manager.cc +++ b/python/py_dof_manager.cc @@ -1,159 +1,166 @@ /* * 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 "py_dof_manager.hh" #include "py_aka_array.hh" #include "py_akantu_pybind11_compatibility.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { - class PySolverCallback : public SolverCallback { + template + class PySolverCallback : public Parent { public: - using SolverCallback::SolverCallback; + using Parent::Parent; /// get the type of matrix needed - MatrixType getMatrixType(const ID & matrix_id) override { + MatrixType getMatrixType(const ID & matrix_id) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(MatrixType, SolverCallback, getMatrixType, matrix_id); } /// callback to assemble a Matrix void assembleMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleMatrix, matrix_id); } /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleLumpedMatrix, matrix_id); } /// callback to assemble the residual (rhs) void assembleResidual() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleResidual); } /* ---------------------------------------------------------------------- */ /* Dynamic simulations part */ /* ---------------------------------------------------------------------- */ /// callback for the predictor (in case of dynamic simulation) void predictor() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, predictor); } /// callback for the corrector (in case of dynamic simulation) void corrector() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, corrector); } void beforeSolveStep() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, beforeSolveStep); } void afterSolveStep(bool converged) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, afterSolveStep, converged); } }; } // namespace /* -------------------------------------------------------------------------- */ void register_dof_manager(py::module & mod) { py::class_(mod, "DOFManager") .def("getMatrix", &DOFManager::getMatrix, py::return_value_policy::reference) .def( "getNewMatrix", [](DOFManager & self, const std::string & name, const std::string & matrix_to_copy_id) -> decltype(auto) { return self.getNewMatrix(name, matrix_to_copy_id); }, py::return_value_policy::reference) .def( "getResidual", [](DOFManager & self) -> decltype(auto) { return self.getResidual(); }, py::return_value_policy::reference) .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs) .def( "hasMatrix", [](DOFManager & self, const ID & name) -> bool { return self.hasMatrix(name); }, py::arg("name")) .def("assembleToResidual", &DOFManager::assembleToResidual, py::arg("dof_id"), py::arg("array_to_assemble"), py::arg("scale_factor") = 1.) .def("assembleToLumpedMatrix", &DOFManager::assembleToLumpedMatrix, py::arg("dof_id"), py::arg("array_to_assemble"), py::arg("lumped_mtx"), py::arg("scale_factor") = 1.) .def("assemblePreassembledMatrix", &DOFManager::assemblePreassembledMatrix, py::arg("matrix_id"), py::arg("terms")); py::class_(mod, "NonLinearSolver") .def( "set", [](NonLinearSolver & self, const std::string & id, const Real & val) { if (id == "max_iterations") { self.set(id, int(val)); } else { self.set(id, val); } }) .def("set", [](NonLinearSolver & self, const std::string & id, const SolveConvergenceCriteria & val) { self.set(id, val); }); - py::class_(mod, "SolverCallback") + py::class_>(mod, + "SolverCallback") .def(py::init_alias()) .def("getMatrixType", &SolverCallback::getMatrixType) .def("assembleMatrix", &SolverCallback::assembleMatrix) .def("assembleLumpedMatrix", &SolverCallback::assembleLumpedMatrix) .def("assembleResidual", [](SolverCallback & self) { self.assembleResidual(); }) .def("predictor", &SolverCallback::predictor) .def("corrector", &SolverCallback::corrector) .def("beforeSolveStep", &SolverCallback::beforeSolveStep) .def("afterSolveStep", &SolverCallback::afterSolveStep) .def_property_readonly("dof_manager", &SolverCallback::getSCDOFManager, py::return_value_policy::reference); + + py::class_>( + mod, "InterceptSolverCallback") + .def(py::init_alias()); } } // namespace akantu diff --git a/src/model/common/solver_callback.hh b/src/model/common/solver_callback.hh index 7d329b4a4..bb42d4082 100644 --- a/src/model/common/solver_callback.hh +++ b/src/model/common/solver_callback.hh @@ -1,112 +1,168 @@ /** * @file solver_callback.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Nov 27 2019 * * @brief Class defining the interface for non_linear_solver callbacks * * * @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" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_CALLBACK_HH_ #define AKANTU_SOLVER_CALLBACK_HH_ namespace akantu { class DOFManager; } namespace akantu { class SolverCallback { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: explicit SolverCallback(DOFManager & dof_manager); explicit SolverCallback(); /* ------------------------------------------------------------------------ */ virtual ~SolverCallback(); protected: void setDOFManager(DOFManager & dof_manager); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get the type of matrix needed - virtual MatrixType getMatrixType(const ID &) = 0; + virtual MatrixType getMatrixType(const ID &) const = 0; /// callback to assemble a Matrix virtual void assembleMatrix(const ID &) = 0; /// callback to assemble a lumped Matrix virtual void assembleLumpedMatrix(const ID &) = 0; /// callback to assemble the residual (rhs) virtual void assembleResidual() = 0; /// callback to assemble the rhs parts, (e.g. internal_forces + /// external_forces) virtual void assembleResidual(const ID & /*residual_part*/) {} /* ------------------------------------------------------------------------ */ /* Dynamic simulations part */ /* ------------------------------------------------------------------------ */ /// callback for the predictor (in case of dynamic simulation) virtual void predictor() {} /// callback for the corrector (in case of dynamic simulation) virtual void corrector() {} /// tells if the residual can be computed in separated parts - virtual bool canSplitResidual() { return false; } + virtual bool canSplitResidual() const { return false; } /* ------------------------------------------------------------------------ */ /* management callbacks */ /* ------------------------------------------------------------------------ */ virtual void beforeSolveStep() {} virtual void afterSolveStep(bool /*converged*/ = true) {} DOFManager & getSCDOFManager() { return *sc_dof_manager; } protected: /// DOFManager prefixed to avoid collision in multiple inheritance cases DOFManager * sc_dof_manager{nullptr}; }; namespace debug { class SolverCallbackResidualPartUnknown : public Exception { public: SolverCallbackResidualPartUnknown(const ID & residual_part) : Exception(residual_part + " is not known here.") {} }; } // namespace debug +/* -------------------------------------------------------------------------- */ +class InterceptSolverCallback : public SolverCallback { +public: + InterceptSolverCallback(SolverCallback & solver_callback) + : solver_callback(solver_callback) {} + + /// get the type of matrix needed + MatrixType getMatrixType(const ID & matrix_id) const override { + return solver_callback.getMatrixType(matrix_id); + } + + /// callback to assemble a Matrix + void assembleMatrix(const ID & matrix_id) override { + solver_callback.assembleMatrix(matrix_id); + } + + /// callback to assemble a lumped Matrix + void assembleLumpedMatrix(const ID & matrix_id) override { + solver_callback.assembleLumpedMatrix(matrix_id); + } + + /// callback to assemble the residual (rhs) + void assembleResidual() override { solver_callback.assembleResidual(); } + + /// callback to assemble the rhs parts, (e.g. internal_forces + + /// external_forces) + void assembleResidual(const ID & residual_part) override { + solver_callback.assembleResidual(residual_part); + } + + /* ------------------------------------------------------------------------ */ + /* Dynamic simulations part */ + /* ------------------------------------------------------------------------ */ + /// callback for the predictor (in case of dynamic simulation) + void predictor() override { solver_callback.predictor(); } + + /// callback for the corrector (in case of dynamic simulation) + void corrector() override { solver_callback.corrector(); } + + /// tells if the residual can be computed in separated parts + bool canSplitResidual() const override { + return solver_callback.canSplitResidual(); + } + + /* ------------------------------------------------------------------------ */ + /* management callbacks */ + /* ------------------------------------------------------------------------ */ + void beforeSolveStep() override { solver_callback.beforeSolveStep(); } + void afterSolveStep(bool converged = true) override { + solver_callback.afterSolveStep(converged); + } + +protected: + SolverCallback & solver_callback; +}; + } // namespace akantu #endif /* AKANTU_SOLVER_CALLBACK_HH_ */ diff --git a/src/model/common/time_step_solvers/time_step_solver.hh b/src/model/common/time_step_solvers/time_step_solver.hh index 8108b9616..f4456e5cb 100644 --- a/src/model/common/time_step_solvers/time_step_solver.hh +++ b/src/model/common/time_step_solvers/time_step_solver.hh @@ -1,167 +1,167 @@ /** * @file time_step_solver.hh * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Sep 08 2020 * * @brief This corresponding to the time step evolution solver * * * @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_array.hh" #include "integration_scheme.hh" #include "parameter_registry.hh" #include "solver_callback.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_TIME_STEP_SOLVER_HH_ #define AKANTU_TIME_STEP_SOLVER_HH_ namespace akantu { class DOFManager; class NonLinearSolver; } // namespace akantu namespace akantu { class TimeStepSolver : public ParameterRegistry, public SolverCallback { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: TimeStepSolver(DOFManager & dof_manager, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, SolverCallback & solver_callback, const ID & id); ~TimeStepSolver() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// solves on step virtual void solveStep(SolverCallback & solver_callback) = 0; /// register an integration scheme for a given dof void setIntegrationScheme(const ID & dof_id, const IntegrationSchemeType & type, IntegrationScheme::SolutionType solution_type = IntegrationScheme::_not_defined); /// register an integration scheme for a given dof void setIntegrationScheme(const ID & dof_id, std::unique_ptr & scheme, IntegrationScheme::SolutionType solution_type = IntegrationScheme::_not_defined); protected: /// register an integration scheme for a given dof virtual std::unique_ptr getIntegrationSchemeInternal(const ID & dof_id, const IntegrationSchemeType & type, IntegrationScheme::SolutionType solution_type = IntegrationScheme::_not_defined) = 0; virtual void setIntegrationSchemeInternal(const ID & dof_id, std::unique_ptr & scheme, IntegrationScheme::SolutionType solution_type = IntegrationScheme::_not_defined) = 0; public: /// replies if a integration scheme has been set virtual bool hasIntegrationScheme(const ID & dof_id) const = 0; /* ------------------------------------------------------------------------ */ /* Solver Callback interface */ /* ------------------------------------------------------------------------ */ public: /// implementation of the SolverCallback::getMatrixType() - MatrixType getMatrixType(const ID & /*unused*/) final { + MatrixType getMatrixType(const ID & /*unused*/) const final { return _mt_not_defined; } /// implementation of the SolverCallback::predictor() void predictor() override; /// implementation of the SolverCallback::corrector() void corrector() override; /// implementation of the SolverCallback::assembleJacobian() void assembleMatrix(const ID & matrix_id) override; /// implementation of the SolverCallback::assembleJacobian() void assembleLumpedMatrix(const ID & matrix_id) override; /// implementation of the SolverCallback::assembleResidual() void assembleResidual() override; /// implementation of the SolverCallback::assembleResidual() void assembleResidual(const ID & residual_part) override; void beforeSolveStep() override; void afterSolveStep(bool converged = true) override; - bool canSplitResidual() override { + bool canSplitResidual() const override { return solver_callback->canSplitResidual(); } /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(TimeStep, time_step, Real); AKANTU_SET_MACRO(TimeStep, time_step, Real); AKANTU_GET_MACRO(NonLinearSolver, non_linear_solver, const NonLinearSolver &); AKANTU_GET_MACRO_NOT_CONST(NonLinearSolver, non_linear_solver, NonLinearSolver &); protected: MatrixType getCommonMatrixType(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: ID id; /// Underlying dof manager containing the dof to treat DOFManager & _dof_manager; /// Type of solver TimeStepSolverType type; /// The time step for this solver Real time_step; /// Temporary storage for solver callback SolverCallback * solver_callback; /// NonLinearSolver used by this tome step solver NonLinearSolver & non_linear_solver; /// List of required matrices std::map needed_matrices; /// specifies if the solvers gives to full solution or just the increment of /// solution bool is_solution_increment{true}; }; } // namespace akantu #endif /* AKANTU_TIME_STEP_SOLVER_HH_ */ diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index a37223369..3b9b48f11 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,713 +1,713 @@ /** * @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" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #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) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, id + ":contact_detector"); registerFEEngineObject("ContactFacetsFEEngine", mesh, Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::~ContactMechanicsModel() = 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); // intepenetration value must be positive for contact mechanics // model to work by default the gap value from detector is negative std::for_each((*gaps).begin(), (*gaps).end(), [](Real & gap) { gap *= -1.; }); if (!contact_elements.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) { +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); } } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ 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; } #endif /* -------------------------------------------------------------------------- */ 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/contact_mechanics/contact_mechanics_model.hh b/src/model/contact_mechanics/contact_mechanics_model.hh index b57d4472f..cecc2e64b 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,374 +1,374 @@ /** * @file contact_mechanics_model.hh * * @author Mohit Pundir * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Jun 24 2021 * * @brief Model of Contact Mechanics * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "contact_detector.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_MECHANICS_MODEL_HH__ #define __AKANTU_CONTACT_MECHANICS_MODEL_HH__ namespace akantu { class Resolution; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class ContactMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: ContactMechanicsModel( Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "contact_mechanics_model", std::shared_ptr dof_manager = nullptr, ModelType model_type = ModelType::_contact_mechanics_model); ~ContactMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// allocate all vectors void initSolver(TimeStepSolverType /*unused*/, NonLinearSolverType /*unused*/) override; /// initialize all internal arrays for resolutions void initResolutions(); /// initialize the modelType void initModel() override; /// call back for the solver, computes the force residual void assembleResidual() override; /// get the type of matrix needed - MatrixType getMatrixType(const ID & matrix_id) override; + MatrixType getMatrixType(const ID & matrix_id) const override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Contact Detection */ /* ------------------------------------------------------------------------ */ public: void search(); void computeNodalAreas(GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Contact Resolution */ /* ------------------------------------------------------------------------ */ public: /// register an empty contact resolution of a given type Resolution & registerNewResolution(const ID & res_name, const ID & res_type, const ID & opt_param); protected: /// register a resolution in the dynamic database Resolution & registerNewResolution(const ParserSection & res_section); /// read the resolution files to instantiate all the resolutions void instantiateResolutions(); /// save the parameters from previous state void savePreviousState(); /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: #if defined(AKANTU_USE_IOHELPER) std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldUInt(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; #endif /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: /// contact detection class friend class ContactDetector; /// contact resolution class friend class Resolution; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the ContactMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the ContactMechanicsModel::increment vector \warn only consistent /// if ContactMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the ContactMechanics::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the ContactMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the ContactMechanics::normal_force vector (normal forces) AKANTU_GET_MACRO(NormalForce, *normal_force, Array &); /// get the ContactMechanics::tangential_force vector (friction forces) AKANTU_GET_MACRO(TangentialForce, *tangential_force, Array &); /// get the ContactMechanics::traction vector (friction traction) AKANTU_GET_MACRO(TangentialTractions, *tangential_tractions, Array &); /// get the ContactMechanics::previous_tangential_tractions vector AKANTU_GET_MACRO(PreviousTangentialTractions, *previous_tangential_tractions, Array &); /// get the ContactMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the ContactMechanics::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the ContactMechanics::gaps (contact gaps) AKANTU_GET_MACRO(Gaps, *gaps, Array &); /// get the ContactMechanics::normals (normals on slave nodes) AKANTU_GET_MACRO(Normals, *normals, Array &); /// get the ContactMechanics::tangents (tangents on slave nodes) AKANTU_GET_MACRO(Tangents, *tangents, Array &); /// get the ContactMechanics::previous_tangents (tangents on slave nodes) AKANTU_GET_MACRO(PreviousTangents, *previous_tangents, Array &); /// get the ContactMechanics::areas (nodal areas) AKANTU_GET_MACRO(NodalArea, *nodal_area, Array &); /// get the ContactMechanics::previous_projections (previous_projections) AKANTU_GET_MACRO(PreviousProjections, *previous_projections, Array &); /// get the ContactMechanics::projections (projections) AKANTU_GET_MACRO(Projections, *projections, Array &); /// get the ContactMechanics::contact_state vector (no_contact/stick/slip /// state) AKANTU_GET_MACRO(ContactState, *contact_state, Array &); /// get the ContactMechanics::previous_master_elements AKANTU_GET_MACRO(PreviousMasterElements, *previous_master_elements, Array &); /// get contact detector AKANTU_GET_MACRO_NOT_CONST(ContactDetector, *detector, ContactDetector &); /// get the contact elements inline const Array & getContactElements() const { return contact_elements; } /// get the current positions of the nodes inline Array & getPositions() { return detector->getPositions(); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// tells if the resolutions are instantiated bool are_resolutions_instantiated; /// displacements array std::unique_ptr> displacement; /// increment of displacement std::unique_ptr> displacement_increment; /// contact forces array std::unique_ptr> internal_force; /// external forces array std::unique_ptr> external_force; /// normal force array std::unique_ptr> normal_force; /// friction force array std::unique_ptr> tangential_force; /// friction traction array std::unique_ptr> tangential_tractions; /// previous friction traction array std::unique_ptr> previous_tangential_tractions; /// boundary vector std::unique_ptr> blocked_dofs; /// array to store gap between slave and master std::unique_ptr> gaps; /// array to store normals from master to slave std::unique_ptr> normals; /// array to store tangents on the master element std::unique_ptr> tangents; /// array to store previous tangents on the master element std::unique_ptr> previous_tangents; /// array to store nodal areas std::unique_ptr> nodal_area; /// array to store stick/slip state : std::unique_ptr> contact_state; /// array to store previous projections in covariant basis std::unique_ptr> previous_projections; // array to store projections in covariant basis std::unique_ptr> projections; /// contact detection std::unique_ptr detector; /// list of contact resolutions std::vector> resolutions; /// mapping between resolution name and resolution internal id std::map resolutions_names_to_id; /// array to store contact elements Array contact_elements; /// array to store previous master elements std::unique_ptr> previous_master_elements; }; } // namespace akantu /* ------------------------------------------------------------------------ */ /* inline functions */ /* ------------------------------------------------------------------------ */ #include "parser.hh" #include "resolution.hh" /* ------------------------------------------------------------------------ */ #endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ diff --git a/src/model/model_couplers/cohesive_contact_solvercallback.cc b/src/model/model_couplers/cohesive_contact_solvercallback.cc index 445fc895c..8d933d710 100644 --- a/src/model/model_couplers/cohesive_contact_solvercallback.cc +++ b/src/model/model_couplers/cohesive_contact_solvercallback.cc @@ -1,141 +1,141 @@ /** * @file cohesive_contact_solvercallback.cc * * @author Mohit Pundir * * @date creation: Sun Jun 06 2021 * @date last modification: Sun Jun 06 2021 * * @brief class for coupling of solid mechanics cohesive and contact mechanics * model via solvercallback * * * @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 "cohesive_contact_solvercallback.hh" /* -------------------------------------------------------------------------- */ namespace akantu { CohesiveContactSolverCallback::CohesiveContactSolverCallback( SolidMechanicsModelCohesive & solid, ContactMechanicsModel & contact, AnalysisMethod & method) : solid(solid), contact(contact), method(method) {} void CohesiveContactSolverCallback::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { solid.assembleStiffnessMatrix(); switch (method) { case _static: case _implicit_dynamic: { contact.assembleStiffnessMatrix(); break; } default: break; } } else if (matrix_id == "M") { solid.assembleMass(); } } /* -------------------------------------------------------------------------- */ void CohesiveContactSolverCallback::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid.assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CohesiveContactSolverCallback::assembleResidual() { // computes the internal forces switch (method) { case _explicit_lumped_mass: { auto & current_positions = contact.getContactDetector().getPositions(); current_positions.copy(solid.getCurrentPosition()); contact.search(); break; } default: break; } solid.assembleInternalForces(); contact.assembleInternalForces(); auto & internal_force = solid.getInternalForce(); auto & external_force = solid.getExternalForce(); auto & contact_force = contact.getInternalForce(); solid.getDOFManager().assembleToResidual("displacement", external_force, 1); solid.getDOFManager().assembleToResidual("displacement", internal_force, 1); solid.getDOFManager().assembleToResidual("displacement", contact_force, 1); } /* -------------------------------------------------------------------------- */ void CohesiveContactSolverCallback::predictor() { auto & solid_model_solver = aka::as_type(solid); solid_model_solver.predictor(); } /* -------------------------------------------------------------------------- */ void CohesiveContactSolverCallback::beforeSolveStep() { auto & solid_model_solver = aka::as_type(solid); solid_model_solver.beforeSolveStep(); } /* -------------------------------------------------------------------------- */ void CohesiveContactSolverCallback::afterSolveStep(bool converged) { auto & solid_model_solver = aka::as_type(solid); solid_model_solver.afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ void CohesiveContactSolverCallback::corrector() { auto & solid_model_solver = aka::as_type(solid); solid_model_solver.corrector(); switch (method) { case _static: case _implicit_dynamic: { auto & current_positions = contact.getContactDetector().getPositions(); current_positions.copy(solid.getCurrentPosition()); contact.search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ MatrixType -CohesiveContactSolverCallback::getMatrixType(const ID & /*matrix_id*/) { +CohesiveContactSolverCallback::getMatrixType(const ID & /*matrix_id*/) const { return _symmetric; } } // namespace akantu diff --git a/src/model/model_couplers/cohesive_contact_solvercallback.hh b/src/model/model_couplers/cohesive_contact_solvercallback.hh index 078e56e19..12ab5dbad 100644 --- a/src/model/model_couplers/cohesive_contact_solvercallback.hh +++ b/src/model/model_couplers/cohesive_contact_solvercallback.hh @@ -1,90 +1,90 @@ /** * @file cohesive_contact_solvercallback.hh * * @author Mohit Pundir * * @date creation: Sun Jun 06 2021 * @date last modification: Wed Jul 28 2021 * * @brief class for coupling of solid mechanics and conatct mechanics * model via solvercallback * * * @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 "contact_mechanics_model.hh" #include "mesh_iterators.hh" #include "non_linear_solver.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COHESIVE_CONTACT_SOLVERCALLBACK_HH__ #define __AKANTU_COHESIVE_CONTACT_SOLVERCALLBACK_HH__ namespace akantu { class CohesiveContactSolverCallback : public SolverCallback { public: CohesiveContactSolverCallback(SolidMechanicsModelCohesive & /*solid*/, ContactMechanicsModel & /*contact*/, AnalysisMethod & /*method*/); public: /// implementation of SolverCallback::assembleMatrix void assembleMatrix(const ID & /*matrix_id*/) override; /// implementation of SolverCallback::assembleResidual void assembleResidual() override; /// implementation of SolverCallback::assembleLumpedMatrix void assembleLumpedMatrix(const ID & /*matrix_id*/) override; /// implementation of SolverCallback::getMatrixType - MatrixType getMatrixType(const ID & /*unused*/) override; + MatrixType getMatrixType(const ID & /*unused*/) const override; /// implementation of SolverCallback::predictor void predictor() override; /// implementation of SolverCallback::corrector void corrector() override; /// implementation of SolverCallback::beforeSolveStep void beforeSolveStep() override; /// implementation of SolverCallback::afterSolveStep void afterSolveStep(bool converged = true) override; private: /// model for the solid mechanics part of the coupling SolidMechanicsModelCohesive & solid; /// model for the contact resoluion of the coupling ContactMechanicsModel & contact; /// Method of resolution for the coupling solver AnalysisMethod & method; }; } // namespace akantu #endif diff --git a/src/model/model_couplers/coupler_solid_contact.hh b/src/model/model_couplers/coupler_solid_contact.hh index a6a99e77e..2744a5562 100644 --- a/src/model/model_couplers/coupler_solid_contact.hh +++ b/src/model/model_couplers/coupler_solid_contact.hh @@ -1,282 +1,282 @@ /** * @file coupler_solid_contact.hh * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Sat Jun 26 2021 * * @brief class for coupling of solid mechanics and conatct mechanics * model in explicit * * * @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 "contact_mechanics_model.hh" #include "solid_mechanics_model.hh" #if defined(AKANTU_COHESIVE_ELEMENT) #include "solid_mechanics_model_cohesive.hh" #endif /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COUPLER_SOLID_CONTACT_HH__ #define __AKANTU_COUPLER_SOLID_CONTACT_HH__ /* ------------------------------------------------------------------------ */ /* Coupling : Solid Mechanics / Contact Mechanics */ /* ------------------------------------------------------------------------ */ namespace akantu { /* -------------------------------------------------------------------------- */ template class CouplerSolidContactTemplate : public Model, public DataAccessor, public DataAccessor { static_assert( std::is_base_of::value, "SolidMechanicsModelType should be derived from SolidMechanicsModel"); /* ------------------------------------------------------------------------ */ /* Constructor/Destructor */ /* ------------------------------------------------------------------------ */ public: CouplerSolidContactTemplate( Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "coupler_solid_contact", std::shared_ptr dof_manager = nullptr, ModelType model_type = std::is_same::value ? ModelType::_coupler_solid_cohesive_contact : ModelType::_coupler_solid_contact); ~CouplerSolidContactTemplate() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); #if defined(AKANTU_COHESIVE_ELEMENT) template ::value> * = nullptr> UInt checkCohesiveStress() { return solid->checkCohesiveStress(); } #endif template inline void applyBC(const FunctorType & func) { solid->applyBC(func); } template inline void applyBC(const FunctorType & func, const std::string & group_name) { solid->applyBC(func, group_name); } template inline void applyBC(const FunctorType & func, const ElementGroup & element_group) { solid->applyBC(func, element_group); } protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; - bool canSplitResidual() override { return true; } + bool canSplitResidual() const override { return true; } /// get the type of matrix needed - MatrixType getMatrixType(const ID & matrix_id) override; + MatrixType getMatrixType(const ID & matrix_id) const override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// callback for the model to instantiate the matricess when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; /* ------------------------------------------------------------------------ */ /* Mass matrix for solid mechanics model */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass; } /* ------------------------------------------------------------------------ */ public: // DataAccessor UInt getNbData(const Array & /*elements*/, const SynchronizationTag & /*tag*/) const override { return 0; } void packData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) const override {} void unpackData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) override {} // DataAccessor nodes UInt getNbData(const Array & /*nodes*/, const SynchronizationTag & /*tag*/) const override { return 0; } void packData(CommunicationBuffer & /*buffer*/, const Array & /*nodes*/, const SynchronizationTag & /*tag*/) const override {} void unpackData(CommunicationBuffer & /*buffer*/, const Array & /*nodes*/, const SynchronizationTag & /*tag*/) override {} /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the solid mechanics model #if defined(AKANTU_COHESIVE_ELEMENT) template ::value> * = nullptr> SolidMechanicsModelCohesive & getSolidMechanicsModelCohesive() { return *solid; } #endif template ::value> * = nullptr> SolidMechanicsModelType & getSolidMechanicsModel() { return *solid; } /// get the contact mechanics model AKANTU_GET_MACRO(ContactMechanicsModel, *contact, ContactMechanicsModel &) /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: #if defined(AKANTU_USE_IOHELPER) std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldUInt(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; #endif void dump(const std::string & dumper_name) override; void dump(const std::string & dumper_name, UInt step) override; void dump(const std::string & dumper_name, Real time, UInt step) override; void dump() override; void dump(UInt step) override; void dump(Real time, UInt step) override; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ private: /// solid mechanics model std::unique_ptr solid; /// contact mechanics model std::unique_ptr contact; UInt step; }; using CouplerSolidContact = CouplerSolidContactTemplate; } // namespace akantu #include "coupler_solid_contact_tmpl.hh" #endif /* __COUPLER_SOLID_CONTACT_HH__ */ diff --git a/src/model/model_couplers/coupler_solid_contact_tmpl.hh b/src/model/model_couplers/coupler_solid_contact_tmpl.hh index 76d14ae35..496c37caf 100644 --- a/src/model/model_couplers/coupler_solid_contact_tmpl.hh +++ b/src/model/model_couplers/coupler_solid_contact_tmpl.hh @@ -1,401 +1,401 @@ /** * @file coupler_solid_contact_tmpl.hh * * @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" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template CouplerSolidContactTemplate< SolidMechanicsModelType>::~CouplerSolidContactTemplate() = default; /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); auto & contact_model_solver = aka::as_type(*contact); contact_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); } /* -------------------------------------------------------------------------- */ template std::tuple CouplerSolidContactTemplate::getDefaultSolverID( const AnalysisMethod & method) { return solid->getDefaultSolverID(method); } /* -------------------------------------------------------------------------- */ template TimeStepSolverType CouplerSolidContactTemplate::getDefaultSolverType() const { return solid->getDefaultSolverType(); } /* -------------------------------------------------------------------------- */ template ModelSolverOptions CouplerSolidContactTemplate::getDefaultSolverOptions( const TimeStepSolverType & type) const { return solid->getDefaultSolverOptions(type); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::assembleResidual() { // computes the internal forces switch (method) { case _explicit_lumped_mass: { auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(solid->getCurrentPosition()); contact->search(); break; } default: break; } this->assembleInternalForces(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", internal_force, 1); this->getDOFManager().assembleToResidual("displacement", contact_force, 1); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::assembleResidual( const ID & residual_part) { AKANTU_DEBUG_IN(); // contact->assembleInternalForces(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", contact_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("displacement", internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::predictor() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.predictor(); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::corrector() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.corrector(); switch (method) { case _static: case _implicit_dynamic: { auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(solid->getCurrentPosition()); contact->search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ template MatrixType CouplerSolidContactTemplate::getMatrixType( - const ID & matrix_id) { + const ID & matrix_id) const { if (matrix_id == "K") { return _symmetric; } if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::assembleMatrix( const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::assembleLumpedMatrix( const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::beforeSolveStep() { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.beforeSolveStep(); auto & contact_solver_callback = aka::as_type(*contact); contact_solver_callback.beforeSolveStep(); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::afterSolveStep( bool converged) { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.afterSolveStep(converged); auto & contact_solver_callback = aka::as_type(*contact); contact_solver_callback.afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate< SolidMechanicsModelType>::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); contact->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate< SolidMechanicsModelType>::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(true); switch (method) { case _static: case _implicit_dynamic: { contact->assembleStiffnessMatrix(); break; } default: break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate< SolidMechanicsModelType>::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::assembleMassLumped( GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::assembleMass( GhostType ghost_type) { solid->assembleMass(ghost_type); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ template std::shared_ptr CouplerSolidContactTemplate::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) { std::shared_ptr field; field = contact->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); if (not field) { field = solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); } return field; } /* -------------------------------------------------------------------------- */ template std::shared_ptr CouplerSolidContactTemplate::createNodalFieldReal( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; field = contact->createNodalFieldReal(field_name, group_name, padding_flag); if (not field) { field = solid->createNodalFieldReal(field_name, group_name, padding_flag); } return field; } /* -------------------------------------------------------------------------- */ template std::shared_ptr CouplerSolidContactTemplate::createNodalFieldUInt( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; field = contact->createNodalFieldUInt(field_name, group_name, padding_flag); if (not field) { field = solid->createNodalFieldUInt(field_name, group_name, padding_flag); } return field; } /* -------------------------------------------------------------------------- */ template std::shared_ptr CouplerSolidContactTemplate::createNodalFieldBool( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; field = contact->createNodalFieldBool(field_name, group_name, padding_flag); if (not field) { field = solid->createNodalFieldBool(field_name, group_name, padding_flag); } return field; } #endif /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::dump( const std::string & dumper_name) { solid->onDump(); Model::dump(dumper_name); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::dump( const std::string & dumper_name, UInt step) { solid->onDump(); Model::dump(dumper_name, step); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::dump( const std::string & dumper_name, Real time, UInt step) { solid->onDump(); Model::dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::dump() { solid->onDump(); Model::dump(); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::dump(UInt step) { solid->onDump(); Model::dump(step); } /* -------------------------------------------------------------------------- */ template void CouplerSolidContactTemplate::dump(Real time, UInt step) { solid->onDump(); Model::dump(time, step); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index a36cb1397..cd6d345fb 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1255 +1,1255 @@ /** * @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" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif #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), material_local_numbering("material local numbering", id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("solid_mechanics_model", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif 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); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* 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) { +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."); // 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")) { 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")) { 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 diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 625f149a4..f9008eae3 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,593 +1,593 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Apr 09 2021 * * @brief Model of Solid Mechanics * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager_callback.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_HH_ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "solid_mechanics_model", std::shared_ptr dof_manager = nullptr, ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; public: /// initialize all internal arrays for materials virtual void initMaterials(); protected: /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(bool need_to_reassemble = false); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; - bool canSplitResidual() override { return true; } + bool canSplitResidual() const override { return true; } /// get the type of matrix needed - MatrixType getMatrixType(const ID & matrix_id) override; + MatrixType getMatrixType(const ID & matrix_id) const override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); public: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); protected: /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(ElementType type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(GhostType ghost_type) override; void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, ElementKind element_kind); //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, ElementKind kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, ElementKind kind, GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(ElementKind kind); std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; void dump(const std::string & dumper_name) override; void dump(const std::string & dumper_name, UInt step) override; void dump(const std::string & dumper_name, Real time, UInt step) override; void dump() override; void dump(UInt step) override; void dump(Real time, UInt step) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement); /// get the SolidMechanicsModel::previous_displacement array AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement); /// get the SolidMechanicsModel::current_position array const Array & getCurrentPosition(); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment); /// get the lumped SolidMechanicsModel::mass array AKANTU_GET_MACRO_DEREF_PTR(Mass, mass); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force); /// get the SolidMechanicsModel::force array (external forces) [[deprecated("Use getExternalForce instead of this function")]] Array & getForce() { return getExternalForce(); } /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force); /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs); /// get an iterable on the materials inline decltype(auto) getMaterials(); /// get an iterable on the materials inline decltype(auto) getMaterials() const; /// get a particular material (by numerical material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by numerical material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /** * @brief Returns the total energy for a given energy type * * Energy types of SolidMechanicsModel expected as argument are: * - `kinetic` * - `external work` * * Other energy types are passed on to the materials. All materials should * define a `potential` energy type. For additional energy types, see material * documentation. */ Real getEnergy(const std::string & energy_id); /// Compute energy for an element type and material index Real getEnergy(const std::string & energy_id, ElementType type, UInt index); /// Compute energy for an individual element Real getEnergy(const std::string & energy_id, const Element & element) { return getEnergy(energy_id, element.type, element.element); } /// Compute energy for an element group Real getEnergy(const ID & energy_id, const ID & group_id); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, // material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, material_selector, std::shared_ptr); void setMaterialSelector(std::shared_ptr material_selector) { this->material_selector = std::move(material_selector); } /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: /// compute the stable time step Real getStableTimeStep(GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// release version of the displacement array UInt displacement_release{0}; /// release version of the current_position array UInt current_position_release{0}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// mapping between material name and material internal id std::map materials_names_to_id; protected: /// conversion coefficient form force/mass to acceleration Real f_m2a{1.0}; /// displacements array std::unique_ptr> displacement; /// displacements array at the previous time step (used in finite deformation) std::unique_ptr> previous_displacement; /// increment of displacement std::unique_ptr> displacement_increment; /// lumped mass array std::unique_ptr> mass; /// velocities array std::unique_ptr> velocity; /// accelerations array std::unique_ptr> acceleration; /// external forces array std::unique_ptr> external_force; /// internal forces array std::unique_ptr> internal_force; /// array specifing if a degree of freedom is blocked or not std::unique_ptr> blocked_dofs; /// array of current position used during update residual std::unique_ptr> current_position; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; /// list of used materials std::vector> materials; /// class defining of to choose a material std::shared_ptr material_selector; using flatten_internal_map = std::map, std::unique_ptr>>; /// tells if the material are instantiated flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; /// tells if the material are instantiated bool are_materials_instantiated{false}; friend class Material; template friend class CouplerSolidContactTemplate; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.hh" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */