diff --git a/python/py_solid_mechanics_model.cc b/python/py_solid_mechanics_model.cc index abb461dab..e18541124 100644 --- a/python/py_solid_mechanics_model.cc +++ b/python/py_solid_mechanics_model.cc @@ -1,113 +1,117 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ -void -register_solid_mechanics_model(py::module & mod) { +void register_solid_mechanics_model(py::module & mod) { py::class_(mod, "SolidMechanicsModelOptions") .def(py::init(), py::arg("_analysis_method") = _explicit_lumped_mass); py::class_(mod, "SolidMechanicsModel", py::multiple_inheritance()) .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "solid_mechanics_model", py::arg("memory_id") = 0, py::arg("model_type") = ModelType::_solid_mechanics_model) .def( "initFull", [](SolidMechanicsModel & self, const SolidMechanicsModelOptions & options) { self.initFull(options); }, py::arg("option") = SolidMechanicsModelOptions()) .def( "initFull", [](SolidMechanicsModel & self, const AnalysisMethod & analysis_method) { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method")) .def_deprecated("applyDirichletBC", "Deprecated: use applyBC") .def("applyBC", [](SolidMechanicsModel & self, BC::Dirichlet::DirichletFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("applyBC", [](SolidMechanicsModel & self, BC::Neumann::NeumannFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("setTimeStep", &SolidMechanicsModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def("getEnergy", py::overload_cast( &SolidMechanicsModel::getEnergy), py::arg("energy_id")) + .def("getEnergy", + py::overload_cast( + &SolidMechanicsModel::getEnergy), + py::arg("energy_id"), py::arg("group_id")) + .def_function(assembleStiffnessMatrix) .def_function(assembleInternalForces) .def_function(assembleMass) .def_function(assembleMassLumped) .def_function(getStableTimeStep) .def_function_nocopy(getExternalForce) .def_function_nocopy(getDisplacement) .def_function_nocopy(getPreviousDisplacement) .def_function_nocopy(getIncrement) .def_function_nocopy(getInternalForce) .def_function_nocopy(getMass) .def_function_nocopy(getVelocity) .def_function_nocopy(getAcceleration) .def_function_nocopy(getInternalForce) .def_function_nocopy(getBlockedDOFs) .def_function_nocopy(getMesh) .def("dump", py::overload_cast<>(&SolidMechanicsModel::dump)) .def("dump", py::overload_cast(&SolidMechanicsModel::dump)) .def("dump", py::overload_cast( &SolidMechanicsModel::dump)) .def("dump", py::overload_cast( &SolidMechanicsModel::dump)) .def("getMaterial", py::overload_cast(&SolidMechanicsModel::getMaterial), py::return_value_policy::reference) .def("getMaterial", py::overload_cast( &SolidMechanicsModel::getMaterial), py::return_value_policy::reference) .def("getMaterialIndex", &SolidMechanicsModel::getMaterialIndex) .def("setMaterialSelector", &SolidMechanicsModel::setMaterialSelector) .def("getMaterialSelector", &SolidMechanicsModel::getMaterialSelector); } } // namespace akantu diff --git a/src/mesh/element_group.hh b/src/mesh/element_group.hh index b1a8eaa48..ed3cc443a 100644 --- a/src/mesh/element_group.hh +++ b/src/mesh/element_group.hh @@ -1,202 +1,205 @@ /** * @file element_group.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Wed Nov 08 2017 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "dumpable.hh" #include "element_type_map.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_GROUP_HH_ #define AKANTU_ELEMENT_GROUP_HH_ namespace akantu { class Mesh; class Element; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class ElementGroup : private Memory, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ElementGroup(const std::string & name, const Mesh & mesh, NodeGroup & node_group, UInt dimension = _all_dimensions, const std::string & id = "element_group", const MemoryID & mem_id = 0); ElementGroup(const ElementGroup & /*unused*/); /* ------------------------------------------------------------------------ */ /* Type definitions */ /* ------------------------------------------------------------------------ */ public: using ElementList = ElementTypeMapArray; using NodeList = Array; /* ------------------------------------------------------------------------ */ /* Element iterator */ /* ------------------------------------------------------------------------ */ using type_iterator = ElementList::type_iterator; [[deprecated("Use elementTypes instead")]] inline type_iterator firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const; [[deprecated("Use elementTypes instead")]] inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const; template inline decltype(auto) elementTypes(pack &&... _pack) const { return elements.elementTypes(_pack...); } using const_element_iterator = Array::const_iterator; inline const_element_iterator begin(ElementType type, GhostType ghost_type = _not_ghost) const; inline const_element_iterator end(ElementType type, GhostType ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// empty the element group void clear(); void clear(ElementType type, GhostType ghost_type = _not_ghost); bool empty() const __attribute__((warn_unused_result)); /// append another group to this group /// BE CAREFUL: it doesn't conserve the element order void append(const ElementGroup & other_group); /// add an element to the group. By default the it does not add the nodes to /// the group inline void add(const Element & el, bool add_nodes = false, bool check_for_duplicate = true); /// \todo fix the default for add_nodes : make it coherent with the other /// method inline void add(ElementType type, UInt element, GhostType ghost_type = _not_ghost, bool add_nodes = true, bool check_for_duplicate = true); inline void addNode(UInt node_id, bool check_for_duplicate = true); inline void removeNode(UInt node_id); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// fill the elements based on the underlying node group. virtual void fillFromNodeGroup(); // sort and remove duplicated values void optimize(); /// change the dimension if needed void addDimension(UInt dimension); private: inline void addElement(ElementType elem_type, UInt elem_id, GhostType ghost_type); friend class GroupManager; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: const Array & getElements(ElementType type, GhostType ghost_type = _not_ghost) const; AKANTU_GET_MACRO(Elements, elements, const ElementTypeMapArray &); AKANTU_GET_MACRO_NOT_CONST(Elements, elements, ElementTypeMapArray &); template auto size(Args &&... pack) const { return elements.size(std::forward(pack)...); } + decltype(auto) getElementsIterable(ElementType type, + GhostType ghost_type = _not_ghost) const; + // AKANTU_GET_MACRO(Nodes, node_group.getNodes(), const Array &); AKANTU_GET_MACRO(NodeGroup, node_group, const NodeGroup &); AKANTU_GET_MACRO_NOT_CONST(NodeGroup, node_group, NodeGroup &); AKANTU_GET_MACRO(Dimension, dimension, UInt); AKANTU_GET_MACRO(Name, name, std::string); inline UInt getNbNodes() const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Mesh to which this group belongs const Mesh & mesh; /// name of the group std::string name; /// list of elements composing the group ElementList elements; /// sub list of nodes which are composing the elements NodeGroup & node_group; /// group dimension UInt dimension{_all_dimensions}; /// empty arry for the iterator to work when an element type not present Array empty_elements; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const ElementGroup & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "element.hh" #include "element_group_inline_impl.hh" #endif /* AKANTU_ELEMENT_GROUP_HH_ */ diff --git a/src/mesh/element_group_inline_impl.hh b/src/mesh/element_group_inline_impl.hh index d45a3ffb9..9ade23bc3 100644 --- a/src/mesh/element_group_inline_impl.hh +++ b/src/mesh/element_group_inline_impl.hh @@ -1,143 +1,146 @@ /** * @file element_group_inline_impl.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Sun Aug 13 2017 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "mesh.hh" - +#include "aka_iterators.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_GROUP_INLINE_IMPL_HH_ #define AKANTU_ELEMENT_GROUP_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline void ElementGroup::add(const Element & el, bool add_nodes, bool check_for_duplicate) { this->add(el.type, el.element, el.ghost_type, add_nodes, check_for_duplicate); } /* -------------------------------------------------------------------------- */ inline void ElementGroup::add(ElementType type, UInt element, GhostType ghost_type, bool add_nodes, bool check_for_duplicate) { addElement(type, element, ghost_type); if (add_nodes) { Array::const_vector_iterator it = mesh.getConnectivity(type, ghost_type) .begin(mesh.getNbNodesPerElement(type)) + element; const Vector & conn = *it; for (UInt i = 0; i < conn.size(); ++i) { addNode(conn[i], check_for_duplicate); } } } /* -------------------------------------------------------------------------- */ inline void ElementGroup::addNode(UInt node_id, bool check_for_duplicate) { node_group.add(node_id, check_for_duplicate); } /* -------------------------------------------------------------------------- */ inline void ElementGroup::removeNode(UInt node_id) { node_group.remove(node_id); } /* -------------------------------------------------------------------------- */ -inline void ElementGroup::addElement(ElementType elem_type, - UInt elem_id, +inline void ElementGroup::addElement(ElementType elem_type, UInt elem_id, GhostType ghost_type) { if (!(elements.exists(elem_type, ghost_type))) { elements.alloc(0, 1, elem_type, ghost_type); } elements(elem_type, ghost_type).push_back(elem_id); this->dimension = UInt( std::max(Int(this->dimension), Int(mesh.getSpatialDimension(elem_type)))); } /* -------------------------------------------------------------------------- */ inline UInt ElementGroup::getNbNodes() const { return node_group.size(); } /* -------------------------------------------------------------------------- */ inline ElementGroup::type_iterator ElementGroup::firstType(UInt dim, GhostType ghost_type, ElementKind kind) const { return elements.elementTypes(dim, ghost_type, kind).begin(); } /* -------------------------------------------------------------------------- */ inline ElementGroup::type_iterator -ElementGroup::lastType(UInt dim, GhostType ghost_type, - ElementKind kind) const { +ElementGroup::lastType(UInt dim, GhostType ghost_type, ElementKind kind) const { return elements.elementTypes(dim, ghost_type, kind).end(); } /* -------------------------------------------------------------------------- */ inline ElementGroup::const_element_iterator -ElementGroup::begin(ElementType type, - GhostType ghost_type) const { +ElementGroup::begin(ElementType type, GhostType ghost_type) const { if (elements.exists(type, ghost_type)) { return elements(type, ghost_type).begin(); } return empty_elements.begin(); } /* -------------------------------------------------------------------------- */ inline ElementGroup::const_element_iterator -ElementGroup::end(ElementType type, - GhostType ghost_type) const { +ElementGroup::end(ElementType type, GhostType ghost_type) const { if (elements.exists(type, ghost_type)) { return elements(type, ghost_type).end(); } return empty_elements.end(); } /* -------------------------------------------------------------------------- */ inline const Array & -ElementGroup::getElements(ElementType type, - GhostType ghost_type) const { +ElementGroup::getElements(ElementType type, GhostType ghost_type) const { if (elements.exists(type, ghost_type)) { return elements(type, ghost_type); } return empty_elements; } /* -------------------------------------------------------------------------- */ +inline decltype(auto) +ElementGroup::getElementsIterable(ElementType type, + GhostType ghost_type) const { + return make_transform_adaptor(this->elements(type, ghost_type), + [type, ghost_type](auto && el) { + return Element{type, el, ghost_type}; + }); +} } // namespace akantu #endif /* AKANTU_ELEMENT_GROUP_INLINE_IMPL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 33fa68ec6..7551fcb83 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1240 +1,1258 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei * @author Guillaume Anciaux * @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: Wed Feb 21 2018 * * @brief Implementation of the SolidMechanicsModel class * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "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 memory_id the id of the memory * @param model_type this is an internal parameter for inheritance purposes */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id), material_index("material index", id, memory_id), material_local_numbering("material local numbering", id, memory_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->initDOFManager(); 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) { // \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() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // Check if materials need to recompute the matrix bool need_to_reassemble = false; 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(), this->getMemoryID()); 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, this->memory_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 a5d896289..00921dc07 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,581 +1,583 @@ /** * @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: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "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", const MemoryID & memory_id = 0, - ModelType model_type = ModelType::_solid_mechanics_model); + SolidMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, + const ID & id = "solid_mechanics_model", + const MemoryID & memory_id = 0, + 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(); /// 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; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) 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(); 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); /// 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; + 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; + 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; + 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); + 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); + 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; + UInt spatial_dimension, ElementKind kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// 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(); /// get the energies Real getEnergy(const std::string & energy_id); - /// compute the energy for energy - Real getEnergy(const std::string & energy_id, ElementType type, - UInt index); + /// compute the energy for an element + Real getEnergy(const std::string & energy_id, ElementType type, UInt index); + + /// compute the energy for an element + Real getEnergy(const std::string & energy_id, const Element & element) { + return getEnergy(energy_id, element.type, element.element); + } + + /// compute the 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, MaterialSelector &); - void setMaterialSelector(std::shared_ptr material_selector) { + 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>>; /// map a registered internals to be flattened for dump purposes 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}; }; /* -------------------------------------------------------------------------- */ 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_ */ diff --git a/third-party/akantu_iterators/include/aka_iterators.hh b/third-party/akantu_iterators/include/aka_iterators.hh index fd5e3dd9c..91b676892 100644 --- a/third-party/akantu_iterators/include/aka_iterators.hh +++ b/third-party/akantu_iterators/include/aka_iterators.hh @@ -1,36 +1,37 @@ /** * @file aka_iterators.hh * * @author Nicolas Richart * * @date creation: Fri Aug 11 2017 * @date last modification: Mon Jan 29 2018 * * @brief iterator interfaces * * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * akantu-iterators 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-iterators 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-iterators. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "iterators/aka_arange_iterator.hh" +#include "iterators/aka_concatenate_iterator.hh" #include "iterators/aka_enumerate_iterator.hh" #include "iterators/aka_filter_iterator.hh" #include "iterators/aka_transform_iterator.hh" #include "iterators/aka_zip_iterator.hh" /* -------------------------------------------------------------------------- */ diff --git a/third-party/akantu_iterators/include/aka_tuple_tools.hh b/third-party/akantu_iterators/include/aka_tuple_tools.hh index 31e5e88c4..b6773b150 100644 --- a/third-party/akantu_iterators/include/aka_tuple_tools.hh +++ b/third-party/akantu_iterators/include/aka_tuple_tools.hh @@ -1,337 +1,409 @@ /** * @file aka_tuple_tools.hh * * @author Nicolas Richart * * @date creation: Fri Aug 11 2017 * @date last modification: Mon Jan 29 2018 * * @brief iterator interfaces * * * Copyright 2019 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_compatibilty_with_cpp_standard.hh" #include "aka_str_hash.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_AKA_TUPLE_TOOLS_HH #define AKANTU_AKA_TUPLE_TOOLS_HH #ifndef AKANTU_ITERATORS_NAMESPACE #define AKANTU_ITERATORS_NAMESPACE akantu #endif namespace AKANTU_ITERATORS_NAMESPACE { namespace tuple { /* ---------------------------------------------------------------------- */ template struct named_tag { using _tag = tag; ///< key using _type = type; ///< value type template < typename T, std::enable_if_t::value> * = nullptr> explicit named_tag(T && value) // NOLINT : _value(std::forward(value)) {} type _value; }; namespace details { /* ---------------------------------------------------------------------- */ #if (defined(__GNUC__) || defined(__GNUG__)) #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Weffc++" #endif template struct named_tag_proxy { using _tag = tag; template decltype(auto) operator=(T && value) { return named_tag<_tag, T>{std::forward(value)}; } }; #if (defined(__GNUC__) || defined(__GNUG__)) #pragma GCC diagnostic pop #endif } // namespace details /* ---------------------------------------------------------------------- */ template struct is_named_tag : public std::false_type {}; template struct is_named_tag> : public std::true_type {}; template struct is_named_tag> : public std::true_type {}; /* ---------------------------------------------------------------------- */ template struct named_tuple : public std::tuple { using Names_t = std::tuple; using parent = std::tuple; named_tuple(Params &&... params) : parent(std::forward(params._value)...) {} named_tuple(typename Params::_type &&... args) : parent(std::forward(args)...) {} private: template * = nullptr> static constexpr std::size_t get_element_index() noexcept { return -1; } template * = nullptr> static constexpr std::size_t get_element_index() noexcept { using _tag = std::tuple_element_t; return (std::is_same<_tag, tag>::value) ? Idx : get_element_index(); } public: template ::value> * = nullptr> constexpr decltype(auto) get(NT && /*unused*/) noexcept { const auto index = get_element_index(); static_assert((index != -1), "wrong named_tag"); return (std::get(*this)); } template ::value> * = nullptr> constexpr decltype(auto) get(NT && /*unused*/) const noexcept { const auto index = get_element_index(); static_assert((index != -1), "wrong named_tag"); return std::get(*this); } }; /* ---------------------------------------------------------------------- */ template struct is_named_tuple : public std::false_type {}; template struct is_named_tuple> : public std::true_type {}; /* ---------------------------------------------------------------------- */ template constexpr decltype(auto) make_named_tuple(Params &&... params) noexcept { return named_tuple(std::forward(params)...); } template constexpr decltype(auto) make_named_tag() noexcept { return details::named_tag_proxy{}; } template constexpr decltype(auto) get() { return make_named_tag>(); } template constexpr decltype(auto) get(Tuple && tuple) { return tuple.get(get()); } template ::value> * = nullptr> constexpr decltype(auto) get(Tuple && tuple) noexcept { return tuple.template get(); } #if defined(__INTEL_COMPILER) // intel warnings here #elif defined(__clang__) // clang warnings here #pragma clang diagnostic push #pragma clang diagnostic ignored "-Wgnu-string-literal-operator-template" #elif (defined(__GNUC__) || defined(__GNUG__)) // gcc warnings here #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wpedantic" #endif /// this is a GNU exstension /// http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3599.html template constexpr decltype(auto) operator"" _n() { return make_named_tag::hash>>(); } #if defined(__INTEL_COMPILER) #elif defined(__clang__) #pragma clang diagnostic pop #elif (defined(__GNUC__) || defined(__GNUG__)) #pragma GCC diagnostic pop #endif /* ------------------------------------------------------------------------ */ namespace details { template struct Foreach { template - static inline bool not_equal(Tuple && a, Tuple && b) { + + static constexpr inline auto not_equal(Tuple && a, Tuple && b) -> bool { if (std::get(std::forward(a)) == std::get(std::forward(b))) { return false; } return Foreach::not_equal(std::forward(a), std::forward(b)); } + + template + static constexpr inline auto find(Tuple && tuple, V && value) + -> std::size_t { + constexpr auto size = std::tuple_size>::value; + if (std::get(std::forward(tuple)) == value) { + return size - N; + } + + return Foreach::find(std::forward(tuple), + std::forward(value)); + } }; /* ---------------------------------------------------------------------- */ template <> struct Foreach<0> { template - static inline bool not_equal(Tuple && a, Tuple && b) { + static constexpr inline auto not_equal(Tuple && a, Tuple && b) -> bool { return std::get<0>(std::forward(a)) != std::get<0>(std::forward(b)); } + + template + static constexpr inline auto find(Tuple && tuple, V && value) + -> std::size_t { + constexpr auto size = std::tuple_size>::value; + if (std::get(std::forward(tuple)) == value) { + return size - 1; + } + + return size; + } }; template decltype(auto) make_tuple_no_decay(Ts &&... args) { return std::tuple(std::forward(args)...); } template - decltype(auto) make_named_tuple_no_decay(std::tuple /*unused*/, - Ts &&... args) { + constexpr decltype(auto) + make_named_tuple_no_decay(std::tuple /*unused*/, Ts &&... args) { return named_tuple...>(std::forward(args)...); } template - void foreach_impl(F && func, Tuple && tuple, - std::index_sequence && /*unused*/) { + constexpr void foreach_impl(F && func, Tuple && tuple, + std::index_sequence && /*unused*/) { (void)std::initializer_list{ (std::forward(func)(std::get(std::forward(tuple))), 0)...}; } template - decltype(auto) transform_impl(F && func, Tuple && tuple, - std::index_sequence && /*unused*/) { + constexpr decltype(auto) + transform_impl(F && func, Tuple && tuple, + std::index_sequence && /*unused*/) { return make_tuple_no_decay( std::forward(func)(std::get(std::forward(tuple)))...); } template - decltype(auto) + constexpr decltype(auto) + transform_impl(F && func, Tuple && tuple1, Tuple && tuple2, + std::index_sequence && /*unused*/) { + return make_tuple_no_decay( + std::forward(func)(std::get(std::forward(tuple1)), + std::get(std::forward(tuple2)))...); + } + + template + constexpr decltype(auto) transform_named_impl(F && func, Tuple && tuple, std::index_sequence && /*unused*/) { return make_named_tuple_no_decay( typename std::decay_t::Names_t{}, std::forward(func)(std::get(std::forward(tuple)))...); } } // namespace details - /* ------------------------------------------------------------------------ */ + /* ------------------------------------------------------------------------ + */ template >::value> * = nullptr> - bool are_not_equal(Tuple && a, Tuple && b) { + constexpr auto are_not_equal(Tuple && a, Tuple && b) -> bool { return details::Foreach>::value>:: not_equal(std::forward(a), std::forward(b)); } template < class Tuple, std::enable_if_t>::value> * = nullptr> - bool are_not_equal(Tuple && a, Tuple && b) { + constexpr auto are_not_equal(Tuple && a, Tuple && b) -> bool { return details::Foreach< std::tuple_size::parent>::value>:: not_equal(std::forward(a), std::forward(b)); } + template + constexpr decltype(auto) find(Tuple && tuple, V && value) { + return details::Foreach>::value>::find( + std::forward(tuple), std::forward(value)); + } + template >::value> * = nullptr> - void foreach (F && func, Tuple && tuple) { + constexpr void foreach (F && func, Tuple && tuple) { return details::foreach_impl( std::forward(func), std::forward(tuple), std::make_index_sequence< std::tuple_size>::value>{}); } template < class F, class Tuple, std::enable_if_t>::value> * = nullptr> - void foreach (F && func, Tuple && tuple) { + constexpr void foreach (F && func, Tuple && tuple) { return details::foreach_impl( std::forward(func), std::forward(tuple), std::make_index_sequence< std::tuple_size::parent>::value>{}); } template >::value> * = nullptr> - decltype(auto) transform(F && func, Tuple && tuple) { + constexpr decltype(auto) transform(F && func, Tuple && tuple) { return details::transform_impl( std::forward(func), std::forward(tuple), std::make_index_sequence< std::tuple_size>::value>{}); } + template >::value> * = + nullptr> + constexpr decltype(auto) transform(F && func, Tuple && tuple1, + Tuple && tuple2) { + return details::transform_impl( + std::forward(func), std::forward(tuple1), + std::forward(tuple2), + std::make_index_sequence< + std::tuple_size>::value>{}); + } + template < class F, class Tuple, std::enable_if_t>::value> * = nullptr> - decltype(auto) transform(F && func, Tuple && tuple) { + constexpr decltype(auto) transform(F && func, Tuple && tuple) { return details::transform_named_impl( std::forward(func), std::forward(tuple), std::make_index_sequence< std::tuple_size::parent>::value>{}); } namespace details { template - decltype(auto) flatten(Tuple && tuples, - std::index_sequence /*unused*/) { + constexpr decltype(auto) flatten(Tuple && tuples, + std::index_sequence /*unused*/) { return std::tuple_cat(std::get(tuples)...); } } // namespace details - template decltype(auto) flatten(Tuple && tuples) { + template constexpr decltype(auto) flatten(Tuple && tuples) { return details::flatten(std::forward(tuples), std::make_index_sequence< std::tuple_size>::value>()); } + + namespace details { + template + decltype(auto) dynamic_get_impl(size_t i, Tuple && tuple) { + constexpr auto size = std::tuple_size>::value; + if (i == n) { + return std::get(tuple); + } else if (n == size - 1) { + throw std::out_of_range("Tuple element out of range."); + } else { + return dynamic_get_impl<(n < size - 1 ? n + 1 : 0)>(i, tuple); + } + } + } // namespace details + + template + constexpr decltype(auto) dynamic_get(std::size_t i, Tuple && tuple) { + return details::dynamic_get_impl<0>(i, std::forward(tuple)); + } } // namespace tuple } // namespace AKANTU_ITERATORS_NAMESPACE /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace std { template struct iterator_traits< ::AKANTU_ITERATORS_NAMESPACE::tuple::named_tag> { using iterator_category = typename type::iterator_category; using value_type = typename type::value_type; using difference_type = typename type::difference_type; using pointer = typename type::pointer; using reference = typename type::reference; }; } // namespace std #endif /* AKANTU_AKA_TUPLE_TOOLS_HH */ diff --git a/third-party/akantu_iterators/include/iterators/aka_concatenate_iterator.hh b/third-party/akantu_iterators/include/iterators/aka_concatenate_iterator.hh new file mode 100644 index 000000000..a7948c5a9 --- /dev/null +++ b/third-party/akantu_iterators/include/iterators/aka_concatenate_iterator.hh @@ -0,0 +1,205 @@ +/** + * @file aka_concatenate_iterator.hh + * + * @author Nicolas Richart + * + * @date creation jeu déc 12 2019 + * + * @brief implementation of arange + * + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * akantu-iterators 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-iterators 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-iterators. If not, see . + * + */ +/* -------------------------------------------------------------------------- */ +#include "aka_compatibilty_with_cpp_standard.hh" +#include "aka_iterator_tools.hh" +#include "aka_tuple_tools.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +/* -------------------------------------------------------------------------- */ + +#ifndef AKA_CONCATENATE_ITERATOR_H +#define AKA_CONCATENATE_ITERATOR_H + +#ifndef AKANTU_ITERATORS_NAMESPACE +#define AKANTU_ITERATORS_NAMESPACE akantu +#endif + +namespace AKANTU_ITERATORS_NAMESPACE { + +namespace iterators { + + /* ------------------------------------------------------------------------ */ + template + class ConcatIterator + : public details::CopyAssignmentEnabler< + aka::conjunction..., + std::is_copy_constructible...>::value>, + public details::MoveAssignmentEnabler< + aka::conjunction..., + std::is_move_constructible...>::value> { + private: + using tuple_t = std::tuple; + + public: + using value_type = + std::tuple::value_type...>; + using difference_type = std::common_type_t< + typename std::iterator_traits::difference_type...>; + using pointer = + std::tuple::pointer...>; + using reference = + std::tuple::reference...>; + using iterator_category = std::input_iterator_tag; + + public: + explicit ConcatIterator(tuple_t iterators, tuple_t end_iterators) + : iterators(std::move(iterators)), + end_iterators(std::move(end_iterators)) {} + + /* ---------------------------------------------------------------------- */ + // input iterator ++it + ConcatIterator & operator++() { + auto && ends = + tuple::transform([](auto && a, auto && b) { return a == b; }, + iterators, end_iterators); + auto && pos = tuple::find(ends, false); + ++(tuple::dynamic_get(pos, iterators)); + return *this; + } + + // input iterator it++ + ConcatIterator operator++(int) { + auto cpy = *this; + this->operator++(); + return cpy; + } + + // input iterator it != other_it + bool operator!=(const ConcatIterator & other) const { + return tuple::are_not_equal(iterators, other.iterators); + } + + // input iterator dereference *it + decltype(auto) operator*() { + auto && ends = + tuple::transform([](auto && a, auto && b) { return a == b; }, + iterators, end_iterators); + auto && pos = tuple::find(ends, false); + return *(tuple::dynamic_get(pos, iterators)); + } + + template < + class iterator_category_ = iterator_category, + std::enable_if_t::value> * = nullptr> + bool operator==(const ConcatIterator & other) const { + return not tuple::are_not_equal(iterators, other.iterators); + } + + private: + tuple_t iterators; + tuple_t end_iterators; + }; + +} // namespace iterators + +/* -------------------------------------------------------------------------- */ +template +decltype(auto) +concat_iterator(std::tuple && iterators_tuple, + std::tuple && end_iterators_tuple) { + auto concat = iterators::ConcatIterator( + std::forward(iterators_tuple), + std::forward(end_iterators_tuple)); + return concat; +} + +/* -------------------------------------------------------------------------- */ +namespace containers { + template class ConcatContainer { + using containers_t = std::tuple; + + public: + explicit ConcatContainer(Containers &&... containers) + : containers(std::forward(containers)...) {} + + decltype(auto) begin() const { + return concat_iterator( + tuple::transform([](auto && c) { return c.begin(); }, + std::forward(containers)), + tuple::transform([](auto && c) { return c.end(); }, + std::forward(containers))); + } + + decltype(auto) end() const { + return concat_iterator( + tuple::transform([](auto && c) { return c.end(); }, + std::forward(containers)), + tuple::transform([](auto && c) { return c.end(); }, + std::forward(containers))); + } + + decltype(auto) begin() { + return concat_iterator( + tuple::transform([](auto && c) { return c.begin(); }, + std::forward(containers)), + tuple::transform([](auto && c) { return c.end(); }, + std::forward(containers))); + } + + decltype(auto) end() { + return concat_iterator( + tuple::transform([](auto && c) { return c.end(); }, + std::forward(containers)), + tuple::transform([](auto && c) { return c.end(); }, + std::forward(containers))); + } + + private: + containers_t containers; + }; + + /* ------------------------------------------------------------------------ */ + template decltype(auto) concat(Containers &&... conts) { + return containers::ConcatContainer( + std::forward(conts)...); + } + +} // namespace containers +} // namespace AKANTU_ITERATORS_NAMESPACE + +namespace std { +template +struct iterator_traits< + ::AKANTU_ITERATORS_NAMESPACE::iterators::ConcatIterator> { +private: + using iterator_type = + typename ::AKANTU_ITERATORS_NAMESPACE::iterators::ConcatIterator; + +public: + using iterator_category = typename iterator_type::iterator_category; + using value_type = typename iterator_type::value_type; + using difference_type = typename iterator_type::difference_type; + using pointer = typename iterator_type::pointer; + using reference = typename iterator_type::reference; +}; +} // namespace std + +#endif // __AKA_CONCATENATE_ITERATOR_H_ diff --git a/third-party/akantu_iterators/include/iterators/aka_iterator_tools.hh b/third-party/akantu_iterators/include/iterators/aka_iterator_tools.hh new file mode 100644 index 000000000..aec788231 --- /dev/null +++ b/third-party/akantu_iterators/include/iterators/aka_iterator_tools.hh @@ -0,0 +1,69 @@ +/** + * @file aka_iterator_tools.hh + * + * @author Nicolas Richart + * + * @date creation jeu déc 12 2019 + * + * @brief A Documented file. + * + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * akantu-iterators 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-iterators 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-iterators. If not, see . + * + */ +/* -------------------------------------------------------------------------- */ + +#ifndef AKA_ITERATOR_TOOLS_H +#define AKA_ITERATOR_TOOLS_H + +#ifndef AKANTU_ITERATORS_NAMESPACE +#define AKANTU_ITERATORS_NAMESPACE akantu +#endif + +namespace AKANTU_ITERATORS_NAMESPACE { + +/* -------------------------------------------------------------------------- */ +namespace iterators { + + namespace details { + template struct CopyAssignmentEnabler {}; + template <> struct CopyAssignmentEnabler { + CopyAssignmentEnabler() = default; + CopyAssignmentEnabler(const CopyAssignmentEnabler &) = default; + CopyAssignmentEnabler(CopyAssignmentEnabler &&) = default; + auto operator=(const CopyAssignmentEnabler &) + -> CopyAssignmentEnabler & = delete; + auto operator=(CopyAssignmentEnabler &&) + -> CopyAssignmentEnabler & = default; + }; + + template struct MoveAssignmentEnabler {}; + template <> struct MoveAssignmentEnabler { + MoveAssignmentEnabler() = default; + MoveAssignmentEnabler(const MoveAssignmentEnabler &) = default; + MoveAssignmentEnabler(MoveAssignmentEnabler &&) = default; + auto operator=(const MoveAssignmentEnabler &) + -> MoveAssignmentEnabler & = delete; + auto operator=(MoveAssignmentEnabler &&) + -> MoveAssignmentEnabler & = default; + }; + + } // namespace details +} // namespace iterators +} // namespace AKANTU_ITERATORS_NAMESPACE + +#endif // AKA_ITERATOR_TOOLS_H diff --git a/third-party/akantu_iterators/include/iterators/aka_zip_iterator.hh b/third-party/akantu_iterators/include/iterators/aka_zip_iterator.hh index 3a82dc04b..30c338996 100644 --- a/third-party/akantu_iterators/include/iterators/aka_zip_iterator.hh +++ b/third-party/akantu_iterators/include/iterators/aka_zip_iterator.hh @@ -1,303 +1,282 @@ /** * @file aka_zip_iterator.hh * * @author Nicolas Richart * * @date creation jeu déc 12 2019 * * @brief A Documented file. * * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * akantu-iterators 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-iterators 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-iterators. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_compatibilty_with_cpp_standard.hh" +#include "aka_iterator_tools.hh" #include "aka_tuple_tools.hh" /* -------------------------------------------------------------------------- */ #include -#include #include /* -------------------------------------------------------------------------- */ #ifndef AKA_ZIP_ITERATOR_HH #define AKA_ZIP_ITERATOR_HH #ifndef AKANTU_ITERATORS_NAMESPACE #define AKANTU_ITERATORS_NAMESPACE akantu #endif namespace AKANTU_ITERATORS_NAMESPACE { /* -------------------------------------------------------------------------- */ namespace iterators { - namespace details { - template struct CopyAssignmentEnabler {}; - template <> struct CopyAssignmentEnabler { - CopyAssignmentEnabler() = default; - CopyAssignmentEnabler(const CopyAssignmentEnabler &) = default; - CopyAssignmentEnabler(CopyAssignmentEnabler &&) = default; - CopyAssignmentEnabler & operator=(const CopyAssignmentEnabler &) = delete; - CopyAssignmentEnabler & operator=(CopyAssignmentEnabler &&) = default; - }; - - template struct MoveAssignmentEnabler {}; - template <> struct MoveAssignmentEnabler { - MoveAssignmentEnabler() = default; - MoveAssignmentEnabler(const MoveAssignmentEnabler &) = default; - MoveAssignmentEnabler(MoveAssignmentEnabler &&) = default; - MoveAssignmentEnabler & operator=(const MoveAssignmentEnabler &) = delete; - MoveAssignmentEnabler & operator=(MoveAssignmentEnabler &&) = default; - }; - - } // namespace details - /* ------------------------------------------------------------------------ */ template