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..470987bd4 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,437 +1,440 @@ /** * @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, + template + inline Stored & operator()(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..cedbb7f98 100644 --- a/src/mesh/element_type_map_tmpl.hh +++ b/src/mesh/element_type_map_tmpl.hh @@ -1,661 +1,662 @@ /** * @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, - const GhostType & ghost_type) { +template +inline Stored & ElementTypeMap::operator()( + 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/material.cc b/src/model/solid_mechanics/material.cc index ac464fc04..51cbada6c 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1553 +1,1553 @@ /** * @file material.cc * * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Nov 24 2015 * * @brief Implementation of the common part of the material class * * @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 "material.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_material, id), is_init(false), fem(model.getFEEngine()), finite_deformation(false), name(""), model(model), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id, this->memory_id), stress("stress", *this), eigengradu("eigen_grad_u", *this), gradu("grad_u", *this), green_strain("green_strain", *this), piola_kirchhoff_2("piola_kirchhoff_2", *this), potential_energy("potential_energy", *this), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse coordinates", *this), interpolation_points_matrices("interpolation points matrices", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(model.getMesh(), _spatial_dimension = spatial_dimension); // model.getMesh().initElementTypeMapArray(element_filter, 1, // spatial_dimension, // false, _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_material, id), - is_init(false), fem(model.getFEEngine()), finite_deformation(false), + is_init(false), fem(fe_engine), finite_deformation(false), name(""), model(model), spatial_dimension(dim), element_filter("element_filter", id, this->memory_id), stress("stress", *this, dim, fe_engine, this->element_filter), eigengradu("eigen_grad_u", *this, dim, fe_engine, this->element_filter), gradu("gradu", *this, dim, fe_engine, this->element_filter), green_strain("green_strain", *this, dim, fe_engine, this->element_filter), - piola_kirchhoff_2("poila_kirchhoff_2", *this, dim, fe_engine, + piola_kirchhoff_2("piola_kirchhoff_2", *this, dim, fe_engine, this->element_filter), potential_energy("potential_energy", *this, dim, fe_engine, this->element_filter), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse_coordinates", *this, dim, fe_engine, this->element_filter), interpolation_points_matrices("interpolation points matrices", *this, dim, fe_engine, this->element_filter) { AKANTU_DEBUG_IN(); element_filter.initialize(mesh, _spatial_dimension = spatial_dimension); // mesh.initElementTypeMapArray(element_filter, 1, spatial_dimension, false, // _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() = default; /* -------------------------------------------------------------------------- */ void Material::initialize() { registerParam("rho", rho, Real(0.), _pat_parsable | _pat_modifiable, "Density"); registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("finite_deformation", finite_deformation, false, _pat_parsable | _pat_readable, "Is finite deformation"); registerParam("inelastic_deformation", inelastic_deformation, false, _pat_internal, "Is inelastic deformation"); /// allocate gradu stress for local elements eigengradu.initialize(spatial_dimension * spatial_dimension); gradu.initialize(spatial_dimension * spatial_dimension); stress.initialize(spatial_dimension * spatial_dimension); potential_energy.initialize(1); this->model.registerEventHandler(*this); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); if (finite_deformation) { this->piola_kirchhoff_2.initialize(spatial_dimension * spatial_dimension); if (use_previous_stress) this->piola_kirchhoff_2.initializeHistory(); this->green_strain.initialize(spatial_dimension * spatial_dimension); } if (use_previous_stress) this->stress.initializeHistory(); if (use_previous_gradu) this->gradu.initializeHistory(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::savePreviousState() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { if (it->second->hasHistory()) it->second->saveCurrentValues(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ // void Material::updateResidual(GhostType ghost_type) { // AKANTU_DEBUG_IN(); // computeAllStresses(ghost_type); // assembleResidual(ghost_type); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); if (!finite_deformation) { auto & internal_force = const_cast &>(model.getInternalForce()); // Mesh & mesh = fem.getMesh(); for (auto && type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) continue; const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by /// @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Array * sigma_dphi_dx = new Array(nb_element * nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); fem.computeBtD(stress(type, ghost_type), *sigma_dphi_dx, type, ghost_type, elem_filter); // Array * shapesd_filtered = // new Array(0, size_of_shapes_derivatives, "filtered shapesd"); // FEEngine::filterElementalData(mesh, shapes_derivatives, // *shapesd_filtered, // *it, ghost_type, elem_filter); // Array & stress_vect = this->stress(*it, ghost_type); // Array::matrix_iterator sigma = // stress_vect.begin(spatial_dimension, spatial_dimension); // Array::matrix_iterator B = // shapesd_filtered->begin(spatial_dimension, nb_nodes_per_element); // Array::matrix_iterator Bt_sigma_it = // sigma_dphi_dx->begin(spatial_dimension, nb_nodes_per_element); // for (UInt q = 0; q < nb_element * nb_quadrature_points; // ++q, ++sigma, ++B, ++Bt_sigma_it) // Bt_sigma_it->mul(*sigma, *B); // delete shapesd_filtered; /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by * @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ Array * int_sigma_dphi_dx = new Array(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); fem.integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, type, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble model.getDOFManager().assembleElementalArrayLocalArray( *int_sigma_dphi_dx, internal_force, type, ghost_type, -1, elem_filter); delete int_sigma_dphi_dx; } } else { switch (spatial_dimension) { case 1: this->assembleInternalForces<1>(ghost_type); break; case 2: this->assembleInternalForces<2>(ghost_type); break; case 3: this->assembleInternalForces<3>(ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the gradu * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (const auto & type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, spatial_dimension, type, ghost_type, elem_filter); gradu_vect -= eigengradu(type, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllCauchyStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(finite_deformation, "The Cauchy stress can only be " "computed if you are working in " "finite deformation."); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { switch (spatial_dimension) { case 1: this->computeCauchyStress<1>(type, ghost_type); break; case 2: this->computeCauchyStress<2>(type, ghost_type); break; case 3: this->computeCauchyStress<3>(type, ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeCauchyStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::matrix_iterator gradu_it = this->gradu(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator gradu_end = this->gradu(el_type, ghost_type).end(dim, dim); Array::matrix_iterator piola_it = this->piola_kirchhoff_2(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator stress_it = this->stress(el_type, ghost_type).begin(dim, dim); Matrix F_tensor(dim, dim); for (; gradu_it != gradu_end; ++gradu_it, ++piola_it, ++stress_it) { Matrix & grad_u = *gradu_it; Matrix & piola = *piola_it; Matrix & sigma = *stress_it; gradUToF(grad_u, F_tensor); this->computeCauchyStressOnQuad(F_tensor, piola, sigma); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & displacement = model.getDisplacement(); // resizeInternalArray(gradu); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(displacement, gradu_vect, spatial_dimension, type, ghost_type, elem_filter); setToSteadyState(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { if (finite_deformation) { switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL<1>(type, ghost_type); assembleStiffnessMatrixL2<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrixNL<2>(type, ghost_type); assembleStiffnessMatrixL2<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrixNL<3>(type, ghost_type); assembleStiffnessMatrixL2<3>(type, ghost_type); break; } } } else { switch (spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(type, ghost_type); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) { AKANTU_DEBUG_OUT(); return; } // const Array & shapes_derivatives = // fem.getShapesDerivatives(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); // Array * shapesd_filtered = new Array( // nb_element, dim * nb_nodes_per_element, "filtered shapesd"); // fem.filterElementalData(fem.getMesh(), shapes_derivatives, // *shapesd_filtered, type, ghost_type, // elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); fem.computeBtDB(*tangent_stiffness_matrix, *bt_d_b, 4, type, ghost_type, elem_filter); // Matrix B(tangent_size, dim * nb_nodes_per_element); // Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); // Array::matrix_iterator shapes_derivatives_filtered_it = // shapesd_filtered->begin(dim, nb_nodes_per_element); // Array::matrix_iterator Bt_D_B_it = // bt_d_b->begin(dim * nb_nodes_per_element, dim * nb_nodes_per_element); // Array::matrix_iterator D_it = // tangent_stiffness_matrix->begin(tangent_size, tangent_size); // Array::matrix_iterator D_end = // tangent_stiffness_matrix->end(tangent_size, tangent_size); // for (; D_it != D_end; ++D_it, ++Bt_D_B_it, // ++shapes_derivatives_filtered_it) { // Matrix & D = *D_it; // Matrix & Bt_D_B = *Bt_D_B_it; // VoigtHelper::transferBMatrixToSymVoigtBMatrix( // *shapes_derivatives_filtered_it, B, nb_nodes_per_element); // Bt_D.mul(B, D); // Bt_D_B.mul(Bt_D, B); // } delete tangent_stiffness_matrix; // delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); // Array & gradu_vect = delta_gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); fem.filterElementalData(fem.getMesh(), shapes_derivatives, *shapes_derivatives_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_s_b_size = dim * nb_nodes_per_element; Array * bt_s_b = new Array(nb_element * nb_quadrature_points, bt_s_b_size * bt_s_b_size, "B^t*D*B"); UInt piola_matrix_size = getCauchyStressMatrixSize(dim); Matrix B(piola_matrix_size, bt_s_b_size); Matrix Bt_S(bt_s_b_size, piola_matrix_size); Matrix S(piola_matrix_size, piola_matrix_size); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); auto Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); auto piola_it = piola_kirchhoff_2(type, ghost_type).begin(dim, dim); for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) { auto & Bt_S_B = *Bt_S_B_it; const auto & Piola_kirchhoff_matrix = *piola_it; setCauchyStressMatrix(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.template mul(B, S); Bt_S_B.template mul(Bt_S, B); } delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_s_b_size * bt_s_b_size, "K_e"); fem.integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, elem_filter); delete bt_s_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); fem.filterElementalData(fem.getMesh(), shapes_derivatives, *shapes_derivatives_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_D_B_it = bt_d_b->begin(bt_d_b_size, bt_d_b_size); auto grad_u_it = gradu_vect.begin(dim, dim); auto D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); auto D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) { const auto & grad_u = *grad_u_it; const auto & D = *D_it; auto & Bt_D_B = *Bt_D_B_it; // transferBMatrixToBL1 (*shapes_derivatives_filtered_it, B, // nb_nodes_per_element); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element); B += B2; Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); Mesh & mesh = fem.getMesh(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array * shapesd_filtered = new Array( nb_element, size_of_shapes_derivatives, "filtered shapesd"); fem.filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); // Set stress vectors UInt stress_size = getTangentStiffnessVoigtSize(dim); // Set matrices B and BNL* UInt bt_s_size = dim * nb_nodes_per_element; auto * bt_s = new Array(nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); auto grad_u_it = this->gradu(type, ghost_type).begin(dim, dim); auto grad_u_end = this->gradu(type, ghost_type).end(dim, dim); auto stress_it = this->piola_kirchhoff_2(type, ghost_type).begin(dim, dim); shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1); Matrix S_vect(stress_size, 1); Matrix B_tensor(stress_size, bt_s_size); Matrix B2_tensor(stress_size, bt_s_size); for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) { auto & grad_u = *grad_u_it; auto & r_it = *bt_s_it; auto & S_it = *stress_it; setCauchyStressArray(S_it, S_vect); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element); B_tensor += B2_tensor; r_it.template mul(B_tensor, S_vect); } delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * r_e = new Array(nb_element, bt_s_size, "r_e"); fem.integrate(*bt_s, *r_e, bt_s_size, type, ghost_type, elem_filter); delete bt_s; model.getDOFManager().assembleElementalArrayLocalArray( *r_e, internal_force, type, ghost_type, -1., elem_filter); delete r_e; } AKANTU_DEBUG_OUT(); } // /* -------------------------------------------------------------------------- */ // void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) { // AKANTU_DEBUG_IN(); // UInt spatial_dimension = model.getSpatialDimension(); // for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { // switch (spatial_dimension) { // case 1: { // computeAllStressesFromTangentModuli<1>(type, ghost_type); // break; // } // case 2: { // computeAllStressesFromTangentModuli<2>(type, ghost_type); // break; // } // case 3: { // computeAllStressesFromTangentModuli<3>(type, ghost_type); // break; // } // } // } // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- */ // template // void Material::computeAllStressesFromTangentModuli(const ElementType & type, // GhostType ghost_type) { // AKANTU_DEBUG_IN(); // const Array & shapes_derivatives = // fem.getShapesDerivatives(type, ghost_type); // Array & elem_filter = element_filter(type, ghost_type); // Array & gradu_vect = gradu(type, ghost_type); // UInt nb_element = elem_filter.size(); // if (nb_element) { // UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); // UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); // gradu_vect.resize(nb_quadrature_points * nb_element); // Array & disp = model.getDisplacement(); // fem.gradientOnIntegrationPoints(disp, gradu_vect, dim, type, ghost_type, // elem_filter); // UInt tangent_moduli_size = getTangentStiffnessVoigtSize(dim); // Array * tangent_moduli_tensors = new Array( // nb_element * nb_quadrature_points, // tangent_moduli_size * tangent_moduli_size, "tangent_moduli_tensors"); // tangent_moduli_tensors->clear(); // computeTangentModuli(type, *tangent_moduli_tensors, ghost_type); // Array * shapesd_filtered = new Array( // nb_element, dim * nb_nodes_per_element, "filtered shapesd"); // FEEngine::filterElementalData(fem.getMesh(), shapes_derivatives, // *shapesd_filtered, type, ghost_type, // elem_filter); // Array filtered_u(nb_element, // nb_nodes_per_element * spatial_dimension); // fem.extractNodalToElementField(fem.getMesh(), disp, filtered_u, type, // ghost_type, elem_filter); // /// compute @f$\mathbf{D} \mathbf{B} \mathbf{u}@f$ // auto shapes_derivatives_filtered_it = // shapesd_filtered->begin(dim, nb_nodes_per_element); // auto D_it = // tangent_moduli_tensors->begin(tangent_moduli_size, tangent_moduli_size); // auto sigma_it = // stress(type, ghost_type).begin(spatial_dimension, spatial_dimension); // auto u_it = filtered_u.begin(spatial_dimension * nb_nodes_per_element); // Matrix B(tangent_moduli_size, // spatial_dimension * nb_nodes_per_element); // Vector Bu(tangent_moduli_size); // Vector DBu(tangent_moduli_size); // for (UInt e = 0; e < nb_element; ++e, ++u_it) { // for (UInt q = 0; q < nb_quadrature_points; // ++q, ++D_it, ++shapes_derivatives_filtered_it, ++sigma_it) { // const auto & u = *u_it; // const auto & D = *D_it; // auto & sigma = *sigma_it; // VoigtHelper::transferBMatrixToSymVoigtBMatrix( // *shapes_derivatives_filtered_it, B, nb_nodes_per_element); // Bu.mul(B, u); // DBu.mul(D, Bu); // // Voigt notation to full symmetric tensor // for (UInt i = 0; i < dim; ++i) // sigma(i, i) = DBu(i); // if (dim == 2) { // sigma(0, 1) = sigma(1, 0) = DBu(2); // } else if (dim == 3) { // sigma(1, 2) = sigma(2, 1) = DBu(3); // sigma(0, 2) = sigma(2, 0) = DBu(4); // sigma(0, 1) = sigma(1, 0) = DBu(5); // } // } // } // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElements() { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { computePotentialEnergy(type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType, GhostType) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElements(); /// integrate the potential energy for each type of elements for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { epot += fem.integrate(potential_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy(ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real epot = 0.; Vector epot_on_quad_points(fem.getNbIntegrationPoints(type)); computePotentialEnergyByElement(type, index, epot_on_quad_points); epot = fem.integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(const std::string & type) { AKANTU_DEBUG_IN(); if (type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(const std::string & energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "potential") return getPotentialEnergy(type, index); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation( const ElementTypeMapArray & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); this->fem.initElementalFieldInterpolationFromIntegrationPoints( interpolation_points_coordinates, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, &(this->element_filter)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type) { this->fem.interpolateElementalFieldFromIntegrationPoints( this->stress, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, result, ghost_type, &(this->element_filter)); } /* -------------------------------------------------------------------------- */ void Material::interpolateStressOnFacets( ElementTypeMapArray & result, ElementTypeMapArray & by_elem_result, const GhostType ghost_type) { interpolateStress(by_elem_result, ghost_type); UInt stress_size = this->stress.getNbComponent(); const Mesh & mesh = this->model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_fil = element_filter(type, ghost_type); Array & by_elem_res = by_elem_result(type, ghost_type); UInt nb_element = elem_fil.size(); UInt nb_element_full = this->model.getMesh().getNbElement(type, ghost_type); UInt nb_interpolation_points_per_elem = by_elem_res.size() / nb_element_full; const Array & facet_to_element = mesh_facets.getSubelementToElement(type, ghost_type); ElementType type_facet = Mesh::getFacetType(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); UInt nb_quad_per_facet = nb_interpolation_points_per_elem / nb_facet_per_elem; Element element_for_comparison{type, 0, ghost_type}; const Array> * element_to_facet = nullptr; GhostType current_ghost_type = _casper; Array * result_vec = nullptr; Array::const_matrix_iterator result_it = by_elem_res.begin_reinterpret( stress_size, nb_interpolation_points_per_elem, nb_element_full); for (UInt el = 0; el < nb_element; ++el) { UInt global_el = elem_fil(el); element_for_comparison.element = global_el; for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element facet_elem = facet_to_element(global_el, f); UInt global_facet = facet_elem.element; if (facet_elem.ghost_type != current_ghost_type) { current_ghost_type = facet_elem.ghost_type; element_to_facet = &mesh_facets.getElementToSubelement( type_facet, current_ghost_type); result_vec = &result(type_facet, current_ghost_type); } bool is_second_element = (*element_to_facet)(global_facet)[0] != element_for_comparison; for (UInt q = 0; q < nb_quad_per_facet; ++q) { Vector result_local(result_vec->storage() + (global_facet * nb_quad_per_facet + q) * result_vec->getNbComponent() + is_second_element * stress_size, stress_size); const Matrix & result_tmp(result_it[global_el]); result_local = result_tmp(f * nb_quad_per_facet + q); } } } } } /* -------------------------------------------------------------------------- */ template const Array & Material::getArray(__attribute__((unused)) const ID & vect_id, __attribute__((unused)) const ElementType & type, __attribute__((unused)) const GhostType & ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template Array & Material::getArray(__attribute__((unused)) const ID & vect_id, __attribute__((unused)) const ElementType & type, __attribute__((unused)) const GhostType & ghost_type) { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template const InternalField & Material::getInternal(__attribute__((unused)) const ID & int_id) const { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template InternalField & Material::getInternal(__attribute__((unused)) const ID & int_id) { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ void Material::addElements(const Array & elements_to_add) { AKANTU_DEBUG_IN(); UInt mat_id = model.getInternalIndexFromID(getID()); Array::const_iterator el_begin = elements_to_add.begin(); Array::const_iterator el_end = elements_to_add.end(); for (; el_begin != el_end; ++el_begin) { const Element & element = *el_begin; Array & mat_indexes = model.getMaterialByElement(element.type, element.ghost_type); Array & mat_loc_num = model.getMaterialLocalNumbering(element.type, element.ghost_type); UInt index = this->addElement(element.type, element.element, element.ghost_type); mat_indexes(element.element) = mat_id; mat_loc_num(element.element) = index; } this->resizeInternals(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::removeElements(const Array & elements_to_remove) { AKANTU_DEBUG_IN(); Array::const_iterator el_begin = elements_to_remove.begin(); Array::const_iterator el_end = elements_to_remove.end(); if (el_begin == el_end) return; ElementTypeMapArray material_local_new_numbering( "remove mat filter elem", getID(), getMemoryID()); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; element.ghost_type = ghost_type; ElementTypeMapArray::type_iterator it = element_filter.firstType(_all_dimensions, ghost_type, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, ghost_type, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; element.type = type; Array & elem_filter = this->element_filter(type, ghost_type); Array & mat_loc_num = this->model.getMaterialLocalNumbering(type, ghost_type); if (!material_local_new_numbering.exists(type, ghost_type)) material_local_new_numbering.alloc(elem_filter.size(), 1, type, ghost_type); Array & mat_renumbering = material_local_new_numbering(type, ghost_type); UInt nb_element = elem_filter.size(); Array elem_filter_tmp; UInt new_id = 0; for (UInt el = 0; el < nb_element; ++el) { element.element = elem_filter(el); if (std::find(el_begin, el_end, element) == el_end) { elem_filter_tmp.push_back(element.element); mat_renumbering(el) = new_id; mat_loc_num(element.element) = new_id; ++new_id; } else { mat_renumbering(el) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.size()); elem_filter.copy(elem_filter_tmp); } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resizeInternals() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::onElementsAdded(const Array &, const NewElementsEvent &) { this->resizeInternals(); } /* -------------------------------------------------------------------------- */ void Material::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { UInt my_num = model.getInternalIndexFromID(getID()); ElementTypeMapArray material_local_new_numbering( "remove mat filter elem", getID(), getMemoryID()); Array::const_iterator el_begin = element_list.begin(); Array::const_iterator el_end = element_list.end(); for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType gt = *g; ElementTypeMapArray::type_iterator it = new_numbering.firstType(_all_dimensions, gt, _ek_not_defined); ElementTypeMapArray::type_iterator end = new_numbering.lastType(_all_dimensions, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if (element_filter.exists(type, gt) && element_filter(type, gt).size()) { Array & elem_filter = element_filter(type, gt); Array & mat_indexes = this->model.getMaterialByElement(*it, gt); Array & mat_loc_num = this->model.getMaterialLocalNumbering(*it, gt); UInt nb_element = this->model.getMesh().getNbElement(type, gt); // all materials will resize of the same size... mat_indexes.resize(nb_element); mat_loc_num.resize(nb_element); if (!material_local_new_numbering.exists(type, gt)) material_local_new_numbering.alloc(elem_filter.size(), 1, type, gt); Array & mat_renumbering = material_local_new_numbering(type, gt); const Array & renumbering = new_numbering(type, gt); Array elem_filter_tmp; UInt ni = 0; Element el{type, 0, gt}; for (UInt i = 0; i < elem_filter.size(); ++i) { el.element = elem_filter(i); if (std::find(el_begin, el_end, el) == el_end) { UInt new_el = renumbering(el.element); AKANTU_DEBUG_ASSERT( new_el != UInt(-1), "A not removed element as been badly renumbered"); elem_filter_tmp.push_back(new_el); mat_renumbering(i) = ni; mat_indexes(new_el) = my_num; mat_loc_num(new_el) = ni; ++ni; } else { mat_renumbering(i) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.size()); elem_filter.copy(elem_filter_tmp); } } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); } /* -------------------------------------------------------------------------- */ void Material::beforeSolveStep() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::afterSolveStep() { for (auto & type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { this->updateEnergies(type, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDamageIteration() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onDamageUpdate() { ElementTypeMapArray::type_iterator it = this->element_filter.firstType( _all_dimensions, _not_ghost, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined); for (; it != end; ++it) { this->updateEnergiesAfterDamage(*it, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDump() { if (this->isFiniteDeformation()) this->computeAllCauchyStresses(_not_ghost); } /* -------------------------------------------------------------------------- */ void Material::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /// extrapolate internal values void Material::extrapolateInternal(const ID & id, const Element & element, __attribute__((unused)) const Matrix & point, Matrix & extrapolated) { if (this->isInternal(id, element.kind())) { UInt nb_element = this->element_filter(element.type, element.ghost_type).size(); const ID name = this->getID() + ":" + id; UInt nb_quads = this->internal_vectors_real[name]->getFEEngine().getNbIntegrationPoints( element.type, element.ghost_type); const Array & internal = this->getArray(id, element.type, element.ghost_type); UInt nb_component = internal.getNbComponent(); Array::const_matrix_iterator internal_it = internal.begin_reinterpret(nb_component, nb_quads, nb_element); Element local_element = this->convertToLocalElement(element); /// instead of really extrapolating, here the value of the first GP /// is copied into the result vector. This works only for linear /// elements /// @todo extrapolate!!!! AKANTU_DEBUG_WARNING("This is a fix, values are not truly extrapolated"); const Matrix & values = internal_it[local_element.element]; UInt index = 0; Vector tmp(nb_component); for (UInt j = 0; j < values.cols(); ++j) { tmp = values(j); if (tmp.norm() > 0) { index = j; break; } } for (UInt i = 0; i < extrapolated.size(); ++i) { extrapolated(i) = values(index); } } else { Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); extrapolated = default_values; } } /* -------------------------------------------------------------------------- */ void Material::applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const GhostType ghost_type) { for (auto && type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { if (!element_filter(type, ghost_type).size()) continue; auto eigen_it = this->eigengradu(type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto eigen_end = this->eigengradu(type, ghost_type) .end(spatial_dimension, spatial_dimension); for (; eigen_it != eigen_end; ++eigen_it) { auto & current_eigengradu = *eigen_it; current_eigengradu = prescribed_eigen_grad_u; } } } } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc deleted file mode 100644 index def0887e8..000000000 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc +++ /dev/null @@ -1,766 +0,0 @@ -/** - * @file material_reinforcement.cc - * - * @author Lucas Frerot - * - * @date creation: Wed Mar 25 2015 - * @date last modification: Tue Dec 08 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 . - * - */ - -/* -------------------------------------------------------------------------- */ - -#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, - this->element_filter), - gradu_embedded("gradu_embedded", *this, 1, fe_engine, - this->element_filter), - directing_cosines("directing_cosines", *this, 1, fe_engine, - this->element_filter), - pre_stress("pre_stress", *this, 1, fe_engine, this->element_filter), - area(1.0), shape_derivatives() { - AKANTU_DEBUG_IN(); - this->initialize(model); - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::initialize(SolidMechanicsModel & a_model) { - 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); - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ - -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() { - 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->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(); - - Mesh & interface_mesh = model->getInterfaceMesh(); - Mesh & mesh = model->getMesh(); - - ghost_type_t::iterator int_ghost_it = ghost_type_t::begin(); - - // 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); - - // 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); - - // 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; - - std::string shaped_id = "embedded_shape_derivatives"; - - if (int_ghost == _ghost) shaped_id += ":ghost"; - - ElementTypeMapArray * shaped_etma = new ElementTypeMapArray(shaped_id, this->name); - - 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); - } - } - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ - -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); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ - -template void MaterialReinforcement::initDirectingCosines() { - AKANTU_DEBUG_IN(); - - Mesh & mesh = model->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); - 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) { - AKANTU_DEBUG_IN(); - - Mesh & interface_mesh = model->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) { - AKANTU_DEBUG_IN(); - - Mesh & interface_mesh = model->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) { - AKANTU_DEBUG_IN(); - - Array & elem_filter = element_filter(type, ghost_type); - UInt nb_element = elem_filter.size(); - UInt nb_quad_points = model->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); - - 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::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) { - AKANTU_DEBUG_IN(); - - Mesh::type_iterator it = model->getInterfaceMesh().firstType(); - Mesh::type_iterator last_type = model->getInterfaceMesh().lastType(); - - for (; it != last_type; ++it) { - computeGradU(*it, ghost_type); - computeStress(*it, ghost_type); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialReinforcement::assembleInternalForces( - const ElementType & type, GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - Mesh & mesh = model->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( - const ElementType & interface_type, const ElementType & background_type, - GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - UInt voigt_size = getTangentStiffnessVoigtSize(dim); - - FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); - - Array & elem_filter = 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); - 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); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ - -template -void MaterialReinforcement::computeDirectingCosines( - const ElementType & type, GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - Mesh & interface_mesh = this->model->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") - .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)); - - 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 - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ - -template -void MaterialReinforcement::assembleStiffnessMatrix( - const ElementType & type, GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - Mesh & mesh = model->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( - const ElementType & interface_type, const ElementType & background_type, - GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - UInt voigt_size = getTangentStiffnessVoigtSize(dim); - - FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); - - Array & elem_filter = 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); - - 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"); - 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); - - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ - -/// In this function, type and ghost_type refer to background elements -template -void MaterialReinforcement::computeBackgroundShapeDerivatives( - const ElementType & type, GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - Mesh & interface_mesh = model->getInterfaceMesh(); - - FEEngine & engine = model->getFEEngine(); - FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); - - Mesh::type_iterator interface_type = interface_mesh.firstType(); - Mesh::type_iterator interface_last = interface_mesh.lastType(); - - for (; interface_type != interface_last; ++interface_type) { - Array & filter = element_filter(*interface_type, ghost_type); - - const UInt nb_elements = filter.size(); - const UInt nb_nodes = Mesh::getNbNodesPerElement(type); - const UInt nb_quad_per_element = - interface_engine.getNbIntegrationPoints(*interface_type); - - 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); - - 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); - - Array::scalar_iterator back_it = background_elements->begin(), - back_end = background_elements->end(); - - Array::matrix_iterator shapesd_it = - background_shapesd.begin(dim, nb_nodes); - - Array::vector_iterator quad_pos_it = quad_pos.begin(dim); - - 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); - } - - delete background_elements; - } - - 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); - - 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; - } - - return epot; - } - - 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); - -/* -------------------------------------------------------------------------- */ - -} // akantu 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..8ff027c8f 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,210 @@ /** * @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 filters for background elements + void initFilters(); + /// 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); + /// Add the prestress to the computed stress + void addPrestress(const ElementType & type, GhostType ghost_type); + + /// Compute displacement gradient in reinforcement + void computeGradU(const ElementType & interface_type, GhostType ghost_type); + /// 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, + void filterInterfaceBackgroundElements(Array & foreground, + Array & background, 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); + /// Get background filter + Array & getBackgroundFilter(const ElementType & fg_type, + const ElementType & bg_type, + GhostType ghost_type) { + return (*background_filter(fg_type, ghost_type))(bg_type, ghost_type); + } + + /// Get foreground filter + Array & getForegroundFilter(const ElementType & fg_type, + const ElementType & bg_type, + GhostType ghost_type) { + return (*foreground_filter(fg_type, ghost_type))(bg_type, ghost_type); + } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Embedded model - EmbeddedInterfaceModel * model; - - /// Stress in the reinforcement - InternalField stress_embedded; + EmbeddedInterfaceModel & emodel; /// 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" + /// Foreground mesh filter (contains segment ids) + CrossMap foreground_filter; + + /// Background element filter (contains bg ids) + 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_tmpl.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh new file mode 100644 index 000000000..09caa4037 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh @@ -0,0 +1,818 @@ +/** + * @file material_reinforcement_tmpl.hh + * + * @author Lucas Frerot + * + * @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( + EmbeddedInterfaceModel & model, const ID & id) + : Mat(model, 1, + model.getInterfaceMesh(), + model.getFEEngine("EmbeddedInterfaceFEEngine"), id), + emodel(model), + gradu_embedded("gradu_embedded", *this, 1, + model.getFEEngine("EmbeddedInterfaceFEEngine"), + this->element_filter), + directing_cosines("directing_cosines", *this, 1, + model.getFEEngine("EmbeddedInterfaceFEEngine"), + this->element_filter), + pre_stress("pre_stress", *this, 1, + model.getFEEngine("EmbeddedInterfaceFEEngine"), + this->element_filter), + area(1.0), shape_derivatives() { + AKANTU_DEBUG_IN(); + this->initialize(); + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialReinforcement::initialize() { + AKANTU_DEBUG_IN(); + + 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->emodel.getInterfaceMesh(), + _spatial_dimension = 1); + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ + +template +MaterialReinforcement::~MaterialReinforcement() { + AKANTU_DEBUG_IN(); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ + +template +void MaterialReinforcement::initMaterial() { + Mat::initMaterial(); + + gradu_embedded.initialize(dim * dim); + pre_stress.initialize(1); + + /// We initialise the stuff that is not going to change during the simulation + this->initFilters(); + this->allocBackgroundShapeDerivatives(); + this->initBackgroundShapeDerivatives(); + this->initDirectingCosines(); +} + +/* -------------------------------------------------------------------------- */ + +namespace detail { + class FilterInitializer : public MeshElementTypeMapArrayInializer { + public: + FilterInitializer(EmbeddedInterfaceModel & emodel, + const GhostType & ghost_type) + : MeshElementTypeMapArrayInializer(emodel.getMesh(), + 1, emodel.getSpatialDimension(), + ghost_type, _ek_regular) {} + + UInt size(const ElementType & /*bgtype*/) const override { return 0; } + }; +} + +/* -------------------------------------------------------------------------- */ +/// Initialize the filter for background elements +template +void MaterialReinforcement::initFilters() { + for (auto gt : ghost_types) { + for (auto && type : emodel.getInterfaceMesh().elementTypes(1, gt)) { + std::string shaped_id = "filter"; + if (gt == _ghost) + shaped_id += ":ghost"; + + auto & background = + background_filter(std::make_unique>( + "bg_" + shaped_id, this->name), + type, gt); + auto & foreground = foreground_filter( + std::make_unique>(shaped_id, this->name), + type, gt); + foreground->initialize(detail::FilterInitializer(emodel, gt), 0, true); + background->initialize(detail::FilterInitializer(emodel, gt), 0, true); + + // Computing filters + for (auto && bg_type : background->elementTypes(dim, gt)) { + filterInterfaceBackgroundElements( + (*foreground)(bg_type), (*background)(bg_type), bg_type, type, gt); + } + } + } +} + +/* -------------------------------------------------------------------------- */ +/// Construct a filter for a (interface_type, background_type) pair +template +void MaterialReinforcement::filterInterfaceBackgroundElements( + Array & foreground, Array & background, + const ElementType & type, const ElementType & interface_type, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + foreground.resize(0); + background.resize(0); + + Array & elements = + emodel.getInterfaceAssociatedElements(interface_type, ghost_type); + Array & elem_filter = this->element_filter(interface_type, ghost_type); + + for (auto & elem_id : elem_filter) { + Element & elem = elements(elem_id); + if (elem.type == type) { + background.push_back(elem.element); + foreground.push_back(elem_id); + } + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ + +namespace detail { + class BackgroundShapeDInitializer : public ElementTypeMapArrayInializer { + public: + BackgroundShapeDInitializer(UInt spatial_dimension, + FEEngine & engine, + const ElementType & foreground_type, + ElementTypeMapArray & filter, + const GhostType & ghost_type) + : ElementTypeMapArrayInializer(spatial_dimension, 0, ghost_type, + _ek_regular) { + auto nb_quad = engine.getNbIntegrationPoints(foreground_type); + // Counting 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() * nb_quad, 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 = shape_derivatives( + std::make_unique>(shaped_id, this->name), + type, gt); + shaped_etma->initialize( + detail::BackgroundShapeDInitializer( + emodel.getSpatialDimension(), + emodel.getFEEngine("EmbeddedInterfaceFEEngine"), type, + *background_filter(type, gt), gt), + 0, true); + } + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ + +template +void MaterialReinforcement::initBackgroundShapeDerivatives() { + AKANTU_DEBUG_IN(); + + 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 nb_nodes_elem_bg = Mesh::getNbNodesPerElement(bg_type); + // 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 & foreground_elements = + (*foreground_filter(interface_type, ghost_type))(bg_type, ghost_type); + + auto shapesd_begin = + background_shapesd.begin(dim, nb_nodes_elem_bg, nb_quads_per_elem); + auto quad_begin = quad_pos.begin(dim, nb_quads_per_elem); + + for (auto && tuple : zip(background_elements, foreground_elements)) { + UInt bg = std::get<0>(tuple), fg = std::get<1>(tuple); + for (UInt i = 0; i < nb_quads_per_elem; ++i) { + Matrix shapesd = Tensor3(shapesd_begin[fg])(i); + Vector quads = Matrix(quad_begin[fg])(i); + + engine.computeShapeDerivatives(quads, bg, bg_type, shapesd, + ghost_type); + } + } +} + +/* -------------------------------------------------------------------------- */ + +template +void MaterialReinforcement::initDirectingCosines() { + AKANTU_DEBUG_IN(); + + 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 = VoigtHelper::size; + directing_cosines.initialize(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) { + AKANTU_DEBUG_IN(); + + 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) { + AKANTU_DEBUG_IN(); + + 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::computeAllStresses(GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + Mesh::type_iterator it = emodel.getInterfaceMesh().firstType(); + Mesh::type_iterator last_type = emodel.getInterfaceMesh().lastType(); + + for (; it != last_type; ++it) { + computeGradU(*it, ghost_type); + this->computeStress(*it, ghost_type); + addPrestress(*it, ghost_type); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialReinforcement::addPrestress(const ElementType & type, + GhostType ghost_type) { + auto & stress = this->stress(type, ghost_type); + auto & sigma_p = this->pre_stress(type, ghost_type); + + for (auto && tuple : zip(stress, sigma_p)) { + std::get<0>(tuple) += std::get<1>(tuple); + } +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialReinforcement::assembleInternalForces( + const ElementType & type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + 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( + const ElementType & interface_type, const ElementType & background_type, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + UInt voigt_size = VoigtHelper::size; + + FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); + + 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::vector_iterator C_it = + directing_cosines(interface_type, ghost_type).begin(voigt_size); + + auto sigma_it = this->stress(interface_type, ghost_type).begin(); + + Matrix Bvoigt(voigt_size, back_dof); + + for (; integrant_it != integrant_end; + ++integrant_it, ++B_it, ++C_it, ++sigma_it) { + VoigtHelper::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt, + nodes_per_background_e); + Vector & C = *C_it; + Vector & BtCt_sigma = *integrant_it; + + BtCt_sigma.mul(Bvoigt, C); + BtCt_sigma *= *sigma_it * 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"); + + auto & filter = + getBackgroundFilter(interface_type, background_type, ghost_type); + + emodel.getDOFManager().assembleElementalArrayLocalArray( + residual_interface, emodel.getInternalForce(), background_type, + ghost_type, -1., filter); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ + +template +void MaterialReinforcement::computeDirectingCosines( + const ElementType & type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + 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 = 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->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(1, 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 + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ + +template +void MaterialReinforcement::assembleStiffnessMatrix( + const ElementType & type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + 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( + const ElementType & interface_type, const ElementType & background_type, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + UInt voigt_size = VoigtHelper::size; + + FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); + + 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"); + 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 DCB(1, 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(1, 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); + + DCB.mul(C, Bvoigt); + DCB *= D * area; + CtDCB.mul(C, DCB); + BtCtDCB.mul(Bvoigt, CtDCB); + } + + tangent_moduli.resize(0); + + 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 + + auto & filter = + getBackgroundFilter(interface_type, background_type, ghost_type); + + emodel.getDOFManager().assembleElementalMatricesToMatrix( + "K", "displacement", K_interface, background_type, ghost_type, _symmetric, + filter); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +Real MaterialReinforcement::getEnergy(const std::string & id) { + AKANTU_DEBUG_IN(); + if (id == "potential") { + Real epot = 0.; + + this->computePotentialEnergyByElements(); + + Mesh::type_iterator it = this->element_filter.firstType( + this->spatial_dimension), + end = this->element_filter.lastType( + this->spatial_dimension); + + 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; + } + + return epot; + } + + AKANTU_DEBUG_OUT(); + return 0; +} + +/* -------------------------------------------------------------------------- */ + +template +void MaterialReinforcement::computeGradU( + const ElementType & interface_type, GhostType ghost_type) { + // Looping over background types + for (auto && bg_type : + background_filter(interface_type, ghost_type)->elementTypes(dim)) { + + const UInt nodes_per_background_e = Mesh::getNbNodesPerElement(bg_type); + const UInt voigt_size = VoigtHelper::size; + + auto & bg_shapesd = + (*shape_derivatives(interface_type, ghost_type))(bg_type, ghost_type); + + auto & filter = getBackgroundFilter(interface_type, bg_type, ghost_type); + + Array disp_per_element(0, dim * nodes_per_background_e, "disp_elem"); + + FEEngine::extractNodalToElementField( + emodel.getMesh(), emodel.getDisplacement(), disp_per_element, bg_type, + ghost_type, filter); + + Matrix concrete_du(dim, dim); + Matrix epsilon(dim, dim); + Vector evoigt(voigt_size); + + for (auto && tuple : + zip(make_view(disp_per_element, dim, nodes_per_background_e), + make_view(bg_shapesd, dim, nodes_per_background_e), + this->gradu(interface_type, ghost_type), + make_view(this->directing_cosines(interface_type, ghost_type), + voigt_size))) { + auto & u = std::get<0>(tuple); + auto & B = std::get<1>(tuple); + auto & du = std::get<2>(tuple); + auto & C = std::get<3>(tuple); + + concrete_du.mul(u, B); + auto epsilon = 0.5 * (concrete_du + concrete_du.transpose()); + strainTensorToVoigtVector(epsilon, evoigt); + du = C.dot(evoigt); + } + } +} + +/* -------------------------------------------------------------------------- */ + +/** + * 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; + + Real sq_length = 0.; + for (UInt i = 0; i < dim; i++) { + sq_length += delta(i) * delta(i); + } + + 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 +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(); +} + +} // 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..81df39886 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,178 @@ /** * @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 + SolidMechanicsModel::initModel(); 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, + for_each_element(getInterfaceMesh(), [&](auto && element) { auto mat_index = (*interface_material_selector)(element); - material_index(element) = mat_index; + // material_index(element) = mat_index; materials[mat_index]->addElement(element); + // this->material_local_numbering(element) = index; }, - _element_filter = filter); + _element_filter = filter, + _spatial_dimension = 1); 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 diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index 4ff13730b..dc3e5e8d2 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,441 +1,442 @@ /** * @file structural_mechanics_model.cc * * @author Fabian Barras * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Thu Oct 15 2015 * * @brief Model implementation for StucturalMechanics elements * * @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 "structural_mechanics_model.hh" #include "dof_manager.hh" -#include "sparse_matrix.hh" #include "integrator_gauss.hh" -#include "shape_structural.hh" #include "mesh.hh" +#include "shape_structural.hh" +#include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER +#include "dumpable_inline_impl.hh" #include "dumper_elemental_field.hh" #include "dumper_iohelper_paraview.hh" -#include "dumpable_inline_impl.hh" #include "group_manager_inline_impl.cc" #endif /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model_inline_impl.cc" /* -------------------------------------------------------------------------- */ - namespace akantu { /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) - : Model(mesh, ModelType::_structural_mechanics_model, dim, id, memory_id), time_step(NAN), f_m2a(1.0), - stress("stress", id, memory_id), + : Model(mesh, ModelType::_structural_mechanics_model, dim, id, memory_id), + time_step(NAN), f_m2a(1.0), stress("stress", id, memory_id), element_material("element_material", id, memory_id), set_ID("beam sets", id, memory_id), rotation_matrix("rotation_matices", id, memory_id) { AKANTU_DEBUG_IN(); registerFEEngineObject("StructuralMechanicsFEEngine", mesh, spatial_dimension); if (spatial_dimension == 2) nb_degree_of_freedom = 3; else if (spatial_dimension == 3) nb_degree_of_freedom = 6; else { AKANTU_DEBUG_TO_IMPLEMENT(); } #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("paraview_all", id, true); #endif this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); this->initDOFManager(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) { // <<<< This is the SolidMechanicsModel implementation for future ref >>>> // material_index.initialize(mesh, _element_kind = _ek_not_defined, - // _default_value = UInt(-1), _with_nb_element = true); + // _default_value = UInt(-1), _with_nb_element = + // true); // material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, // _with_nb_element = true); // Model::initFullImpl(options); // // initialize pbc // if (this->pbc_pair.size() != 0) // this->initPBC(); // // initialize the materials // if (this->parser.getLastParsedFile() != "") { // this->instantiateMaterials(); // } // this->initMaterials(); - // this->initBC(*this, *displacement, *displacement_increment, *external_force); + // this->initBC(*this, *displacement, *displacement_increment, + // *external_force); // <<<< END >>>> Model::initFullImpl(options); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFEEngineBoundary() { /// TODO: this function should not be reimplemented /// we're just avoiding a call to Model::initFEEngineBoundary() } /* -------------------------------------------------------------------------- */ // void StructuralMechanicsModel::setTimeStep(Real time_step) { // this->time_step = time_step; // #if defined(AKANTU_USE_IOHELPER) // this->mesh.getDumper().setTimeStep(time_step); // #endif // } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType) { AKANTU_DEBUG_IN(); this->allocNodalField(displacement_rotation, nb_degree_of_freedom, "displacement"); this->allocNodalField(external_force_momentum, nb_degree_of_freedom, "external_force"); this->allocNodalField(internal_force_momentum, nb_degree_of_freedom, "internal_force"); this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs"); auto & dof_manager = this->getDOFManager(); if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *displacement_rotation, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); } if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(velocity, spatial_dimension, "velocity"); this->allocNodalField(acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { - for (auto && type : - mesh.elementTypes(_element_kind = _ek_structural)) { + for (auto && type : mesh.elementTypes(_element_kind = _ek_structural)) { // computeRotationMatrix(type); element_material.alloc(mesh.getNbElement(type), 1, type); } getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); getDOFManager().getMatrix("K").clear(); - for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { + for (auto & type : + mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define ASSEMBLE_STIFFNESS_MATRIX(type) assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); - for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { + for (auto & type : + mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define COMPUTE_STRESS_ON_QUAD(type) computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::updateResidual() { AKANTU_DEBUG_IN(); auto & K = getDOFManager().getMatrix("K"); internal_force_momentum->clear(); K.matVecMul(*displacement_rotation, *internal_force_momentum); *internal_force_momentum *= -1.; getDOFManager().assembleToResidual("displacement", *external_force_momentum, 1.); getDOFManager().assembleToResidual("displacement", *internal_force_momentum, 1.); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeRotationMatrix(const ElementType & type) { Mesh & mesh = getFEEngine().getMesh(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); if (!rotation_matrix.exists(type)) { rotation_matrix.alloc(nb_element, nb_degree_of_freedom * nb_nodes_per_element * - nb_degree_of_freedom * nb_nodes_per_element, + nb_degree_of_freedom * nb_nodes_per_element, type); } else { rotation_matrix(type).resize(nb_element); } rotation_matrix(type).clear(); Array rotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom); rotations.clear(); #define COMPUTE_ROTATION_MATRIX(type) computeRotationMatrix(rotations); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX); #undef COMPUTE_ROTATION_MATRIX auto R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); auto T_it = rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element, nb_degree_of_freedom * nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) { auto & T = *T_it; auto & R = *R_it; for (UInt k = 0; k < nb_nodes_per_element; ++k) { for (UInt i = 0; i < nb_degree_of_freedom; ++i) for (UInt j = 0; j < nb_degree_of_freedom; ++j) T(k * nb_degree_of_freedom + i, k * nb_degree_of_freedom + j) = R(i, j); } } } /* -------------------------------------------------------------------------- */ dumper::Field * StructuralMechanicsModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; dumper::Field * field = NULL; field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * StructuralMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { UInt n; if (spatial_dimension == 2) { n = 2; } else { n = 3; } if (field_name == "displacement") { return mesh.createStridedNodalField(displacement_rotation, group_name, n, 0, padding_flag); } if (field_name == "rotation") { return mesh.createStridedNodalField(displacement_rotation, group_name, nb_degree_of_freedom - n, n, padding_flag); } if (field_name == "force") { return mesh.createStridedNodalField(external_force_momentum, group_name, n, 0, padding_flag); } if (field_name == "momentum") { return mesh.createStridedNodalField(external_force_momentum, group_name, nb_degree_of_freedom - n, n, padding_flag); } if (field_name == "internal_force") { return mesh.createStridedNodalField(internal_force_momentum, group_name, n, 0, padding_flag); } if (field_name == "internal_momentum") { return mesh.createStridedNodalField(internal_force_momentum, group_name, nb_degree_of_freedom - n, n, padding_flag); } return nullptr; } /* -------------------------------------------------------------------------- */ dumper::Field * StructuralMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool, const ElementKind & kind, const std::string &) { dumper::Field * field = NULL; if (field_name == "element_index_by_material") field = mesh.createElementalField( field_name, group_name, this->spatial_dimension, kind); return field; } /* -------------------------------------------------------------------------- */ /* Virtual methods from SolverCallback */ /* -------------------------------------------------------------------------- */ /// get the type of matrix needed MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) { return _symmetric; } /// callback to assemble a Matrix void StructuralMechanicsModel::assembleMatrix(const ID & id) { if (id == "K") assembleStiffnessMatrix(); } /// callback to assemble a lumped Matrix void StructuralMechanicsModel::assembleLumpedMatrix(const ID & /*id*/) {} /// callback to assemble the residual StructuralMechanicsModel::(rhs) void StructuralMechanicsModel::assembleResidual() { this->getDOFManager().assembleToResidual("displacement", *this->external_force_momentum, 1); } /* -------------------------------------------------------------------------- */ /* Virtual methods from MeshEventHandler */ /* -------------------------------------------------------------------------- */ /// function to implement to react on akantu::NewNodesEvent void StructuralMechanicsModel::onNodesAdded(const Array & /*nodes_list*/, const NewNodesEvent & /*event*/) {} /// function to implement to react on akantu::RemovedNodesEvent void StructuralMechanicsModel::onNodesRemoved( const Array & /*nodes_list*/, const Array & /*new_numbering*/, const RemovedNodesEvent & /*event*/) {} /// function to implement to react on akantu::NewElementsEvent void StructuralMechanicsModel::onElementsAdded( const Array & /*elements_list*/, const NewElementsEvent & /*event*/) {} /// function to implement to react on akantu::RemovedElementsEvent void StructuralMechanicsModel::onElementsRemoved( const Array & /*elements_list*/, const ElementTypeMapArray & /*new_numbering*/, const RemovedElementsEvent & /*event*/) {} /// function to implement to react on akantu::ChangedElementsEvent void StructuralMechanicsModel::onElementsChanged( const Array & /*old_elements_list*/, const Array & /*new_elements_list*/, const ElementTypeMapArray & /*new_numbering*/, const ChangedElementsEvent & /*event*/) {} /* -------------------------------------------------------------------------- */ /* Virtual methods from Model */ /* -------------------------------------------------------------------------- */ /// get some default values for derived classes -std::tuple StructuralMechanicsModel::getDefaultSolverID( - const AnalysisMethod & method) { +std::tuple +StructuralMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _static: { return std::make_tuple("static", _tsst_static); } case _implicit_dynamic: { return std::make_tuple("implicit", _tsst_dynamic); } default: return std::make_tuple("unknown", _tsst_not_defined); } } /* ------------------------------------------------------------------------ */ ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_static: { options.non_linear_solver_type = _nls_linear; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ - } // namespace akantu diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc index 6ad371486..be0eeabd1 100644 --- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc +++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc @@ -1,87 +1,96 @@ /** * @file test_embedded_element_matrix.cc * * @author Lucas Frerot * * @date creation: Wed Mar 25 2015 * @date last modification: Wed May 13 2015 * * @brief test of the class EmbeddedInterfaceModel * * @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 "sparse_solver.hh" #include "sparse_matrix_aij.hh" using namespace akantu; int main(int argc, char * argv[]) { debug::setDebugLevel(dblWarning); initialize("embedded_element.dat", argc, argv); constexpr UInt dim = 2; constexpr ElementType type = _segment_2; - const Real height = 0.5; + const Real height = 0.4; Mesh mesh(dim); mesh.read("triangle.msh"); Mesh reinforcement_mesh(dim, "reinforcement_mesh"); auto & nodes = reinforcement_mesh.getNodes(); nodes.push_back(Vector({0, height})); nodes.push_back(Vector({1, height})); reinforcement_mesh.addConnectivityType(type); auto & connectivity = reinforcement_mesh.getConnectivity(type); connectivity.push_back(Vector({0, 1})); Array names_vec(1, 1, "reinforcement", "reinforcement_names"); reinforcement_mesh.registerData("physical_names").alloc(1, 1, type); reinforcement_mesh.getData("physical_names")(type).copy(names_vec); EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim); model.initFull(_analysis_method = _static); if (model.getInterfaceMesh().getNbElement(type) != 1) return EXIT_FAILURE; if (model.getInterfaceMesh().getSpatialDimension() != 2) return EXIT_FAILURE; - model.assembleStiffnessMatrix(); + try { // matrix should be singular + model.solveStep(); + } catch (debug::SingularMatrixException & e) { + std::cerr << "Matrix is singular, relax, everything is fine :)" << std::endl; + } catch (debug::Exception & e) { + std::cerr << "Unexpceted error: " << e.what() << std::endl; + throw e; + } SparseMatrixAIJ & K = dynamic_cast(model.getDOFManager().getMatrix("K")); + K.saveMatrix("stiffness.mtx"); Math::setTolerance(1e-8); // Testing the assembled stiffness matrix if (!Math::are_float_equal(K(0, 0), 1. - height) || !Math::are_float_equal(K(0, 2), height - 1.) || !Math::are_float_equal(K(2, 0), height - 1.) || !Math::are_float_equal(K(2, 2), 1. - height)) return EXIT_FAILURE; return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc index 76a80f295..116108dc1 100644 --- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc +++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc @@ -1,222 +1,224 @@ /** * @file test_embedded_interface_model_prestress.cc * * @author Lucas Frerot * * @date creation: Tue Apr 28 2015 * @date last modification: Thu Oct 15 2015 * * @brief Embedded model test for prestressing (bases on stress norm) * * @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 "embedded_interface_model.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; #define YG 0.483644859 #define I_eq 0.012488874 #define A_eq (1e-2 + 1. / 7. * 1.) /* -------------------------------------------------------------------------- */ struct StressSolution : public BC::Neumann::FromHigherDim { Real M; Real I; Real yg; Real pre_stress; StressSolution(UInt dim, Real M, Real I, Real yg = 0, Real pre_stress = 0) : BC::Neumann::FromHigherDim(Matrix(dim, dim)), M(M), I(I), yg(yg), pre_stress(pre_stress) {} virtual ~StressSolution() {} void operator()(const IntegrationPoint & /*quad_point*/, Vector & dual, const Vector & coord, const Vector & normals) const { UInt dim = coord.size(); if (dim < 2) AKANTU_DEBUG_ERROR("Solution not valid for 1D"); Matrix stress(dim, dim); stress.clear(); stress(0, 0) = this->stress(coord(1)); dual.mul(stress, normals); } inline Real stress(Real height) const { return -M / I * (height - yg) + pre_stress; } inline Real neutral_axis() const { return -I * pre_stress / M + yg; } }; /* -------------------------------------------------------------------------- */ int main (int argc, char * argv[]) { initialize("prestress.dat", argc, argv); debug::setDebugLevel(dblError); Math::setTolerance(1e-6); const UInt dim = 2; /* -------------------------------------------------------------------------- */ Mesh mesh(dim); mesh.read("embedded_mesh_prestress.msh"); - mesh.createGroupsFromMeshData("physical_names"); + // mesh.createGroupsFromMeshData("physical_names"); Mesh reinforcement_mesh(dim, "reinforcement_mesh"); try { reinforcement_mesh.read("embedded_mesh_prestress_reinforcement.msh"); } catch(debug::Exception & e) {} - reinforcement_mesh.createGroupsFromMeshData("physical_names"); + // reinforcement_mesh.createGroupsFromMeshData("physical_names"); EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim); model.initFull(EmbeddedInterfaceModelOptions(_static)); /* -------------------------------------------------------------------------- */ /* Computation of analytical residual */ /* -------------------------------------------------------------------------- */ /* * q = 1000 N/m * L = 20 m * a = 1 m */ - Real steel_area = - model.getMaterial("reinforcement").get("area"); - Real pre_stress = 1e6; + Real steel_area = model.getMaterial("reinforcement").get("area"); + Real pre_stress = model.getMaterial("reinforcement").get("pre_stress"); Real stress_norm = 0.; - StressSolution * concrete_stress = NULL, * steel_stress = NULL; + StressSolution * concrete_stress = nullptr, * steel_stress = nullptr; Real pre_force = pre_stress * steel_area; Real pre_moment = -pre_force * (YG - 0.25); Real neutral_axis = YG - I_eq / A_eq * pre_force / pre_moment; concrete_stress = new StressSolution(dim, pre_moment, 7. * I_eq, YG, -pre_force / (7. * A_eq)); steel_stress = new StressSolution(dim, pre_moment, I_eq, YG, pre_stress - pre_force / A_eq); stress_norm = std::abs(concrete_stress->stress(1)) * (1 - neutral_axis) * 0.5 + std::abs(concrete_stress->stress(0)) * neutral_axis * 0.5 + std::abs(steel_stress->stress(0.25)) * steel_area; model.applyBC(*concrete_stress, "XBlocked"); - ElementGroup & end_node = mesh.getElementGroup("EndNode"); - NodeGroup & end_node_group = end_node.getNodeGroup(); - NodeGroup::const_node_iterator end_node_it = end_node_group.begin(); + auto end_node = *mesh.getElementGroup("EndNode").getNodeGroup().begin(); - Vector end_node_force = model.getForce().begin(dim)[*end_node_it]; + Vector end_node_force = model.getForce().begin(dim)[end_node]; end_node_force(0) += steel_stress->stress(0.25) * steel_area; Array analytical_residual(mesh.getNbNodes(), dim, "analytical_residual"); analytical_residual.copy(model.getForce()); model.getForce().clear(); delete concrete_stress; delete steel_stress; /* -------------------------------------------------------------------------- */ model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "XBlocked"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "YBlocked"); try { model.solveStep(); } catch (debug::Exception e) { std::cerr << e.what() << std::endl; return EXIT_FAILURE; } /* -------------------------------------------------------------------------- */ /* Computation of FEM residual norm */ /* -------------------------------------------------------------------------- */ ElementGroup & xblocked = mesh.getElementGroup("XBlocked"); NodeGroup & boundary_nodes = xblocked.getNodeGroup(); NodeGroup::const_node_iterator nodes_it = boundary_nodes.begin(), nodes_end = boundary_nodes.end(); - Array::vector_iterator com_res = model.getInternalForce().begin(dim); - Array::vector_iterator ana_res = analytical_residual.begin(dim); + model.assembleInternalForces(); + Array residual(mesh.getNbNodes(), dim, "my_residual"); + residual.copy(model.getInternalForce()); + residual -= model.getForce(); + + Array::vector_iterator com_res = residual.begin(dim); Array::vector_iterator position = mesh.getNodes().begin(dim); Real res_sum = 0.; UInt lower_node = -1; UInt upper_node = -1; Real lower_dist = 1; Real upper_dist = 1; for (; nodes_it != nodes_end ; ++nodes_it) { UInt node_number = *nodes_it; const Vector res = com_res[node_number]; - const Vector ana = ana_res[node_number]; const Vector pos = position[node_number]; if (!Math::are_float_equal(pos(1), 0.25)) { if ((std::abs(pos(1) - 0.25) < lower_dist) && (pos(1) < 0.25)) { lower_dist = std::abs(pos(1) - 0.25); lower_node = node_number; } if ((std::abs(pos(1) - 0.25) < upper_dist) && (pos(1) > 0.25)) { upper_dist = std::abs(pos(1) - 0.25); upper_node = node_number; } } for (UInt i = 0 ; i < dim ; i++) { if (!Math::are_float_equal(pos(1), 0.25)) { res_sum += std::abs(res(i)); } } } const Vector upper_res = com_res[upper_node], lower_res = com_res[lower_node]; - const Vector end_node_res = com_res[*end_node_it]; + const Vector end_node_res = com_res[end_node]; Vector delta = upper_res - lower_res; delta *= lower_dist / (upper_dist + lower_dist); Vector concrete_residual = lower_res + delta; Vector steel_residual = end_node_res - concrete_residual; for (UInt i = 0 ; i < dim ; i++) { res_sum += std::abs(concrete_residual(i)); res_sum += std::abs(steel_residual(i)); } Real relative_error = std::abs(res_sum - stress_norm) / stress_norm; - if (relative_error > 1e-3) + if (relative_error > 1e-3) { + std::cerr << "Relative error = " << relative_error << std::endl; return EXIT_FAILURE; + } finalize(); return 0; }