diff --git a/packages/embedded.cmake b/packages/embedded.cmake index b223dfbb5..a4599d8ee 100644 --- a/packages/embedded.cmake +++ b/packages/embedded.cmake @@ -1,60 +1,59 @@ #=============================================================================== # @file embedded.cmake # # @author Lucas Frerot # @author Nicolas Richart # # @date creation: Tue Oct 16 2012 # @date last modification: Mon Dec 14 2015 # # @brief package descrition for embedded model use # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 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 . # #=============================================================================== package_declare(embedded DESCRIPTION "Add support for the embedded solid mechanics model" DEPENDS CGAL) package_declare_sources(embedded model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.cc model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.hh model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.hh model/solid_mechanics/materials/material_embedded/embedded_internal_field.hh model/solid_mechanics/materials/material_embedded/material_embedded_includes.hh - model/solid_mechanics/materials/material_embedded/material_reinforcement.cc + model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh model/solid_mechanics/materials/material_embedded/material_reinforcement.hh - model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc model/solid_mechanics/materials/material_embedded/material_reinforcement_template.hh model/solid_mechanics/materials/material_embedded/material_reinforcement_template_tmpl.hh ) package_declare_material_infos(embedded LIST AKANTU_EMBEDDED_MATERIAL_LIST INCLUDE material_embedded_includes.hh ) package_declare_documentation(embedded "This package allows the use of the embedded model in solid mechanics. This package depends on the CGAL package." ) #add_example(embedded "Example on how to run embedded model simulation" PACKAGE embedded) diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh index c6fc7a798..cab03849f 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,437 +1,447 @@ /** * @file element_type_map.hh * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Fri Oct 02 2015 * * @brief storage class by element type * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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_memory.hh" #include "aka_named_argument.hh" #include "element.hh" /* -------------------------------------------------------------------------- */ #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(with_nb_element); DECLARE_NAMED_ARGUMENT(with_nb_nodes_per_element); DECLARE_NAMED_ARGUMENT(spatial_dimension); DECLARE_NAMED_ARGUMENT(do_not_default); } // 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, const 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, const 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, const 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, const 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. */ inline Stored & operator()(const Stored & insert, const SupportType & type, const GhostType & ghost_type = _not_ghost); + /// Insert into map by move + inline Stored & operator()(Stored && insert, const SupportType & type, + const GhostType & ghost_type = _not_ghost); +protected: + /// Insert function with universal reference + template + inline Stored & insert(U && insertee, const SupportType & type, + const 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; class type_iterator : private std::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 & other); 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, pack &&... _pack) : ElementTypesIteratorHelper( container, OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), OPTIONAL_NAMED_ARG(ghost_type, _not_ghost), OPTIONAL_NAMED_ARG(element_kind, _ek_regular)) {} ElementTypesIteratorHelper(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(ElementTypesIteratorHelper &&) = default; iterator begin() { return container.get().firstType(dim, ghost_type, kind); } iterator end() { return container.get().lastType(dim, ghost_type, kind); } 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_regular) const; template ElementTypesIteratorHelper elementTypesImpl(const use_named_args_t & /*unused*/, pack &&... _pack) const; public: 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)...); } public: /*! 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*/ 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 */ inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; protected: /*! 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 Memory { public: using 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 void operator=(const ElementTypeMapArray &) = delete; ElementTypeMapArray(const ElementTypeMapArray &) = delete; /*! Constructor * @param id optional: identifier (string) * @param parent_id optional: parent identifier. for organizational purposes * only * @param memory_id optional: choose a specific memory, defaults to memory 0 */ ElementTypeMapArray(const ID & id = "by_element_type_array", const ID & parent_id = "no_parent", const MemoryID & memory_id = 0) : parent(), Memory(parent_id + ":" + id, memory_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 * @return a reference to the allocated array */ inline Array & alloc(UInt size, UInt nb_component, const SupportType & type, const 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*/ 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, const 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); /* 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, const 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, const GhostType & ghost_type, const Array & vect); /*! frees all memory related to the data*/ inline void free(); /*! set all values in the ElementTypeMap to zero*/ inline void clear(); /*! 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; } ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const { ElementTypeMap nb_components; 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); /// name of the elment type map: e.g. connectivity, grad_u ID name; }; /// 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; /// Map of data of type UInt stored in a mesh using UIntDataMap = std::map *>; using ElementTypeMapUIntDataMap = ElementTypeMap; } // 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 af2a331d6..5888081e1 100644 --- a/src/mesh/element_type_map_tmpl.hh +++ b/src/mesh/element_type_map_tmpl.hh @@ -1,661 +1,676 @@ /** * @file element_type_map_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Fri Oct 02 2015 * * @brief implementation of template functions of the ElementTypeMap and * ElementTypeMapArray classes * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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" /* -------------------------------------------------------------------------- */ #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, const GhostType & ghost_type) { std::stringstream sstr; sstr << "(" << ghost_type << ":" << type << ")"; return sstr.str(); } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::exists( const SupportType & type, const 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, const 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, const GhostType & ghost_type) { return this->getData(ghost_type)[type]; } /* -------------------------------------------------------------------------- */ template inline Stored & ElementTypeMap:: -operator()(const Stored & insert, const SupportType & type, +operator()(const Stored & insertee, const SupportType & type, const GhostType & ghost_type) { + return insert(insertee, type, ghost_type); +} + +template +inline Stored & ElementTypeMap:: +operator()(Stored && insertee, const SupportType & type, + const GhostType & ghost_type) { + return insert(std::move(insertee), type, ghost_type); +} + +/* -------------------------------------------------------------------------- */ +template +template +inline Stored & ElementTypeMap::insert( + U && insertee, const SupportType & type, const 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::pair(type, insert)); + 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; for (Int i = 0; i < indent; i++, space += 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 inline Array & ElementTypeMapArray::alloc( UInt size, UInt nb_component, const SupportType & type, const GhostType & ghost_type, const T & default_value) { std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; Array * tmp; auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) { std::stringstream sstr; sstr << this->id << ":" << type << ghost_id; tmp = &(Memory::alloc(sstr.str(), size, nb_component, default_value)); std::stringstream sstrg; sstrg << ghost_type; // tmp->setTag(sstrg.str()); this->getData(ghost_type)[type] = tmp; } else { AKANTU_DEBUG_INFO( "The vector " << this->id << this->printType(type, ghost_type) << " already exists, it is resized instead of allocated."); tmp = it->second; it->second->resize(size); } return *tmp; } /* -------------------------------------------------------------------------- */ 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); for (auto & pair : data) { dealloc(pair.second->getID()); } 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 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, const 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, const 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, const 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.size() == 0) 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)) { memcpy(tmp.storage() + new_i * nb_component, vect.storage() + i * nb_component, nb_component * sizeof(T)); ++new_size; } } tmp.resize(new_size); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ template void ElementTypeMapArray::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += 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 typename ElementTypeMap::ElementTypesIteratorHelper ElementTypeMap::elementTypesImpl(UInt dim, GhostType ghost_type, ElementKind kind) const { return ElementTypesIteratorHelper(*this, dim, ghost_type, kind); } /* -------------------------------------------------------------------------- */ template template typename ElementTypeMap::ElementTypesIteratorHelper ElementTypeMap::elementTypesImpl( const use_named_args_t & unused, pack &&... _pack) const { return ElementTypesIteratorHelper(*this, unused, _pack...); } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator ElementTypeMap::firstType(UInt dim, GhostType ghost_type, ElementKind kind) const { typename DataMap::const_iterator b, e; b = getData(ghost_type).begin(); e = 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 typename ElementTypeMap::type_iterator(b, e, dim, kind); } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator ElementTypeMap::lastType(UInt dim, GhostType ghost_type, ElementKind kind) const { 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 ElementTypeMapArrayInializer { public: ElementTypeMapArrayInializer(UInt spatial_dimension = _all_dimensions, UInt nb_component = 1, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular) : spatial_dimension(spatial_dimension), nb_component(nb_component), ghost_type(ghost_type), element_kind(element_kind) {} const GhostType & ghostType() const { return ghost_type; } protected: UInt spatial_dimension; UInt nb_component; GhostType ghost_type; ElementKind element_kind; }; /* -------------------------------------------------------------------------- */ class MeshElementTypeMapArrayInializer : public ElementTypeMapArrayInializer { public: MeshElementTypeMapArrayInializer( const Mesh & mesh, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular, bool with_nb_element = false, bool with_nb_nodes_per_element = false) : ElementTypeMapArrayInializer(spatial_dimension, nb_component, ghost_type, element_kind), mesh(mesh), with_nb_element(with_nb_element), with_nb_nodes_per_element(with_nb_nodes_per_element) {} decltype(auto) elementTypes() const { return mesh.elementTypes(this->spatial_dimension, this->ghost_type, this->element_kind); } virtual UInt size(const ElementType & type) const { if (with_nb_element) return mesh.getNbElement(type, this->ghost_type); return 0; } UInt nbComponent(const ElementType & type) const { if (with_nb_nodes_per_element) return (this->nb_component * mesh.getNbNodesPerElement(type)); return this->nb_component; } protected: const Mesh & mesh; bool with_nb_element; bool with_nb_nodes_per_element; }; /* -------------------------------------------------------------------------- */ class FEEngineElementTypeMapArrayInializer : public MeshElementTypeMapArrayInializer { public: FEEngineElementTypeMapArrayInializer( const FEEngine & fe_engine, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular); UInt size(const 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) { for (auto & type : f.elementTypes()) { if (not this->exists(type, f.ghostType())) if (do_not_default) { auto & array = this->alloc(0, f.nbComponent(type), type, f.ghostType()); array.resize(f.size(type)); } else { this->alloc(f.size(type), f.nbComponent(type), type, f.ghostType(), default_value); } else { auto & array = this->operator()(type, f.ghostType()); if (not do_not_default) array.resize(f.size(type), default_value); else array.resize(f.size(type)); } } } /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const Mesh & mesh, pack &&... _pack) { GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _casper); 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; this->initialize( MeshElementTypeMapArrayInializer( mesh, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, mesh.getSpatialDimension()), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_regular), OPTIONAL_NAMED_ARG(with_nb_element, false), OPTIONAL_NAMED_ARG(with_nb_nodes_per_element, false)), OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const FEEngine & fe_engine, pack &&... _pack) { bool all_ghost_types = OPTIONAL_NAMED_ARG(all_ghost_types, true); GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _not_ghost); for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; this->initialize(FEEngineElementTypeMapArrayInializer( fe_engine, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, UInt(-2)), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_regular)), 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); } } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh index 09925013b..d35c099e0 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh @@ -1,192 +1,206 @@ /** * @file material_reinforcement.hh * * @author Lucas Frerot * * @date creation: Fri Mar 13 2015 * @date last modification: Tue Nov 24 2015 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 2015 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_REINFORCEMENT_HH__ #define __AKANTU_MATERIAL_REINFORCEMENT_HH__ #include "aka_common.hh" #include "material.hh" #include "embedded_interface_model.hh" #include "embedded_internal_field.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /** * @brief Material used to represent embedded reinforcements * * This class is used for computing the reinforcement stiffness matrix * along with the reinforcement residual. Room is made for constitutive law, * but actual use of contitutive laws is made in MaterialReinforcementTemplate. * * Be careful with the dimensions in this class : * - this->spatial_dimension is always 1 * - the template parameter dim is the dimension of the problem */ -template class MaterialReinforcement : virtual public Material { + +template class MaterialReinforcement : public Mat { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Constructor - MaterialReinforcement(SolidMechanicsModel & model, UInt spatial_dimension, - const Mesh & mesh, FEEngine & fe_engine, - const ID & id = ""); + MaterialReinforcement(EmbeddedInterfaceModel & model, const ID & id = ""); /// Destructor ~MaterialReinforcement() override; protected: - void initialize(SolidMechanicsModel & a_model); + void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Init the material void initMaterial() override; + /// Init the filter for background elements + void initBackgroundFilter(); + /// Init the background shape derivatives void initBackgroundShapeDerivatives(); /// Init the cosine matrices void initDirectingCosines(); /// Assemble stiffness matrix void assembleStiffnessMatrix(GhostType ghost_type) override; /// Compute all the stresses ! void computeAllStresses(GhostType ghost_type) override; /// Compute energy Real getEnergy(const std::string & id) override; /// Assemble the residual of one type of element (typically _segment_2) void assembleInternalForces(GhostType ghost_type) override; // virtual ElementTypeMap getInternalDataPerElem(const ID & field_name, // const ElementKind & // kind, // const ID & // fe_engine_id) const; // /// Reimplementation of Material's function to accomodate for interface // mesh // virtual void flattenInternal(const std::string & field_id, // ElementTypeMapArray & internal_flat, // const GhostType ghost_type = _not_ghost, // ElementKind element_kind = _ek_not_defined) // const; /* ------------------------------------------------------------------------ */ /* Protected methods */ /* ------------------------------------------------------------------------ */ protected: /// Allocate the background shape derivatives void allocBackgroundShapeDerivatives(); /// Compute the directing cosines matrix for one element type void computeDirectingCosines(const ElementType & type, GhostType ghost_type); /// Compute the directing cosines matrix on quadrature points. inline void computeDirectingCosinesOnQuad(const Matrix & nodes, Matrix & cosines); /// Assemble the stiffness matrix for an element type (typically _segment_2) void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// Assemble the stiffness matrix for background & interface types void assembleStiffnessMatrixInterface(const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type); /// Compute the background shape derivatives for a type void computeBackgroundShapeDerivatives(const ElementType & type, GhostType ghost_type); + /// Compute the background shape derivatives for a type pair + void computeBackgroundShapeDerivatives(const ElementType & interface_type, + const ElementType & bg_type, + GhostType ghost_type, + const Array & filter); + /// Filter elements crossed by interface of a type void filterInterfaceBackgroundElements(Array & filter, const ElementType & type, const ElementType & interface_type, GhostType ghost_type); /// Assemble the residual of one type of element (typically _segment_2) void assembleInternalForces(const ElementType & type, GhostType ghost_type); /// Assemble the residual for a pair of elements void assembleInternalForcesInterface(const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type); // TODO figure out why voigt size is 4 in 2D inline void stressTensorToVoigtVector(const Matrix & tensor, Vector & vector); inline void strainTensorToVoigtVector(const Matrix & tensor, Vector & vector); /// Compute gradu on the interface quadrature points - virtual void computeGradU(const ElementType & type, GhostType ghost_type); + void computeGradU(const ElementType & type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Embedded model - EmbeddedInterfaceModel * model; + EmbeddedInterfaceModel & emodel; /// Stress in the reinforcement InternalField stress_embedded; /// Gradu of concrete on reinforcement InternalField gradu_embedded; /// C matrix on quad InternalField directing_cosines; /// Prestress on quad InternalField pre_stress; /// Cross-sectional area Real area; + template + using CrossMap = ElementTypeMap>>; + /// Background mesh shape derivatives - ElementTypeMap *> shape_derivatives; -}; + CrossMap shape_derivatives; -#include "material_reinforcement_inline_impl.cc" + /// Background element filter (contains segment-bg pairs) + CrossMap background_filter; +}; } // akantu +#include "material_reinforcement_tmpl.hh" + #endif // __AKANTU_MATERIAL_REINFORCEMENT_HH__ diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc deleted file mode 100644 index cf23faf42..000000000 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc +++ /dev/null @@ -1,124 +0,0 @@ -/** - * @file material_reinforcement_inline_impl.cc - * - * @author Lucas Frerot - * - * @date creation: Fri Jun 18 2010 - * @date last modification: Tue Jul 14 2015 - * - * @brief Reinforcement material - * - * @section LICENSE - * - * Copyright (©) 2010-2012, 2014, 2015 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 . - * - */ - -/* -------------------------------------------------------------------------- */ - -/** - * The structure of the directing cosines matrix is : - * \f{eqnarray*}{ - * C_{1,\cdot} & = & (l^2, m^2, n^2, mn, ln, lm) \\ - * C_{i,j} & = & 0 - * \f} - * - * with : - * \f[ - * (l, m, n) = \frac{1}{\|\frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}\|} \cdot \frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s} - * \f] - */ -template -inline void MaterialReinforcement::computeDirectingCosinesOnQuad(const Matrix & nodes, - Matrix & cosines) { - AKANTU_DEBUG_IN(); - - AKANTU_DEBUG_ASSERT(nodes.cols() == 2, "Higher order reinforcement elements not implemented"); - - const Vector & a = nodes(0), b = nodes(1); - Vector delta = b - a; - - cosines.clear(); - - Real sq_length = 0.; - for (UInt i = 0 ; i < dim ; i++) { - sq_length += Math::pow<2, Real>(delta(i)); - } - - if (dim == 2) { - cosines(0, 0) = Math::pow<2, Real>(delta(0)); // l^2 - cosines(0, 1) = Math::pow<2, Real>(delta(1)); // m^2 - cosines(0, 2) = delta(0) * delta(1); // lm - } else if (dim == 3) { - cosines(0, 0) = Math::pow<2, Real>(delta(0)); // l^2 - cosines(0, 1) = Math::pow<2, Real>(delta(1)); // m^2 - cosines(0, 2) = Math::pow<2, Real>(delta(2)); // n^2 - - cosines(0, 3) = delta(1) * delta(2); // mn - cosines(0, 4) = delta(0) * delta(2); // ln - cosines(0, 5) = delta(0) * delta(1); // lm - } - - cosines /= sq_length; - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ - -template -inline void MaterialReinforcement::stressTensorToVoigtVector(const Matrix & tensor, - Vector & vector) { - AKANTU_DEBUG_IN(); - - for (UInt i = 0; i < dim; i++) { - vector(i) = tensor(i, i); - } - - if (dim == 2) { - vector(2) = tensor(0, 1); - } else if (dim == 3) { - vector(3) = tensor(1, 2); - vector(4) = tensor(0, 2); - vector(5) = tensor(0, 1); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ - -template -inline void MaterialReinforcement::strainTensorToVoigtVector(const Matrix & tensor, - Vector & vector) { - AKANTU_DEBUG_IN(); - - for (UInt i = 0; i < dim; i++) { - vector(i) = tensor(i, i); - } - - if (dim == 2) { - vector(2) = 2 * tensor(0, 1); - } else if (dim == 3) { - vector(3) = 2 * tensor(1, 2); - vector(4) = 2 * tensor(0, 2); - vector(5) = 2 * tensor(0, 1); - } - - AKANTU_DEBUG_OUT(); -} diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh similarity index 51% rename from src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc rename to src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh index def0887e8..961fd7a82 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh @@ -1,766 +1,830 @@ /** - * @file material_reinforcement.cc + * @file material_reinforcement_tmpl.hh * * @author Lucas Frerot * - * @date creation: Wed Mar 25 2015 - * @date last modification: Tue Dec 08 2015 + * @date creation: Thu Feb 1 2018 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 2015 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_voigthelper.hh" #include "material_reinforcement.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -template -MaterialReinforcement::MaterialReinforcement(SolidMechanicsModel & model, - UInt /*spatial_dimension*/, - const Mesh & mesh, - FEEngine & fe_engine, - const ID & id) - : Material(model, dim, mesh, fe_engine, - id), // /!\ dim, not spatial_dimension ! - model(NULL), - stress_embedded("stress_embedded", *this, 1, fe_engine, +template +MaterialReinforcement::MaterialReinforcement( + EmbeddedInterfaceModel & model, const ID & id) + : Mat(model, model.getInterfaceMesh().getSpatialDimension(), + model.getInterfaceMesh(), + model.getFEEngine("EmbeddedInterfaceFEEngine"), id), + emodel(model), + stress_embedded("stress_embedded", *this, 1, + model.getFEEngine("EmbeddedInterfaceFEEngine"), this->element_filter), - gradu_embedded("gradu_embedded", *this, 1, fe_engine, + gradu_embedded("gradu_embedded", *this, 1, + model.getFEEngine("EmbeddedInterfaceFEEngine"), this->element_filter), - directing_cosines("directing_cosines", *this, 1, fe_engine, + directing_cosines("directing_cosines", *this, 1, + model.getFEEngine("EmbeddedInterfaceFEEngine"), this->element_filter), - pre_stress("pre_stress", *this, 1, fe_engine, this->element_filter), + pre_stress("pre_stress", *this, 1, + model.getFEEngine("EmbeddedInterfaceFEEEngine"), + this->element_filter), area(1.0), shape_derivatives() { AKANTU_DEBUG_IN(); - this->initialize(model); + this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::initialize(SolidMechanicsModel & a_model) { +template +void MaterialReinforcement::initialize() { AKANTU_DEBUG_IN(); - this->model = dynamic_cast(&a_model); - AKANTU_DEBUG_ASSERT(this->model != NULL, - "MaterialReinforcement needs an EmbeddedInterfaceModel"); this->registerParam("area", area, _pat_parsable | _pat_modifiable, "Reinforcement cross-sectional area"); this->registerParam("pre_stress", pre_stress, _pat_parsable | _pat_modifiable, "Uniform pre-stress"); this->unregisterInternal(this->stress); // Fool the AvgHomogenizingFunctor // stress.initialize(dim * dim); // Reallocate the element filter - this->element_filter.initialize(this->model->getInterfaceMesh(), - _spatial_dimension = 1); + this->element_filter.initialize(this->emodel.getInterfaceMesh(), + _spatial_dimension = 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template MaterialReinforcement::~MaterialReinforcement() { +template +MaterialReinforcement::~MaterialReinforcement() { AKANTU_DEBUG_IN(); - ElementTypeMap *>::type_iterator - it = shape_derivatives.firstType(), - end = shape_derivatives.lastType(); - - for (; it != end; ++it) { - delete shape_derivatives(*it, _not_ghost); - delete shape_derivatives(*it, _ghost); - } - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template void MaterialReinforcement::initMaterial() { +template +void MaterialReinforcement::initMaterial() { Material::initMaterial(); stress_embedded.initialize(dim * dim); gradu_embedded.initialize(dim * dim); pre_stress.initialize(1); /// We initialise the stuff that is not going to change during the simulation + this->initBackgroundFilter(); this->allocBackgroundShapeDerivatives(); this->initBackgroundShapeDerivatives(); this->initDirectingCosines(); } /* -------------------------------------------------------------------------- */ -/** - * Background shape derivatives need to be stored per background element - * types but also per embedded element type, which is why they are stored - * in an ElementTypeMap *>. The outer ElementTypeMap - * refers to the embedded types, and the inner refers to the background types. - */ -template -void MaterialReinforcement::allocBackgroundShapeDerivatives() { - AKANTU_DEBUG_IN(); +namespace detail { + class BackgroundFilterInitializer : public MeshElementTypeMapArrayInializer { + public: + BackgroundFilterInitializer(EmbeddedInterfaceModel & emodel, + const GhostType & ghost_type) + : MeshElementTypeMapArrayInializer( + emodel.getMesh(), emodel.getSpatialDimension(), 2, // pairs + ghost_type, _ek_regular) {} + + UInt size(const ElementType & /*bgtype*/) const override { + return 0; + } + }; +} - Mesh & interface_mesh = model->getInterfaceMesh(); - Mesh & mesh = model->getMesh(); +/* -------------------------------------------------------------------------- */ +/// Initialize the filter for background elements +template +void MaterialReinforcement::initBackgroundFilter() { + for (auto gt : ghost_types) { + for (auto && type : emodel.getInterfaceMesh().elementTypes(1, gt)) { + std::string shaped_id = "embedded_shape_derivatives"; + if (gt == _ghost) + shaped_id += ":ghost"; - ghost_type_t::iterator int_ghost_it = ghost_type_t::begin(); + auto filter = + std::make_unique>(shaped_id, this->name); + background_filter(std::move(filter), type, gt); - // Loop over interface ghosts - for (; int_ghost_it != ghost_type_t::end(); ++int_ghost_it) { - Mesh::type_iterator interface_type_it = - interface_mesh.firstType(1, *int_ghost_it); - Mesh::type_iterator interface_type_end = - interface_mesh.lastType(1, *int_ghost_it); + filter->initialize( + detail::BackgroundFilterInitializer(emodel, gt), 0, true); - // Loop over interface types - for (; interface_type_it != interface_type_end; ++interface_type_it) { - Mesh::type_iterator background_type_it = - mesh.firstType(dim, *int_ghost_it); - Mesh::type_iterator background_type_end = - mesh.lastType(dim, *int_ghost_it); + // Computing filters + for (auto && bg_type : filter->elementTypes(dim, gt)) { + filterInterfaceBackgroundElements((*filter)(bg_type), bg_type, type, + gt); + } + } + } +} - // Loop over background types - for (; background_type_it != background_type_end ; ++background_type_it) { - const ElementType & int_type = *interface_type_it; - const ElementType & back_type = *background_type_it; - const GhostType & int_ghost = *int_ghost_it; +/* -------------------------------------------------------------------------- */ +/// Construct a filter for a (interface_type, background_type) pair +template +void MaterialReinforcement::filterInterfaceBackgroundElements( + Array & filter, const ElementType & type, + const ElementType & interface_type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); - std::string shaped_id = "embedded_shape_derivatives"; + filter.resize(0); + + AKANTU_DEBUG_ASSERT(filter.getNbComponent() == 2, + "background filter should have 2 components"); + + Array & elements = + emodel.getInterfaceAssociatedElements(interface_type, ghost_type); + Array & elem_filter = this->element_filter(interface_type, ghost_type); - if (int_ghost == _ghost) shaped_id += ":ghost"; + for (auto & elem_id : elem_filter) { + Element & elem = elements(elem_id); + if (elem.type == type) + filter.push_back(Vector({elem_id, elem.element})); + } - ElementTypeMapArray * shaped_etma = new ElementTypeMapArray(shaped_id, this->name); + AKANTU_DEBUG_OUT(); +} - UInt nb_points = Mesh::getNbNodesPerElement(back_type); - UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine").getNbIntegrationPoints(int_type); - UInt nb_elements = element_filter(int_type, int_ghost).size(); - // Alloc the background ElementTypeMapArray - shaped_etma->alloc(nb_elements * nb_quad_points, - dim * nb_points, - back_type); +/* -------------------------------------------------------------------------- */ - // Insert the background ElementTypeMapArray in the interface - // ElementTypeMap - shape_derivatives(shaped_etma, int_type, int_ghost); +namespace detail { + class BackgroundShapeDInitializer : public ElementTypeMapArrayInializer { + public: + BackgroundShapeDInitializer(UInt spatial_dimension, + ElementTypeMapArray & filter, + const GhostType & ghost_type) + : ElementTypeMapArrayInializer(spatial_dimension, 0, ghost_type, + _ek_regular) { + // Conting how many background elements are affected by elements of + // interface_type + for (auto type : filter.elementTypes(this->spatial_dimension)) { + // Inserting size + array_size_per_bg_type(filter(type).size(), type, this->ghost_type); } } + + auto elementTypes() const -> decltype(auto){ + return array_size_per_bg_type.elementTypes(); + } + + UInt size(const ElementType & bgtype) const { + return array_size_per_bg_type(bgtype, this->ghost_type); + } + + UInt nbComponent(const ElementType & bgtype) const { + return ShapeFunctions::getShapeDerivativesSize(bgtype); + } + + protected: + ElementTypeMap array_size_per_bg_type; + }; +} + +/** + * Background shape derivatives need to be stored per background element + * types but also per embedded element type, which is why they are stored + * in an ElementTypeMap *>. The outer ElementTypeMap + * refers to the embedded types, and the inner refers to the background types. + */ +template +void MaterialReinforcement::allocBackgroundShapeDerivatives() { + AKANTU_DEBUG_IN(); + + for (auto gt : ghost_types) { + for (auto && type : emodel.getInterfaceMesh().elementTypes(1, gt)) { + std::string shaped_id = "embedded_shape_derivatives"; + if (gt == _ghost) + shaped_id += ":ghost"; + + auto shaped_etma = + std::make_unique>(shaped_id, this->name); + shape_derivatives(std::move(shaped_etma), type, gt); + + shaped_etma->initialize( + detail::BackgroundShapeDInitializer(emodel.getSpatialDimension(), + *background_filter(type), gt), + 0, true); + } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::initBackgroundShapeDerivatives() { +template +void MaterialReinforcement::initBackgroundShapeDerivatives() { AKANTU_DEBUG_IN(); - Mesh & mesh = model->getMesh(); - - Mesh::type_iterator type_it = mesh.firstType(dim, _not_ghost); - Mesh::type_iterator type_end = mesh.lastType(dim, _not_ghost); - - for (; type_it != type_end; ++type_it) { - computeBackgroundShapeDerivatives(*type_it, _not_ghost); - // computeBackgroundShapeDerivatives(*type_it, _ghost); + for (auto interface_type : + this->element_filter.elementTypes(this->spatial_dimension)) { + for (auto type : background_filter(interface_type)->elementTypes(dim)) { + computeBackgroundShapeDerivatives(interface_type, type, _not_ghost, + this->element_filter(interface_type)); + } } AKANTU_DEBUG_OUT(); } + +/* -------------------------------------------------------------------------- */ +template +void MaterialReinforcement::computeBackgroundShapeDerivatives( + const ElementType & interface_type, const ElementType & bg_type, + GhostType ghost_type, const Array & filter) { + auto & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); + auto & engine = emodel.getFEEngine(); + auto & interface_mesh = emodel.getInterfaceMesh(); + + const auto ndof = dim * Mesh::getNbNodesPerElement(bg_type); + const auto nb_stress = + ShapeFunctions::getShapeDerivativesSize(bg_type) / ndof; + // const auto nb_strss = VoigtHelper::size; + const auto nb_quads_per_elem = + interface_engine.getNbIntegrationPoints(interface_type); + + Array quad_pos(0, dim, "interface_quad_pos"); + interface_engine.interpolateOnIntegrationPoints(interface_mesh.getNodes(), + quad_pos, dim, interface_type, + ghost_type, filter); + auto & background_shapesd = + (*shape_derivatives(interface_type, ghost_type))(bg_type, ghost_type); + auto & background_elements = + (*background_filter(interface_type, ghost_type))(bg_type, ghost_type); + + auto shapesd_begin = + background_shapesd.begin(nb_stress, ndof, nb_quads_per_elem); + auto quad_begin = quad_pos.begin(dim, Mesh::getNbNodesPerElement(bg_type)); + + for (const auto & elem : make_view(background_elements, 2)) { + for (UInt i = 0; i < nb_quads_per_elem; ++i) { + Matrix shapesd = Tensor3(shapesd_begin[elem(0)])(i); + Vector quads = Matrix(quad_begin[elem(0)])(i); + + engine.computeShapeDerivatives(quads, elem(1), bg_type, + shapesd, ghost_type); + } + } +} + /* -------------------------------------------------------------------------- */ -template void MaterialReinforcement::initDirectingCosines() { +template +void MaterialReinforcement::initDirectingCosines() { AKANTU_DEBUG_IN(); - Mesh & mesh = model->getInterfaceMesh(); + Mesh & mesh = emodel.getInterfaceMesh(); Mesh::type_iterator type_it = mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = mesh.lastType(1, _not_ghost); - const UInt voigt_size = getTangentStiffnessVoigtSize(dim); + const UInt voigt_size = VoigtHelper::size; directing_cosines.initialize(voigt_size * voigt_size); for (; type_it != type_end; ++type_it) { computeDirectingCosines(*type_it, _not_ghost); // computeDirectingCosines(*type_it, _ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::assembleStiffnessMatrix(GhostType ghost_type) { +template +void MaterialReinforcement::assembleStiffnessMatrix( + GhostType ghost_type) { AKANTU_DEBUG_IN(); - Mesh & interface_mesh = model->getInterfaceMesh(); + Mesh & interface_mesh = emodel.getInterfaceMesh(); Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); for (; type_it != type_end; ++type_it) { assembleStiffnessMatrix(*type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::assembleInternalForces(GhostType ghost_type) { +template +void MaterialReinforcement::assembleInternalForces( + GhostType ghost_type) { AKANTU_DEBUG_IN(); - Mesh & interface_mesh = model->getInterfaceMesh(); + Mesh & interface_mesh = emodel.getInterfaceMesh(); Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); for (; type_it != type_end; ++type_it) { this->assembleInternalForces(*type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::computeGradU(const ElementType & type, - GhostType ghost_type) { +template +void MaterialReinforcement::computeGradU(const ElementType & type, + GhostType ghost_type) { AKANTU_DEBUG_IN(); - Array & elem_filter = element_filter(type, ghost_type); + Array & elem_filter = this->element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); - UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine") + UInt nb_quad_points = emodel.getFEEngine("EmbeddedInterfaceFEEngine") .getNbIntegrationPoints(type); Array & gradu_vec = gradu_embedded(type, ghost_type); - Mesh::type_iterator back_it = model->getMesh().firstType(dim, ghost_type); - Mesh::type_iterator back_end = model->getMesh().lastType(dim, ghost_type); + Mesh::type_iterator back_it = emodel.getMesh().firstType(dim, ghost_type); + Mesh::type_iterator back_end = emodel.getMesh().lastType(dim, ghost_type); for (; back_it != back_end; ++back_it) { UInt nodes_per_background_e = Mesh::getNbNodesPerElement(*back_it); Array & shapesd = shape_derivatives(type, ghost_type)->operator()(*back_it, ghost_type); Array * background_filter = new Array(nb_element, 1, "background_filter"); filterInterfaceBackgroundElements(*background_filter, *back_it, type, ghost_type); - Array * disp_per_element = new Array(0, dim * nodes_per_background_e, "disp_elem"); - FEEngine::extractNodalToElementField(model->getMesh(), - model->getDisplacement(), - *disp_per_element, - *back_it, ghost_type, *background_filter); + Array * disp_per_element = + new Array(0, dim * nodes_per_background_e, "disp_elem"); + FEEngine::extractNodalToElementField( + emodel.getMesh(), emodel.getDisplacement(), *disp_per_element, *back_it, + ghost_type, *background_filter); Array::matrix_iterator disp_it = disp_per_element->begin(dim, nodes_per_background_e); Array::matrix_iterator disp_end = disp_per_element->end(dim, nodes_per_background_e); Array::matrix_iterator shapes_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator grad_u_it = gradu_vec.begin(dim, dim); for (; disp_it != disp_end; ++disp_it) { for (UInt i = 0; i < nb_quad_points; i++, ++shapes_it, ++grad_u_it) { Matrix & B = *shapes_it; Matrix & du = *grad_u_it; Matrix & u = *disp_it; du.mul(u, B); } } delete background_filter; delete disp_per_element; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::computeAllStresses(GhostType ghost_type) { +template +void MaterialReinforcement::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); - Mesh::type_iterator it = model->getInterfaceMesh().firstType(); - Mesh::type_iterator last_type = model->getInterfaceMesh().lastType(); + Mesh::type_iterator it = emodel.getInterfaceMesh().firstType(); + Mesh::type_iterator last_type = emodel.getInterfaceMesh().lastType(); for (; it != last_type; ++it) { computeGradU(*it, ghost_type); - computeStress(*it, ghost_type); + this->computeStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::assembleInternalForces( +template +void MaterialReinforcement::assembleInternalForces( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - Mesh & mesh = model->getMesh(); + Mesh & mesh = emodel.getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); for (; type_it != type_end; ++type_it) { assembleInternalForcesInterface(type, *type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Computes and assemble the residual. Residual in reinforcement is computed as: * * \f[ * \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s} * \f] */ -template -void MaterialReinforcement::assembleInternalForcesInterface( +template +void MaterialReinforcement::assembleInternalForcesInterface( const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - UInt voigt_size = getTangentStiffnessVoigtSize(dim); + UInt voigt_size = VoigtHelper::size; - FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); + FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); - Array & elem_filter = element_filter(interface_type, ghost_type); + Array & elem_filter = this->element_filter(interface_type, ghost_type); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type); UInt nb_element = elem_filter.size(); UInt back_dof = dim * nodes_per_background_e; Array & shapesd = (*shape_derivatives(interface_type, ghost_type))( background_type, ghost_type); Array integrant(nb_quadrature_points * nb_element, back_dof, "integrant"); Array::vector_iterator integrant_it = integrant.begin(back_dof); Array::vector_iterator integrant_end = integrant.end(back_dof); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator C_it = directing_cosines(interface_type, ghost_type) .begin(voigt_size, voigt_size); Array::matrix_iterator sigma_it = stress_embedded(interface_type, ghost_type).begin(dim, dim); Vector sigma(voigt_size); Matrix Bvoigt(voigt_size, back_dof); Vector Ct_sigma(voigt_size); - for (; integrant_it != integrant_end ; ++integrant_it, - ++B_it, - ++C_it, - ++sigma_it) { - VoigtHelper::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt, nodes_per_background_e); + for (; integrant_it != integrant_end; + ++integrant_it, ++B_it, ++C_it, ++sigma_it) { + VoigtHelper::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt, + nodes_per_background_e); Matrix & C = *C_it; Vector & BtCt_sigma = *integrant_it; stressTensorToVoigtVector(*sigma_it, sigma); Ct_sigma.mul(C, sigma); BtCt_sigma.mul(Bvoigt, Ct_sigma); BtCt_sigma *= area; } Array residual_interface(nb_element, back_dof, "residual_interface"); interface_engine.integrate(integrant, residual_interface, back_dof, interface_type, ghost_type, elem_filter); integrant.resize(0); Array background_filter(nb_element, 1, "background_filter"); filterInterfaceBackgroundElements(background_filter, background_type, interface_type, ghost_type); - model->getDOFManager().assembleElementalArrayLocalArray( - residual_interface, model->getInternalForce(), background_type, ghost_type, -1., - background_filter); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::filterInterfaceBackgroundElements( - Array & filter, const ElementType & type, - const ElementType & interface_type, GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - filter.resize(0); - filter.clear(); - - Array & elements = - model->getInterfaceAssociatedElements(interface_type, ghost_type); - Array & elem_filter = element_filter(interface_type, ghost_type); - - Array::scalar_iterator filter_it = elem_filter.begin(), - filter_end = elem_filter.end(); - - for (; filter_it != filter_end; ++filter_it) { - Element & elem = elements(*filter_it); - if (elem.type == type) - filter.push_back(elem.element); - } + emodel.getDOFManager().assembleElementalArrayLocalArray( + residual_interface, emodel.getInternalForce(), background_type, + ghost_type, -1., background_filter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::computeDirectingCosines( +template +void MaterialReinforcement::computeDirectingCosines( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - Mesh & interface_mesh = this->model->getInterfaceMesh(); + Mesh & interface_mesh = emodel.getInterfaceMesh(); const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const UInt steel_dof = dim * nb_nodes_per_element; - const UInt voigt_size = getTangentStiffnessVoigtSize(dim); - const UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine") + const UInt voigt_size = VoigtHelper::size; + const UInt nb_quad_points = emodel.getFEEngine("EmbeddedInterfaceFEEngine") .getNbIntegrationPoints(type, ghost_type); Array node_coordinates(this->element_filter(type, ghost_type).size(), steel_dof); - this->model->getFEEngine().template extractNodalToElementField(interface_mesh, - interface_mesh.getNodes(), - node_coordinates, - type, - ghost_type, - this->element_filter(type, ghost_type)); + this->emodel.getFEEngine().template extractNodalToElementField( + interface_mesh, interface_mesh.getNodes(), node_coordinates, type, + ghost_type, this->element_filter(type, ghost_type)); Array::matrix_iterator directing_cosines_it = directing_cosines(type, ghost_type).begin(voigt_size, voigt_size); Array::matrix_iterator node_coordinates_it = node_coordinates.begin(dim, nb_nodes_per_element); Array::matrix_iterator node_coordinates_end = node_coordinates.end(dim, nb_nodes_per_element); for (; node_coordinates_it != node_coordinates_end; ++node_coordinates_it) { for (UInt i = 0; i < nb_quad_points; i++, ++directing_cosines_it) { Matrix & nodes = *node_coordinates_it; Matrix & cosines = *directing_cosines_it; computeDirectingCosinesOnQuad(nodes, cosines); } } - // Mauro: the directing_cosines internal is defined on the quadrature points of each element + // Mauro: the directing_cosines internal is defined on the quadrature points + // of each element AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::assembleStiffnessMatrix( +template +void MaterialReinforcement::assembleStiffnessMatrix( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - Mesh & mesh = model->getMesh(); + Mesh & mesh = emodel.getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); for (; type_it != type_end; ++type_it) { assembleStiffnessMatrixInterface(type, *type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001) * \f[ * \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T * \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}} * \f] */ -template -void MaterialReinforcement::assembleStiffnessMatrixInterface( +template +void MaterialReinforcement::assembleStiffnessMatrixInterface( const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - UInt voigt_size = getTangentStiffnessVoigtSize(dim); + UInt voigt_size = VoigtHelper::size; - FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); + FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); - Array & elem_filter = element_filter(interface_type, ghost_type); + Array & elem_filter = this->element_filter(interface_type, ghost_type); Array & grad_u = gradu_embedded(interface_type, ghost_type); UInt nb_element = elem_filter.size(); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type); UInt back_dof = dim * nodes_per_background_e; UInt integrant_size = back_dof; grad_u.resize(nb_quadrature_points * nb_element); Array tangent_moduli(nb_element * nb_quadrature_points, 1, "interface_tangent_moduli"); - computeTangentModuli(interface_type, tangent_moduli, ghost_type); + this->computeTangentModuli(interface_type, tangent_moduli, ghost_type); Array & shapesd = (*shape_derivatives(interface_type, ghost_type))( background_type, ghost_type); Array integrant(nb_element * nb_quadrature_points, integrant_size * integrant_size, "B^t*C^t*D*C*B"); /// Temporary matrices for integrant product Matrix Bvoigt(voigt_size, back_dof); Matrix DC(voigt_size, voigt_size); Matrix DCB(voigt_size, back_dof); Matrix CtDCB(voigt_size, back_dof); Array::scalar_iterator D_it = tangent_moduli.begin(); Array::scalar_iterator D_end = tangent_moduli.end(); Array::matrix_iterator C_it = directing_cosines(interface_type, ghost_type) .begin(voigt_size, voigt_size); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator integrant_it = integrant.begin(integrant_size, integrant_size); for (; D_it != D_end; ++D_it, ++C_it, ++B_it, ++integrant_it) { Real & D = *D_it; Matrix & C = *C_it; Matrix & B = *B_it; Matrix & BtCtDCB = *integrant_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix(B, Bvoigt, nodes_per_background_e); DC.clear(); DC(0, 0) = D * area; DC *= C; DCB.mul(DC, Bvoigt); CtDCB.mul(C, DCB); BtCtDCB.mul(Bvoigt, CtDCB); } tangent_moduli.resize(0); - Array K_interface( - nb_element, integrant_size * integrant_size, "K_interface"); + Array K_interface(nb_element, integrant_size * integrant_size, + "K_interface"); interface_engine.integrate(integrant, K_interface, integrant_size * integrant_size, interface_type, ghost_type, elem_filter); integrant.resize(0); // Mauro: Here K_interface contains the local stiffness matrices, // directing_cosines contains the information about the orientation // of the reinforcements, any rotation of the local stiffness matrix // can be done here Array background_filter(nb_element, 1, "background_filter"); filterInterfaceBackgroundElements(background_filter, background_type, interface_type, ghost_type); - model->getDOFManager().assembleElementalMatricesToMatrix( - "K", "displacement", K_interface, background_type, ghost_type, - _symmetric, background_filter); - + emodel.getDOFManager().assembleElementalMatricesToMatrix( + "K", "displacement", K_interface, background_type, ghost_type, _symmetric, + background_filter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ - -/// In this function, type and ghost_type refer to background elements -template -void MaterialReinforcement::computeBackgroundShapeDerivatives( - const ElementType & type, GhostType ghost_type) { +template +Real MaterialReinforcement::getEnergy(const std::string & id) { AKANTU_DEBUG_IN(); + if (id == "potential") { + Real epot = 0.; - Mesh & interface_mesh = model->getInterfaceMesh(); - - FEEngine & engine = model->getFEEngine(); - FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); + this->computePotentialEnergyByElements(); - Mesh::type_iterator interface_type = interface_mesh.firstType(); - Mesh::type_iterator interface_last = interface_mesh.lastType(); + Mesh::type_iterator it = this->element_filter.firstType( + this->spatial_dimension), + end = this->element_filter.lastType( + this->spatial_dimension); - for (; interface_type != interface_last; ++interface_type) { - Array & filter = element_filter(*interface_type, ghost_type); + for (; it != end; ++it) { + FEEngine & interface_engine = + emodel.getFEEngine("EmbeddedInterfaceFEEngine"); + epot += interface_engine.integrate( + this->potential_energy(*it, _not_ghost), *it, _not_ghost, + this->element_filter(*it, _not_ghost)); + epot *= area; + } - const UInt nb_elements = filter.size(); - const UInt nb_nodes = Mesh::getNbNodesPerElement(type); - const UInt nb_quad_per_element = - interface_engine.getNbIntegrationPoints(*interface_type); + return epot; + } - Array quad_pos(nb_quad_per_element * nb_elements, dim, - "interface_quad_points"); - quad_pos.resize(nb_quad_per_element * nb_elements); - interface_engine.interpolateOnIntegrationPoints( - interface_mesh.getNodes(), quad_pos, dim, *interface_type, ghost_type, - filter); + AKANTU_DEBUG_OUT(); + return 0; +} - Array & background_shapesd = - shape_derivatives(*interface_type, ghost_type) - -> - operator()(type, ghost_type); - background_shapesd.clear(); +/* -------------------------------------------------------------------------- */ - Array * background_elements = new Array( - nb_elements, 1, "computeBackgroundShapeDerivatives:background_filter"); - filterInterfaceBackgroundElements(*background_elements, type, - *interface_type, ghost_type); +/** + * The structure of the directing cosines matrix is : + * \f{eqnarray*}{ + * C_{1,\cdot} & = & (l^2, m^2, n^2, mn, ln, lm) \\ + * C_{i,j} & = & 0 + * \f} + * + * with : + * \f[ + * (l, m, n) = \frac{1}{\|\frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}\|} \cdot + * \frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s} + * \f] + */ +template +inline void MaterialReinforcement::computeDirectingCosinesOnQuad( + const Matrix & nodes, Matrix & cosines) { + AKANTU_DEBUG_IN(); - Array::scalar_iterator back_it = background_elements->begin(), - back_end = background_elements->end(); + AKANTU_DEBUG_ASSERT(nodes.cols() == 2, + "Higher order reinforcement elements not implemented"); - Array::matrix_iterator shapesd_it = - background_shapesd.begin(dim, nb_nodes); + const Vector a = nodes(0), b = nodes(1); + Vector delta = b - a; - Array::vector_iterator quad_pos_it = quad_pos.begin(dim); + cosines.clear(); - for (; back_it != back_end; ++back_it) { - for (UInt i = 0; i < nb_quad_per_element; - i++, ++shapesd_it, ++quad_pos_it) - engine.computeShapeDerivatives(*quad_pos_it, *back_it, type, - *shapesd_it, ghost_type); - } + Real sq_length = 0.; + for (UInt i = 0; i < dim; i++) { + sq_length += delta(i) * delta(i); + } - delete background_elements; + if (dim == 2) { + cosines(0, 0) = delta(0) * delta(0); // l^2 + cosines(0, 1) = delta(1) * delta(1); // m^2 + cosines(0, 2) = delta(0) * delta(1); // lm + } else if (dim == 3) { + cosines(0, 0) = delta(0) * delta(0); // l^2 + cosines(0, 1) = delta(1) * delta(1); // m^2 + cosines(0, 2) = delta(2) * delta(2); // n^2 + + cosines(0, 3) = delta(1) * delta(2); // mn + cosines(0, 4) = delta(0) * delta(2); // ln + cosines(0, 5) = delta(0) * delta(1); // lm } + cosines /= sq_length; + AKANTU_DEBUG_OUT(); } -template -Real MaterialReinforcement::getEnergy(const std::string & id) { - AKANTU_DEBUG_IN(); - if (id == "potential") { - Real epot = 0.; - - computePotentialEnergyByElements(); +/* -------------------------------------------------------------------------- */ - Mesh::type_iterator it = element_filter.firstType(spatial_dimension), - end = element_filter.lastType(spatial_dimension); +template +inline void MaterialReinforcement::stressTensorToVoigtVector( + const Matrix & tensor, Vector & vector) { + AKANTU_DEBUG_IN(); - for (; it != end; ++it) { - FEEngine & interface_engine = - model->getFEEngine("EmbeddedInterfaceFEEngine"); - epot += interface_engine.integrate(potential_energy(*it, _not_ghost), *it, - _not_ghost, - element_filter(*it, _not_ghost)); - epot *= area; - } + for (UInt i = 0; i < dim; i++) { + vector(i) = tensor(i, i); + } - return epot; + if (dim == 2) { + vector(2) = tensor(0, 1); + } else if (dim == 3) { + vector(3) = tensor(1, 2); + vector(4) = tensor(0, 2); + vector(5) = tensor(0, 1); } AKANTU_DEBUG_OUT(); - return 0; } -// /* -------------------------------------------------------------------------- -// */ -// template -// ElementTypeMap MaterialReinforcement::getInternalDataPerElem(const -// ID & field_name, -// const -// ElementKind -// & -// kind, -// const -// ID & -// fe_engine_id) -// const -// { -// return Material::getInternalDataPerElem(field_name, kind, -// "EmbeddedInterfaceFEEngine"); -// } - -// /* -------------------------------------------------------------------------- -// */ -// // Author is Guillaume Anciaux, see material.cc -// template -// void MaterialReinforcement::flattenInternal(const std::string & -// field_id, -// ElementTypeMapArray & -// internal_flat, -// const GhostType ghost_type, -// ElementKind element_kind) -// const { -// AKANTU_DEBUG_IN(); - -// if (field_id == "stress_embedded" || field_id == "inelastic_strain") { -// Material::flattenInternalIntern(field_id, -// internal_flat, -// 1, -// ghost_type, -// _ek_not_defined, -// &(this->element_filter), -// &(this->model->getInterfaceMesh())); -// } - -// AKANTU_DEBUG_OUT(); -// } - /* -------------------------------------------------------------------------- */ -INSTANTIATE_MATERIAL_ONLY(MaterialReinforcement); +template +inline void MaterialReinforcement::strainTensorToVoigtVector( + const Matrix & tensor, Vector & vector) { + AKANTU_DEBUG_IN(); -/* -------------------------------------------------------------------------- */ + for (UInt i = 0; i < dim; i++) { + vector(i) = tensor(i, i); + } + + if (dim == 2) { + vector(2) = 2 * tensor(0, 1); + } else if (dim == 3) { + vector(3) = 2 * tensor(1, 2); + vector(4) = 2 * tensor(0, 2); + vector(5) = 2 * tensor(0, 1); + } + + AKANTU_DEBUG_OUT(); +} } // akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc index f3418e008..d774a127a 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc @@ -1,174 +1,175 @@ /** * @file embedded_interface_model.cc * * @author Lucas Frerot * * @date creation: Fri Mar 13 2015 * @date last modification: Mon Dec 14 2015 * * @brief Model of Solid Mechanics with embedded interfaces * * @section LICENSE * * Copyright (©) 2015 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 "embedded_interface_model.hh" -#include "material_reinforcement_template.hh" +#include "material_reinforcement.hh" +#include "material_elastic.hh" #include "mesh_iterators.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER # include "dumper_iohelper_paraview.hh" # include "dumpable_inline_impl.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ EmbeddedInterfaceModel::EmbeddedInterfaceModel(Mesh & mesh, Mesh & primitive_mesh, UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, spatial_dimension, id, memory_id), intersector(mesh, primitive_mesh), interface_mesh(nullptr), primitive_mesh(primitive_mesh), interface_material_selector(nullptr) { this->model_type = ModelType::_embedded_model; // This pointer should be deleted by ~SolidMechanicsModel() auto mat_sel_pointer = std::make_shared>("physical_names", *this); this->setMaterialSelector(mat_sel_pointer); interface_mesh = &(intersector.getInterfaceMesh()); // Create 1D FEEngine on the interface mesh registerFEEngineObject("EmbeddedInterfaceFEEngine", *interface_mesh, 1); // Registering allocator for material reinforcement MaterialFactory::getInstance().registerAllocator( "reinforcement", [&](UInt dim, const ID & constitutive, SolidMechanicsModel &, const ID & id) -> std::unique_ptr { if (constitutive == "elastic") { using mat = MaterialElastic<1>; switch (dim) { case 2: - return std::unique_ptr>{ - new MaterialReinforcementTemplate<2, mat>(*this, id)}; + return std::unique_ptr>{ + new MaterialReinforcement(*this, id)}; case 3: - return std::unique_ptr>{ - new MaterialReinforcementTemplate<3, mat>(*this, id)}; + return std::unique_ptr>{ + new MaterialReinforcement(*this, id)}; default: AKANTU_EXCEPTION("Dimension 1 is invalid for reinforcements"); } } else { AKANTU_EXCEPTION("Reinforcement type" << constitutive << " is not recognized"); } }); } /* -------------------------------------------------------------------------- */ EmbeddedInterfaceModel::~EmbeddedInterfaceModel() { delete interface_material_selector; } /* -------------------------------------------------------------------------- */ void EmbeddedInterfaceModel::initFullImpl(const ModelOptions & options) { const auto & eim_options = dynamic_cast(options); // Do no initialize interface_mesh if told so if (eim_options.has_intersections) intersector.constructData(); SolidMechanicsModel::initFullImpl(options); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("reinforcement", id); this->mesh.addDumpMeshToDumper("reinforcement", *interface_mesh, 1, _not_ghost, _ek_regular); #endif } void EmbeddedInterfaceModel::initModel() { // Initialize interface FEEngine FEEngine & engine = getFEEngine("EmbeddedInterfaceFEEngine"); engine.initShapeFunctions(_not_ghost); engine.initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void EmbeddedInterfaceModel::assignMaterialToElements( const ElementTypeMapArray * filter) { delete interface_material_selector; interface_material_selector = new InterfaceMeshDataMaterialSelector("physical_names", *this); for_each_element(mesh, [&](auto && element) { auto mat_index = (*interface_material_selector)(element); material_index(element) = mat_index; materials[mat_index]->addElement(element); }, _element_filter = filter); SolidMechanicsModel::assignMaterialToElements(filter); } /* -------------------------------------------------------------------------- */ void EmbeddedInterfaceModel::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { #ifdef AKANTU_USE_IOHELPER dumper::Field * field = NULL; // If dumper is reinforcement, create a 1D elemental field if (dumper_name == "reinforcement") field = this->createElementalField(field_id, group_name, padding_flag, 1, element_kind); else { try { SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, element_kind, padding_flag); } catch (...) {} } if (field) { DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name,group_name); Model::addDumpGroupFieldToDumper(field_id,field,dumper); } #endif } } // akantu