diff --git a/python/py_fe_engine.cc b/python/py_fe_engine.cc index 2330ab24f..eab1dd274 100644 --- a/python/py_fe_engine.cc +++ b/python/py_fe_engine.cc @@ -1,291 +1,298 @@ /** * @file py_fe_engine.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Wed Nov 27 2019 * @date last modification: Sat Dec 12 2020 * * @brief pybind11 interface to FEEngine * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" #include "py_aka_common.hh" /* -------------------------------------------------------------------------- */ -#include #include #include -#include +#include +#ifdef AKANTU_COHESIVE_ELEMENT +#include +#endif /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { void register_fe_engine(py::module & mod) { py::class_(mod, "Element") .def(py::init([](ElementType type, UInt id) { return new Element{type, id, _not_ghost}; })) .def(py::init([](ElementType type, UInt id, GhostType ghost_type) { return new Element{type, id, ghost_type}; })) .def("__lt__", [](Element & self, const Element & other) { return (self < other); }) .def("__repr__", [](Element & self) { return std::to_string(self); }); mod.attr("ElementNull") = ElementNull; py::class_(mod, "FEEngine") .def( "getNbIntegrationPoints", [](FEEngine & fem, const ElementType & type, const GhostType & ghost_type) { return fem.getNbIntegrationPoints(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "gradientOnIntegrationPoints", [](FEEngine & fem, const Array & u, Array & nablauq, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type, const Array * filter_elements) { if (filter_elements == nullptr) { // This is due to the ArrayProxy that looses the // empty_filter information filter_elements = &empty_filter; } fem.gradientOnIntegrationPoints(u, nablauq, nb_degree_of_freedom, type, ghost_type, *filter_elements); }, py::arg("u"), py::arg("nablauq"), py::arg("nb_degree_of_freedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) .def( "gradientOnIntegrationPoints", [](FEEngine & fem, const Array & u, Array & nablauq, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type) { fem.gradientOnIntegrationPoints(u, nablauq, nb_degree_of_freedom, type, ghost_type); }, py::arg("u"), py::arg("nablauq"), py::arg("nb_degree_of_freedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "interpolateOnIntegrationPoints", [](FEEngine & self, Array & u, Array & uq, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type) { self.interpolateOnIntegrationPoints(u, uq, nb_degree_of_freedom, type, ghost_type); }, py::arg("u"), py::arg("uq"), py::arg("nb_degree_of_freedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "interpolateOnIntegrationPoints", [](FEEngine & self, Array & u, Array & uq, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type, const Array * filter_elements) { if (filter_elements == nullptr) { // This is due to the ArrayProxy that looses the // empty_filter information filter_elements = &empty_filter; } self.interpolateOnIntegrationPoints(u, uq, nb_degree_of_freedom, type, ghost_type, *filter_elements); }, py::arg("u"), py::arg("uq"), py::arg("nb_degree_of_freedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) .def( "interpolateOnIntegrationPoints", [](FEEngine & self, const Array & u, ElementTypeMapArray & uq, const ElementTypeMapArray * filter_elements) { self.interpolateOnIntegrationPoints(u, uq, filter_elements); }, py::arg("u"), py::arg("uq"), py::arg("filter_elements") = nullptr) .def( "computeIntegrationPointsCoordinates", [](FEEngine & self, ElementTypeMapArray & coordinates, const ElementTypeMapArray * filter_elements) -> decltype(auto) { return self.computeIntegrationPointsCoordinates(coordinates, filter_elements); }, py::arg("coordinates"), py::arg("filter_elements") = nullptr) .def( "computeIntegrationPointsCoordinates", [](FEEngine & self, Array & coordinates, ElementType type, GhostType ghost_type, const Array & filter_elements) -> decltype(auto) { return self.computeIntegrationPointsCoordinates( coordinates, type, ghost_type, filter_elements); }, py::arg("coordinates"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) .def( "assembleFieldLumped", [](FEEngine & fem, const std::function &, const Element &)> & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, ElementType type, GhostType ghost_type) { fem.assembleFieldLumped(field_funct, matrix_id, dof_id, dof_manager, type, ghost_type); }, py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"), py::arg("dof_manager"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "assembleFieldMatrix", [](FEEngine & fem, const std::function &, const Element &)> & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, ElementType type, GhostType ghost_type = _not_ghost) { fem.assembleFieldMatrix(field_funct, matrix_id, dof_id, dof_manager, type, ghost_type); }, py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"), py::arg("dof_manager"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("getElementInradius", [](FEEngine & self, const Element & element) { return self.getElementInradius(element); }) .def("getNormalsOnIntegrationPoints", &FEEngine::getNormalsOnIntegrationPoints, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "computeNtb", [](FEEngine & fem, const Array & bs, Array & Ntbs, ElementType type, GhostType ghost_type) { fem.computeNtb(bs, Ntbs, type, ghost_type); }, py::arg("bs"), py::arg("Ntbs"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "computeNtb", [](FEEngine & fem, const Array & bs, Array & Ntbs, ElementType type, GhostType ghost_type, const Array * filter_elements) { if (filter_elements == nullptr) { // This is due to the ArrayProxy that looses the // empty_filter information filter_elements = &empty_filter; } fem.computeNtb(bs, Ntbs, type, ghost_type, *filter_elements); }, py::arg("Ds"), py::arg("BtDs"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) .def( "computeBtD", [](FEEngine & fem, const Array & Ds, Array & BtDs, ElementType type, GhostType ghost_type) { fem.computeBtD(Ds, BtDs, type, ghost_type); }, py::arg("Ds"), py::arg("BtDs"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "computeBtD", [](FEEngine & fem, const Array & Ds, Array & BtDs, ElementType type, GhostType ghost_type, const Array * filter_elements) { if (filter_elements == nullptr) { // This is due to the ArrayProxy that looses the // empty_filter information filter_elements = &empty_filter; } fem.computeBtD(Ds, BtDs, type, ghost_type, *filter_elements); }, py::arg("Ds"), py::arg("BtDs"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) - // .def("getShapeFunctions", &FEEngine::getShapeFunctionsInterface, - // py::return_value_policy::reference) + .def("getShapeFunctions", &FEEngine::getShapeFunctionsInterface, + py::return_value_policy::reference) .def("getShapes", &FEEngine::getShapes, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("id") = 0, py::return_value_policy::reference) .def("getShapesDerivatives", &FEEngine::getShapesDerivatives, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("id") = 0, py::return_value_policy::reference) .def( "integrate", [](FEEngine & fem, const Array & f, Array & intf, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type) { fem.integrate(f, intf, nb_degree_of_freedom, type, ghost_type); }, py::arg("f"), py::arg("intf"), py::arg("nb_degree_of_feedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "integrate", [](FEEngine & fem, const Array & f, Array & intf, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type, const Array * filter_elements) { if (filter_elements == nullptr) { // This is due to the ArrayProxy that looses the // empty_filter information filter_elements = &empty_filter; } fem.integrate(f, intf, nb_degree_of_freedom, type, ghost_type, *filter_elements); }, py::arg("f"), py::arg("intf"), py::arg("nb_degree_of_feedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) .def("getCohesiveElementType", &FEEngine::getCohesiveElementType); py::class_(mod, "ShapeFunctions"); py::class_(mod, "ShapeLagrangeBase"); + +#ifdef AKANTU_COHESIVE_ELEMENT py::class_, ShapeLagrangeBase>( - mod, "ShapeLagrangeCohesive"); - // .def( - // "interpolateOnIntegrationPointsReduceFunctionDifference", - // [](ShapeLagrange<_ek_cohesive> & self, const Array & field, - // Array & field_difference, UInt nb_degree_of_freedom, - // ElementType type, GhostType ghost_type) { + mod, "ShapeLagrangeCohesive") + .def( + "interpolateOnIntegrationPointsReduceFunctionDifference", + [](ShapeLagrange<_ek_cohesive> & self, const Array & field, + Array & field_difference, UInt nb_degree_of_freedom, + ElementType type, GhostType ghost_type) { + +#define COMPUTE_DIFFERENCE(type) \ + self.interpolateOnIntegrationPoints( \ + field, field_difference, nb_degree_of_freedom, ghost_type); + AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_DIFFERENCE); - // #define COMPUTE_DIFFERENCE(type) \ -// self.interpolateOnIntegrationPoints( \ -// field, field_difference, nb_degree_of_freedom, ghost_type); - // AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_DIFFERENCE); - // #undef COMPUTE_OPENING - // }, - // py::arg("field"), py::arg("field_difference"), - // py::arg("nb_degree_of_freedom"), py::arg("type"), - // py::arg("ghost_type") = _not_ghost); +#undef COMPUTE_OPENING + }, + py::arg("field"), py::arg("field_difference"), + py::arg("nb_degree_of_freedom"), py::arg("type"), + py::arg("ghost_type") = _not_ghost); +#endif + py::class_(mod, "IntegrationPoint"); } } // namespace akantu diff --git a/python/py_heat_transfer_interface_model.cc b/python/py_heat_transfer_interface_model.cc index 735b7da2a..c97c07340 100644 --- a/python/py_heat_transfer_interface_model.cc +++ b/python/py_heat_transfer_interface_model.cc @@ -1,105 +1,102 @@ /** * @file py_heat_transfer_interface_model.cc * * @author Emil Gallyamov * * @date creation: Tue May 16 2023 * @date last modification: Tue May 16 2023 * * @brief pybind11 interface to HeatTransferInterfaceModel * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include -#include /* -------------------------------------------------------------------------- */ // #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, \ [](HeatTransferInterfaceModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](HeatTransferInterfaceModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ void register_heat_transfer_interface_model(py::module & mod) { py::class_( mod, "HeatTransferInterfaceModel") .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "heat_transfer_model") .def( "initFull", [](HeatTransferInterfaceModel & self, const HeatTransferModelOptions & options) { self.initFull(options); }, py::arg("_analysis_method") = HeatTransferModelOptions()) .def( "initFull", [](HeatTransferInterfaceModel & self, const AnalysisMethod & _analysis_method) { self.initFull(HeatTransferModelOptions(_analysis_method)); }, py::arg("_analysis_method")) .def("setTimeStep", &HeatTransferInterfaceModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") - .def("getShapeFunctionsCohesive", - &HeatTransferInterfaceModel::getShapeFunctionsCohesive) .def("getTransversalConductivityOnQpoints", &HeatTransferInterfaceModel::getTransversalConductivityOnQpoints, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def("getLongitudinalConductivityOnQpoints", &HeatTransferInterfaceModel::getLongitudinalConductivityOnQpoints, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getOpening", [](HeatTransferInterfaceModel & self, ElementType & el_type, GhostType & ghost_type) { self.getOpening(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference); } } // namespace akantu diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh index 187f901c4..9b2d299dd 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,491 +1,497 @@ /** * @file element_type_map.hh * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Thu Mar 11 2021 * * @brief storage class by element type * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_named_argument.hh" #include "element.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_TYPE_MAP_HH_ #define AKANTU_ELEMENT_TYPE_MAP_HH_ namespace akantu { class FEEngine; } // namespace akantu namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(all_ghost_types); DECLARE_NAMED_ARGUMENT(default_value); DECLARE_NAMED_ARGUMENT(element_kind); DECLARE_NAMED_ARGUMENT(ghost_type); DECLARE_NAMED_ARGUMENT(nb_component); DECLARE_NAMED_ARGUMENT(nb_component_functor); DECLARE_NAMED_ARGUMENT(with_nb_element); DECLARE_NAMED_ARGUMENT(with_nb_nodes_per_element); DECLARE_NAMED_ARGUMENT(spatial_dimension); DECLARE_NAMED_ARGUMENT(do_not_default); DECLARE_NAMED_ARGUMENT(element_filter); } // namespace template class ElementTypeMap; /* -------------------------------------------------------------------------- */ /* ElementTypeMapBase */ /* -------------------------------------------------------------------------- */ /// Common non templated base class for the ElementTypeMap class class ElementTypeMapBase { public: virtual ~ElementTypeMapBase() = default; }; /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template class ElementTypeMap : public ElementTypeMapBase { public: ElementTypeMap(); ~ElementTypeMap() override; inline static std::string printType(const SupportType & type, GhostType ghost_type); /*! Tests whether a type is present in the object * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return true if the type is present. */ inline bool exists(const SupportType & type, GhostType ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline const Stored & operator()(const SupportType & type, GhostType ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline Stored & operator()(const SupportType & type, GhostType ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. THIS METHOD IS * NOT ARRAY SAFE, when using ElementTypeMapArray, use setArray instead * @param data to insert * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ template inline Stored & operator()(U && insertee, const SupportType & type, GhostType ghost_type = _not_ghost); public: /// print helper virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ /*! iterator allows to iterate over type-data pairs of the map. The interface * expects the SupportType to be ElementType. */ using DataMap = std::map; /// helper class to use in range for constructions class type_iterator { public: using value_type = const SupportType; using pointer = const SupportType *; using reference = const SupportType &; protected: using DataMapIterator = typename ElementTypeMap::DataMap::const_iterator; public: type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek); type_iterator(const type_iterator & it); type_iterator() = default; inline reference operator*(); inline reference operator*() const; inline type_iterator & operator++(); type_iterator operator++(int); inline bool operator==(const type_iterator & other) const; inline bool operator!=(const type_iterator & other) const; type_iterator & operator=(const type_iterator & it); private: DataMapIterator list_begin; DataMapIterator list_end; UInt dim; ElementKind kind; }; /// helper class to use in range for constructions class ElementTypesIteratorHelper { public: using Container = ElementTypeMap; using iterator = typename Container::type_iterator; ElementTypesIteratorHelper(const Container & container, UInt dim, GhostType ghost_type, ElementKind kind) : container(std::cref(container)), dim(dim), ghost_type(ghost_type), kind(kind) {} template ElementTypesIteratorHelper(const Container & container, use_named_args_t /*unused*/, pack &&... _pack) : ElementTypesIteratorHelper( container, OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), OPTIONAL_NAMED_ARG(ghost_type, _not_ghost), OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined)) {} ElementTypesIteratorHelper(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(ElementTypesIteratorHelper &&) noexcept = default; + template ::value> * = nullptr> iterator begin(); + template ::value> * = nullptr> iterator end(); private: std::reference_wrapper container; UInt dim; GhostType ghost_type; ElementKind kind; }; private: ElementTypesIteratorHelper elementTypesImpl(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; template ElementTypesIteratorHelper elementTypesImpl(const use_named_args_t & /*unused*/, pack &&... _pack) const; public: /*! * \param _pack * \parblock * represent optional parameters: * \li \c _spatial_dimension filter for elements of given spatial * dimension * \li \c _ghost_type filter for a certain ghost_type * \li \c _element_kind filter for elements of given kind * \endparblock */ template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(use_named_args, std::forward(_pack)...); } template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(std::forward(_pack)...); } /*! Get an iterator to the beginning of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the first stored data matching the filters * or an iterator to the end of the map if none match*/ [[deprecated("Use elementTypes instead")]] inline type_iterator firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; /*! Get an iterator to the end of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the last stored data matching the filters * or an iterator to the end of the map if none match */ [[deprecated("Use elementTypes instead")]] inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline DataMap & getData(GhostType ghost_type); /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline const DataMap & getData(GhostType ghost_type) const; /* ------------------------------------------------------------------------ */ protected: DataMap data; DataMap ghost_data; }; /* -------------------------------------------------------------------------- */ /* Some typedefs */ /* -------------------------------------------------------------------------- */ template class ElementTypeMapArray : public ElementTypeMap>, SupportType> { public: using value_type = T; using array_type = Array; protected: using parent = ElementTypeMap>, SupportType>; using DataMap = typename parent::DataMap; public: using type_iterator = typename parent::type_iterator; /// standard assigment (copy) operator auto operator=(const ElementTypeMapArray & other) -> ElementTypeMapArray &; ElementTypeMapArray(const ElementTypeMapArray & other); /// explicit copy void copy(const ElementTypeMapArray & other); /*! Constructor * @param id optional: identifier (string) * @param parent_id optional: parent identifier. for organizational purposes * only */ ElementTypeMapArray(const ID & id = "by_element_type_array", const ID & parent_id = "no_parent") : parent(), id(parent_id + ":" + id), name(id){}; /*! allocate memory for a new array * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map * @param ghost_type whether to add the field to the data map or the * ghost_data map * @param default_value the default value to use to fill the array * @return a reference to the allocated array */ inline Array & alloc(UInt size, UInt nb_component, const SupportType & type, GhostType ghost_type, const T & default_value = T()); /*! allocate memory for a new array in both the data and the ghost_data map * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map * @param default_value the default value to use to fill the array */ inline void alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value = T()); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a reference to the array */ inline const Array & operator()(const SupportType & type, GhostType ghost_type = _not_ghost) const; /// access the data of an element, this combine the map and array accessor inline const T & operator()(const Element & element, UInt component = 0) const; /// access the data of an element, this combine the map and array accessor inline T & operator()(const Element & element, UInt component = 0); /// access the data of an element, this combine the map and array accessor inline decltype(auto) get(const Element & element); inline decltype(auto) get(const Element & element) const; /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a const reference to the array */ inline Array & operator()(const SupportType & type, GhostType ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @param vect the vector to include into the map * @return stored data corresponding to type. */ inline void setArray(const SupportType & type, GhostType ghost_type, const Array & vect); /*! frees all memory related to the data*/ inline void free(); inline void clear(); inline bool empty() const __attribute__((warn_unused_result)); /*! set all values in the ElementTypeMap to zero*/ inline void zero() { this->set(T()); } /*! set all values in the ElementTypeMap to value */ template inline void set(const ST & value); /*! deletes and reorders entries in the stored arrays * @param new_numbering a ElementTypeMapArray of new indices. UInt(-1) * indicates * deleted entries. */ inline void onElementsRemoved(const ElementTypeMapArray & new_numbering); /// text output helper void printself(std::ostream & stream, int indent = 0) const override; /*! set the id * @param id the new name */ inline void setID(const ID & id) { this->id = id; } /// return the id inline auto getID() const -> ID { return this->id; } ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType requested_ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const { ElementTypeMap nb_components; bool all_ghost_types = requested_ghost_type == _casper; for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) { continue; } for (auto & type : this->elementTypes(dim, ghost_type, kind)) { UInt nb_comp = (*this)(type, ghost_type).getNbComponent(); nb_components(type, ghost_type) = nb_comp; } } return nb_components; } /* ------------------------------------------------------------------------ */ /* more evolved allocators */ /* ------------------------------------------------------------------------ */ public: /// initialize the arrays in accordance to a functor template void initialize(const Func & f, const T & default_value, bool do_not_default); /// initialize with sizes and number of components in accordance of a mesh /// content template void initialize(const Mesh & mesh, pack &&... _pack); /// initialize with sizes and number of components in accordance of a fe /// engine content (aka integration points) template void initialize(const FEEngine & fe_engine, pack &&... _pack); /* ------------------------------------------------------------------------ */ /* Accesssors */ /* ------------------------------------------------------------------------ */ public: /// get the name of the internal field AKANTU_GET_MACRO(Name, name, ID); /** * get the size of the ElementTypeMapArray * @param[in] _pack * \parblock * optional arguments can be any of: * \li \c _spatial_dimension the dimension to consider (default: * _all_dimensions) * \li \c _ghost_type (default: _not_ghost) * \li \c _element_kind (default: _ek_not_defined) * \li \c _all_ghost_types (default: false) * \endparblock **/ template UInt size(pack &&... _pack) const; bool isNodal() const { return is_nodal; } void isNodal(bool is_nodal) { this->is_nodal = is_nodal; } private: UInt sizeImpl(UInt spatial_dimension, GhostType ghost_type, ElementKind kind) const; private: ID id; protected: /// name of the element type map: e.g. connectivity, grad_u ID name; /// Is the data stored by node of the element bool is_nodal{false}; }; /// to store data Array by element type using ElementTypeMapReal = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapInt = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapUInt = ElementTypeMapArray; } // namespace akantu #endif /* AKANTU_ELEMENT_TYPE_MAP_HH_ */ diff --git a/src/mesh/element_type_map_tmpl.hh b/src/mesh/element_type_map_tmpl.hh index eb6ac82c5..8e872a66c 100644 --- a/src/mesh/element_type_map_tmpl.hh +++ b/src/mesh/element_type_map_tmpl.hh @@ -1,867 +1,870 @@ /** * @file element_type_map_tmpl.hh * * @author Guillaume Anciaux * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Thu Mar 11 2021 * * @brief implementation of template functions of the ElementTypeMap and * ElementTypeMapArray classes * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_static_if.hh" #include "element_type_map.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include "element_type_conversion.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_TYPE_MAP_TMPL_HH_ #define AKANTU_ELEMENT_TYPE_MAP_TMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template inline std::string ElementTypeMap::printType(const SupportType & type, GhostType ghost_type) { std::stringstream sstr; sstr << "(" << ghost_type << ":" << type << ")"; return sstr.str(); } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::exists(const SupportType & type, GhostType ghost_type) const { return this->getData(ghost_type).find(type) != this->getData(ghost_type).end(); } /* -------------------------------------------------------------------------- */ template inline const Stored & ElementTypeMap::operator()(const SupportType & type, GhostType ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) { AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMap::printType(type, ghost_type) << " in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); } return it->second; } /* -------------------------------------------------------------------------- */ template inline Stored & ElementTypeMap::operator()(const SupportType & type, GhostType ghost_type) { return this->getData(ghost_type)[type]; } /* -------------------------------------------------------------------------- */ template template inline Stored & ElementTypeMap::operator()( U && insertee, const SupportType & type, GhostType ghost_type) { auto it = this->getData(ghost_type).find(type); if (it != this->getData(ghost_type).end()) { AKANTU_SILENT_EXCEPTION("Element of type " << ElementTypeMap::printType(type, ghost_type) << " already in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); } else { auto & data = this->getData(ghost_type); const auto & res = data.insert(std::make_pair(type, std::forward(insertee))); it = res.first; } return it->second; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) { if (ghost_type == _not_ghost) { return data; } return ghost_data; } /* -------------------------------------------------------------------------- */ template inline const typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) const { if (ghost_type == _not_ghost) { return data; } return ghost_data; } /* -------------------------------------------------------------------------- */ /// Works only if stored is a pointer to a class with a printself method template void ElementTypeMap::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> [" << std::endl; for (auto && gt : ghost_types) { const DataMap & data = getData(gt); for (auto & pair : data) { stream << space << space << ElementTypeMap::printType(pair.first, gt) << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template ElementTypeMap::ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ template ElementTypeMap::~ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ /* ElementTypeMapArray */ /* -------------------------------------------------------------------------- */ template void ElementTypeMapArray::copy( const ElementTypeMapArray & other) { for (auto ghost_type : ghost_types) { - for (auto type : - this->elementTypes(_all_dimensions, ghost_type, _ek_not_defined)) { - const auto & array_to_copy = other(type, ghost_type); + const DataMap & data = other.getData(ghost_type); + for (auto && [type, array_to_copy] : data) { auto & array = - this->alloc(0, array_to_copy.getNbComponent(), type, ghost_type); - array.copy(array_to_copy); + this->alloc(0, array_to_copy->getNbComponent(), type, ghost_type); + array.copy(*array_to_copy); } } } /* -------------------------------------------------------------------------- */ template ElementTypeMapArray::ElementTypeMapArray( const ElementTypeMapArray & other) : parent(), id(other.id + "_copy"), name(other.name + "_copy") { this->copy(other); } /* -------------------------------------------------------------------------- */ template auto ElementTypeMapArray::operator=( const ElementTypeMapArray & other) -> ElementTypeMapArray & { if (this != &other) { AKANTU_DEBUG_WARNING("You are copying the ElementTypeMapArray " << this->id << " are you sure it is on purpose"); this->id = other.id + "_copy"; this->name = other.name + "_copy"; this->copy(other); } return *this; } /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray::alloc( UInt size, UInt nb_component, const SupportType & type, GhostType ghost_type, const T & default_value) { std::string ghost_id; if (ghost_type == _ghost) { ghost_id = ":ghost"; } auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) { auto id = this->id + ":" + std::to_string(type) + ghost_id; this->getData(ghost_type)[type] = std::make_unique>(size, nb_component, default_value, id); return *(this->getData(ghost_type)[type]); } AKANTU_DEBUG_INFO("The vector " << this->id << this->printType(type, ghost_type) << " already exists, it is resized instead of allocated."); auto && array = *(it->second); array.resize(size); return array; } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value) { this->alloc(size, nb_component, type, _not_ghost, default_value); this->alloc(size, nb_component, type, _ghost, default_value); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::free() { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { auto & data = this->getData(gt); data.clear(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::clear() { for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { vect.second->clear(); } } } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMapArray::empty() const { bool is_empty = true; for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { is_empty &= vect.second->empty(); if (not is_empty) { return false; } } } return is_empty; } /* -------------------------------------------------------------------------- */ template template inline void ElementTypeMapArray::set(const ST & value) { for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { vect.second->set(value); } } } /* -------------------------------------------------------------------------- */ template inline const Array & ElementTypeMapArray::operator()(const SupportType & type, GhostType ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) { AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this const ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class(\"" << this->id << "\")"); } return *(it->second); } /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray::operator()(const SupportType & type, GhostType ghost_type) { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) { AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class (\"" << this->id << "\")"); } return *(it->second); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::setArray( const SupportType & type, GhostType ghost_type, const Array & vect) { auto it = this->getData(ghost_type).find(type); if (AKANTU_DEBUG_TEST(dblWarning) && it != this->getData(ghost_type).end() && it->second != &vect) { AKANTU_DEBUG_WARNING( "The Array " << this->printType(type, ghost_type) << " is already registred, this call can lead to a memory leak."); } this->getData(ghost_type)[type] = &(const_cast &>(vect)); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::onElementsRemoved( const ElementTypeMapArray & new_numbering) { for (auto gt : ghost_types) { for (auto && type : new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) { auto support_type = convertType(type); if (this->exists(support_type, gt)) { const auto & renumbering = new_numbering(type, gt); if (renumbering.empty()) { continue; } auto & vect = this->operator()(support_type, gt); auto nb_component = vect.getNbComponent(); Array tmp(renumbering.size(), nb_component); UInt new_size = 0; for (UInt i = 0; i < vect.size(); ++i) { UInt new_i = renumbering(i); if (new_i != UInt(-1)) { std::copy_n(vect.storage() + i * nb_component, nb_component, tmp.storage() + new_i * nb_component); ++new_size; } } tmp.resize(new_size); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ template void ElementTypeMapArray::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; for (UInt g = _not_ghost; g <= _ghost; ++g) { auto gt = (GhostType)g; const DataMap & data = this->getData(gt); typename DataMap::const_iterator it; for (it = data.begin(); it != data.end(); ++it) { stream << space << space << ElementTypeMapArray::printType(it->first, gt) << " [" << std::endl; it->second->printself(stream, indent + 3); stream << space << space << " ]" << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /* SupportType Iterator */ /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek) : list_begin(list_begin), list_end(list_end), dim(dim), kind(ek) {} /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( const type_iterator & it) : list_begin(it.list_begin), list_end(it.list_end), dim(it.dim), kind(it.kind) {} /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator::operator=( const type_iterator & it) { if (this != &it) { list_begin = it.list_begin; list_end = it.list_end; dim = it.dim; kind = it.kind; } return *this; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() const { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator::operator++() { ++list_begin; while ((list_begin != list_end) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(list_begin->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(list_begin->first))))) { ++list_begin; } return *this; } /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator ElementTypeMap::type_iterator::operator++(int) { type_iterator tmp(*this); operator++(); return tmp; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator::operator==( const type_iterator & other) const { return this->list_begin == other.list_begin; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator::operator!=( const type_iterator & other) const { return this->list_begin != other.list_begin; } /* -------------------------------------------------------------------------- */ template +template ::value> *> auto ElementTypeMap::ElementTypesIteratorHelper::begin() -> iterator { auto b = container.get().getData(ghost_type).begin(); auto e = container.get().getData(ghost_type).end(); // loop until the first valid type while ((b != e) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(b->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(b->first))))) { ++b; } return iterator(b, e, dim, kind); } template +template ::value> *> auto ElementTypeMap::ElementTypesIteratorHelper::end() -> iterator { auto e = container.get().getData(ghost_type).end(); return iterator(e, e, dim, kind); } /* -------------------------------------------------------------------------- */ template auto ElementTypeMap::elementTypesImpl( UInt dim, GhostType ghost_type, ElementKind kind) const -> ElementTypesIteratorHelper { return ElementTypesIteratorHelper(*this, dim, ghost_type, kind); } /* -------------------------------------------------------------------------- */ template template auto ElementTypeMap::elementTypesImpl( const use_named_args_t & unused, pack &&... _pack) const -> ElementTypesIteratorHelper { return ElementTypesIteratorHelper(*this, unused, _pack...); } /* -------------------------------------------------------------------------- */ template inline auto ElementTypeMap::firstType( UInt dim, GhostType ghost_type, ElementKind kind) const -> type_iterator { return elementTypes(dim, ghost_type, kind).begin(); } /* -------------------------------------------------------------------------- */ template inline auto ElementTypeMap::lastType( UInt dim, GhostType ghost_type, ElementKind kind) const -> type_iterator { typename DataMap::const_iterator e; e = getData(ghost_type).end(); return typename ElementTypeMap::type_iterator(e, e, dim, kind); } /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const ElementTypeMap & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ class ElementTypeMapArrayInitializer { protected: using CompFunc = std::function; public: ElementTypeMapArrayInitializer(const CompFunc & comp_func, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) : comp_func(comp_func), spatial_dimension(spatial_dimension), ghost_type(ghost_type), element_kind(element_kind) {} GhostType ghostType() const { return ghost_type; } virtual UInt nbComponent(ElementType type) const { return comp_func(type, ghostType()); } virtual bool isNodal() const { return false; } protected: CompFunc comp_func; UInt spatial_dimension; GhostType ghost_type; ElementKind element_kind; }; /* -------------------------------------------------------------------------- */ class MeshElementTypeMapArrayInitializer : public ElementTypeMapArrayInitializer { using CompFunc = ElementTypeMapArrayInitializer::CompFunc; public: MeshElementTypeMapArrayInitializer( const Mesh & mesh, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, bool with_nb_element = false, bool with_nb_nodes_per_element = false, const ElementTypeMapArray * filter = nullptr) : MeshElementTypeMapArrayInitializer( mesh, [nb_component](ElementType /*unused*/, GhostType /*unused*/) -> UInt { return nb_component; }, spatial_dimension, ghost_type, element_kind, with_nb_element, with_nb_nodes_per_element, filter) {} MeshElementTypeMapArrayInitializer( const Mesh & mesh, const CompFunc & comp_func, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, bool with_nb_element = false, bool with_nb_nodes_per_element = false, const ElementTypeMapArray * filter = nullptr) : ElementTypeMapArrayInitializer(comp_func, spatial_dimension, ghost_type, element_kind), mesh(mesh), with_nb_element(with_nb_element), with_nb_nodes_per_element(with_nb_nodes_per_element), filter(filter) {} decltype(auto) elementTypes() const { if (filter != nullptr) { return filter->elementTypes(this->spatial_dimension, this->ghost_type, this->element_kind); } return mesh.elementTypes(this->spatial_dimension, this->ghost_type, this->element_kind); } virtual UInt size(ElementType type) const { if (with_nb_element) { if (filter != nullptr) { return (*filter)(type, this->ghost_type).size(); } return mesh.getNbElement(type, this->ghost_type); } return 0; } UInt nbComponent(ElementType type) const override { UInt res = ElementTypeMapArrayInitializer::nbComponent(type); if (with_nb_nodes_per_element) { return (res * Mesh::getNbNodesPerElement(type)); } return res; } bool isNodal() const override { return with_nb_nodes_per_element; } protected: const Mesh & mesh; bool with_nb_element{false}; bool with_nb_nodes_per_element{false}; const ElementTypeMapArray * filter{nullptr}; }; /* -------------------------------------------------------------------------- */ class FEEngineElementTypeMapArrayInitializer : public MeshElementTypeMapArrayInitializer { public: FEEngineElementTypeMapArrayInitializer( const FEEngine & fe_engine, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined); FEEngineElementTypeMapArrayInitializer( const FEEngine & fe_engine, const ElementTypeMapArrayInitializer::CompFunc & nb_component, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined); UInt size(ElementType type) const override; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; ElementTypesIteratorHelper elementTypes() const; protected: const FEEngine & fe_engine; }; /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const Func & f, const T & default_value, bool do_not_default) { this->is_nodal = f.isNodal(); auto ghost_type = f.ghostType(); for (auto & type : f.elementTypes()) { if (not this->exists(type, ghost_type)) { if (do_not_default) { auto & array = this->alloc(0, f.nbComponent(type), type, ghost_type); array.resize(f.size(type)); } else { this->alloc(f.size(type), f.nbComponent(type), type, ghost_type, default_value); } } else { auto & array = this->operator()(type, ghost_type); if (not do_not_default) { array.resize(f.size(type), default_value); } else { array.resize(f.size(type)); } } } } /* -------------------------------------------------------------------------- */ /** * All parameters are named optionals * \param _nb_component a functor giving the number of components per * (ElementType, GhostType) pair or a scalar giving a unique number of * components * regardless of type * \param _spatial_dimension a filter for the elements of a specific dimension * \param _element_kind filter with element kind (_ek_regular, _ek_structural, * ...) * \param _with_nb_element allocate the arrays with the number of elements for * each * type in the mesh * \param _with_nb_nodes_per_element multiply the number of components by the * number of nodes per element * \param _default_value default inital value * \param _do_not_default do not initialize the allocated arrays * \param _ghost_type filter a type of ghost */ template template void ElementTypeMapArray::initialize(const Mesh & mesh, pack &&... _pack) { GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _casper); bool all_ghost_types = OPTIONAL_NAMED_ARG(all_ghost_types, requested_ghost_type == _casper); for (GhostType ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) { continue; } auto functor = MeshElementTypeMapArrayInitializer( mesh, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, mesh.getSpatialDimension()), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined), OPTIONAL_NAMED_ARG(with_nb_element, false), OPTIONAL_NAMED_ARG(with_nb_nodes_per_element, false), OPTIONAL_NAMED_ARG(element_filter, nullptr)); this->initialize(functor, OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ /** * All parameters are named optionals * \param _nb_component a functor giving the number of components per * (ElementType, GhostType) pair or a scalar giving a unique number of * components * regardless of type * \param _spatial_dimension a filter for the elements of a specific dimension * \param _element_kind filter with element kind (_ek_regular, _ek_structural, * ...) * \param _default_value default inital value * \param _do_not_default do not initialize the allocated arrays * \param _ghost_type filter a specific ghost type * \param _all_ghost_types get all ghost types */ template template void ElementTypeMapArray::initialize(const FEEngine & fe_engine, pack &&... _pack) { GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _casper); bool all_ghost_types = OPTIONAL_NAMED_ARG(all_ghost_types, requested_ghost_type == _casper); for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) { continue; } auto functor = FEEngineElementTypeMapArrayInitializer( fe_engine, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, UInt(-2)), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined)); this->initialize(functor, OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ template inline T & ElementTypeMapArray::operator()(const Element & element, UInt component) { return this->operator()(element.type, element.ghost_type)(element.element, component); } /* -------------------------------------------------------------------------- */ template inline const T & ElementTypeMapArray::operator()(const Element & element, UInt component) const { return this->operator()(element.type, element.ghost_type)(element.element, component); } /* -------------------------------------------------------------------------- */ template inline decltype(auto) ElementTypeMapArray::get(const Element & element) { auto & array = operator()(element.type, element.ghost_type); auto it = array.begin(array.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ template inline decltype(auto) ElementTypeMapArray::get(const Element & element) const { auto & array = operator()(element.type, element.ghost_type); auto it = array.begin(array.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ template UInt ElementTypeMapArray::sizeImpl(UInt spatial_dimension, GhostType ghost_type, ElementKind kind) const { UInt size = 0; for (auto && type : this->elementTypes(spatial_dimension, ghost_type, kind)) { size += this->operator()(type, ghost_type).size(); } return size; } /* -------------------------------------------------------------------------- */ template template UInt ElementTypeMapArray::size(pack &&... _pack) const { UInt size = 0; GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _casper); bool all_ghost_types = OPTIONAL_NAMED_ARG(all_ghost_types, requested_ghost_type == _casper); for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) { continue; } size += sizeImpl(OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined)); } return size; } } // namespace akantu #endif /* AKANTU_ELEMENT_TYPE_MAP_TMPL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc index 710c2cf1b..ecba9e81c 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc @@ -1,558 +1,558 @@ /** * @file material_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Thu Jan 14 2021 * * @brief Specialization of the material class for cohesive elements * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "aka_random_generator.hh" #include "dof_synchronizer.hh" #include "fe_engine_template.hh" #include "integrator_gauss.hh" #include "shape_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model, id), facet_filter("facet_filter", id), fem_cohesive( model.getFEEngineClass("CohesiveFEEngine")), reversible_energy("reversible_energy", *this), total_energy("total_energy", *this), opening("opening", *this), tractions("tractions", *this), contact_tractions("contact_tractions", *this), contact_opening("contact_opening", *this), delta_max("delta_max", *this), use_previous_delta_max(false), use_previous_opening(false), damage("damage", *this), sigma_c("sigma_c", *this), normal(0, spatial_dimension, "normal") { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable, "Critical stress"); this->registerParam("delta_c", delta_c, Real(0.), _pat_parsable | _pat_readable, "Critical displacement"); this->element_filter.initialize(this->model->getMesh(), _spatial_dimension = spatial_dimension, _element_kind = _ek_cohesive); // this->model->getMesh().initElementTypeMapArray( // this->element_filter, 1, spatial_dimension, false, _ek_cohesive); if (this->model->getIsExtrinsic()) { this->facet_filter.initialize(this->model->getMeshFacets(), _spatial_dimension = spatial_dimension - 1, _element_kind = _ek_regular); } // this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1, // spatial_dimension - // 1); this->reversible_energy.initialize(1); this->total_energy.initialize(1); this->tractions.initialize(spatial_dimension); this->tractions.initializeHistory(); this->contact_tractions.initialize(spatial_dimension); this->contact_opening.initialize(spatial_dimension); this->opening.initialize(spatial_dimension); this->opening.initializeHistory(); this->delta_max.initialize(1); this->damage.initialize(1); if (this->model->getIsExtrinsic()) { this->sigma_c.initialize(1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() = default; /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); if (this->use_previous_delta_max) { this->delta_max.initializeHistory(); } if (this->use_previous_opening) { this->opening.initializeHistory(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & internal_force = const_cast &>(model->getInternalForce()); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } const auto & shapes = fem_cohesive.getShapes(type, ghost_type); auto & traction = tractions(type, ghost_type); auto & contact_traction = contact_tractions(type, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); /// compute @f$t_i N_a@f$ auto * traction_cpy = new Array(nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); auto traction_it = traction.begin(spatial_dimension, 1); auto contact_traction_it = contact_traction.begin(spatial_dimension, 1); auto shapes_filtered_begin = shapes.begin(1, size_of_shapes); auto traction_cpy_it = traction_cpy->begin(spatial_dimension, size_of_shapes); Matrix traction_tmp(spatial_dimension, 1); for (UInt el = 0; el < nb_element; ++el) { UInt current_quad = elem_filter(el) * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++contact_traction_it, ++current_quad, ++traction_cpy_it) { const Matrix & shapes_filtered = shapes_filtered_begin[current_quad]; traction_tmp.copy(*traction_it); traction_tmp += *contact_traction_it; traction_cpy_it->mul(traction_tmp, shapes_filtered); } } /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ auto * partial_int_t_N = new Array( nb_element, spatial_dimension * size_of_shapes, "int_t_N"); fem_cohesive.integrate(*traction_cpy, *partial_int_t_N, spatial_dimension * size_of_shapes, type, ghost_type, elem_filter); delete traction_cpy; auto * int_t_N = new Array( nb_element, 2 * spatial_dimension * size_of_shapes, "int_t_N"); Real * int_t_N_val = int_t_N->storage(); Real * partial_int_t_N_val = partial_int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension, int_t_N_val); std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension, int_t_N_val + size_of_shapes * spatial_dimension); for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n) { int_t_N_val[n] *= -1.; } int_t_N_val += nb_nodes_per_element * spatial_dimension; partial_int_t_N_val += size_of_shapes * spatial_dimension; } delete partial_int_t_N; /// assemble model->getDOFManager().assembleElementalArrayLocalArray( *int_t_N, internal_force, type, ghost_type, 1, elem_filter); delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & shapes = fem_cohesive.getShapes(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0U) { continue; } UInt size_of_shapes = shapes.getNbComponent(); auto * shapes_filtered = new Array(nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_filtered_val = shapes_filtered->storage(); UInt * elem_filter_val = elem_filter.storage(); for (UInt el = 0; el < nb_element; ++el) { auto * shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } Matrix A(spatial_dimension * size_of_shapes, spatial_dimension * nb_nodes_per_element); for (UInt i = 0; i < spatial_dimension * size_of_shapes; ++i) { A(i, i) = 1; A(i, i + spatial_dimension * size_of_shapes) = -1; } /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} /// @f$ auto * tangent_stiffness_matrix = new Array( nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix"); // Array * normal = new Array(nb_element * // nb_quadrature_points, spatial_dimension, "normal"); normal.resize(nb_quadrature_points); computeNormal(model->getCurrentPosition(), normal, type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ // computeOpening(model->getDisplacement(), opening(type, ghost_type), type, // ghost_type); tangent_stiffness_matrix->zero(); computeTangentTraction(type, *tangent_stiffness_matrix, normal, ghost_type); // delete normal; UInt size_at_nt_d_n_a = spatial_dimension * nb_nodes_per_element * spatial_dimension * nb_nodes_per_element; auto * at_nt_d_n_a = new Array(nb_element * nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Array::iterator> shapes_filt_it = shapes_filtered->begin(size_of_shapes); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Array::matrix_iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Array::matrix_iterator At_Nt_D_N_A_end = at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Matrix N(spatial_dimension, spatial_dimension * size_of_shapes); Matrix N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); Matrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for (; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { N.zero(); /** * store the shapes in voigt notations matrix @f$\mathbf{N} = * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$ **/ for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt n = 0; n < size_of_shapes; ++n) { N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n); } } /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} *{\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ N_A.mul(N, A); D_N_A.mul(*D_it, N_A); (*At_Nt_D_N_A_it).mul(D_N_A, N_A); } delete tangent_stiffness_matrix; delete shapes_filtered; auto * K_e = new Array(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive.integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, type, ghost_type, elem_filter); delete at_nt_d_n_a; model->getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _unsymmetric, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); for (const auto & type : element_filter.elementTypes( spatial_dimension, ghost_type, _ek_cohesive)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } UInt nb_quadrature_points = nb_element * fem_cohesive.getNbIntegrationPoints(type, ghost_type); normal.resize(nb_quadrature_points); /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normal, type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(type, ghost_type), type, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normal, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); normal.zero(); #define COMPUTE_NORMAL(type) \ fem_cohesive.getShapeFunctions() \ .computeNormalsOnIntegrationPoints( \ position, normal, ghost_type, element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); #define COMPUTE_OPENING(type) \ fem_cohesive.getShapeFunctions() \ - .interpolateOnIntegrationPoints( \ + .interpolateOnIntegrationPoints( \ displacement, opening, spatial_dimension, ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::updateEnergies(ElementType type) { AKANTU_DEBUG_IN(); if (Mesh::getKind(type) != _ek_cohesive) { return; } Vector b(spatial_dimension); Vector h(spatial_dimension); auto erev = reversible_energy(type).begin(); auto etot = total_energy(type).begin(); auto traction_it = tractions(type).begin(spatial_dimension); auto traction_old_it = tractions.previous(type).begin(spatial_dimension); auto opening_it = opening(type).begin(spatial_dimension); auto opening_old_it = opening.previous(type).begin(spatial_dimension); auto traction_end = tractions(type).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration b = *opening_it; b -= *opening_old_it; h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } /// update old values AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate reversible energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { erev += fem_cohesive.integrate(reversible_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate dissipated energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { Array dissipated_energy(total_energy(type, _not_ghost)); dissipated_energy -= reversible_energy(type, _not_ghost); edis += fem_cohesive.integrate(dissipated_energy, type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getContactEnergy() { AKANTU_DEBUG_IN(); Real econ = 0.; /// integrate contact energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { auto & el_filter = element_filter(type, _not_ghost); UInt nb_quad_per_el = fem_cohesive.getNbIntegrationPoints(type, _not_ghost); UInt nb_quad_points = el_filter.size() * nb_quad_per_el; Array contact_energy(nb_quad_points); auto contact_traction_it = contact_tractions(type, _not_ghost).begin(spatial_dimension); auto contact_opening_it = contact_opening(type, _not_ghost).begin(spatial_dimension); /// loop on each quadrature point for (UInt q = 0; q < nb_quad_points; ++contact_traction_it, ++contact_opening_it, ++q) { contact_energy(q) = .5 * contact_traction_it->dot(*contact_opening_it); } econ += fem_cohesive.integrate(contact_energy, type, _not_ghost, el_filter); } AKANTU_DEBUG_OUT(); return econ; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(const std::string & type) { if (type == "reversible") { return getReversibleEnergy(); } if (type == "dissipated") { return getDissipatedEnergy(); } if (type == "cohesive contact") { return getContactEnergy(); } return 0.; } /* -------------------------------------------------------------------------- */ } // namespace akantu