diff --git a/src/common/aka_extern.cc b/src/common/aka_extern.cc index ca6b5a5e7..78ae76cfa 100644 --- a/src/common/aka_extern.cc +++ b/src/common/aka_extern.cc @@ -1,109 +1,109 @@ /** * @file aka_extern.cc * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Thu Nov 19 2015 * * @brief initialisation of all global variables * to insure the order of creation * * @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_common.hh" #include "aka_math.hh" #include "aka_named_argument.hh" #include "aka_random_generator.hh" #include "communication_tag.hh" #include "cppargparse.hh" #include "parser.hh" #include "solid_mechanics_model.hh" #if defined(AKANTU_COHESIVE_ELEMENT) #include "solid_mechanics_model_cohesive.hh" #endif /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #endif namespace akantu { /* -------------------------------------------------------------------------- */ /* error.hpp variables */ /* -------------------------------------------------------------------------- */ namespace debug { /** \todo write function to get this * values from the environment or a config file */ /// standard output for debug messages std::ostream * _akantu_debug_cout = &std::cerr; /// standard output for normal messages std::ostream & _akantu_cout = std::cout; /// parallel context used in debug messages std::string _parallel_context = ""; Debugger debugger; #if defined(AKANTU_DEBUG_TOOLS) DebugElementManager element_manager; #endif } // namespace debug /* -------------------------------------------------------------------------- */ /// list of ghost iterable types ghost_type_t ghost_types(_casper); /* -------------------------------------------------------------------------- */ /// Paser for commandline arguments ::cppargparse::ArgumentParser static_argparser; /// Parser containing the information parsed by the input file given to initFull Parser static_parser; bool Parser::permissive_parser = false; /* -------------------------------------------------------------------------- */ -Real Math::tolerance = std::numeric_limits::epsilon(); +Real Math::tolerance = 1e2 * std::numeric_limits::epsilon(); /* -------------------------------------------------------------------------- */ const UInt _all_dimensions = UInt(-1); /* -------------------------------------------------------------------------- */ const Array empty_filter(0, 1, "empty_filter"); /* -------------------------------------------------------------------------- */ template <> long int RandomGenerator::_seed = 5489u; template <> std::default_random_engine RandomGenerator::generator(5489u); /* -------------------------------------------------------------------------- */ int Tag::max_tag = 0; /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh index 470987bd4..d813498d7 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,440 +1,443 @@ /** * @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. */ 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; + /// explicit copy + void copy(const ElementTypeMapArray & other); + /*! 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 cedbb7f98..659ef91cb 100644 --- a/src/mesh/element_type_map_tmpl.hh +++ b/src/mesh/element_type_map_tmpl.hh @@ -1,662 +1,677 @@ /** * @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 template -inline Stored & ElementTypeMap::operator()( - U && insertee, const SupportType & type, const GhostType & ghost_type) { +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::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 +void ElementTypeMapArray::copy( + const ElementTypeMapArray & other) { + for (auto && ghost_type : ghost_types) { + for (auto && type : + this->elementTypes(_all_dimensions, ghost_types, _ek_not_defined)) { + const auto & array_to_copy = other(type, ghost_type); + auto & array = + this->alloc(0, array_to_copy.getNbComponent(), type, ghost_type); + array.copy(array_to_copy); + } + } +} + /* -------------------------------------------------------------------------- */ 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 +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/mesh/mesh.cc b/src/mesh/mesh.cc index 4091c694f..49cd44069 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,569 +1,571 @@ /** * @file mesh.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Jan 22 2016 * * @brief class handling meshes * * @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_config.hh" /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "group_manager_inline_impl.cc" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "mesh_utils_distribution.hh" #include "node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumper_field.hh" #include "dumper_internal_material_field.hh" #endif /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id, Communicator & communicator) : Memory(id, memory_id), GroupManager(*this, id + ":group_manager", memory_id), nodes_type(0, 1, id + ":nodes_type"), connectivities("connectivities", id, memory_id), ghosts_counters("ghosts_counters", id, memory_id), normals("normals", id, memory_id), spatial_dimension(spatial_dimension), lower_bounds(spatial_dimension, 0.), upper_bounds(spatial_dimension, 0.), size(spatial_dimension, 0.), local_lower_bounds(spatial_dimension, 0.), local_upper_bounds(spatial_dimension, 0.), mesh_data("mesh_data", id, memory_id), communicator(&communicator) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, communicator) { AKANTU_DEBUG_IN(); this->nodes = std::make_shared>(0, spatial_dimension, id + ":coordinates"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id, memory_id) {} /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, Communicator::getStaticCommunicator()) { this->nodes = nodes; this->nb_global_nodes = this->nodes->size(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getBarycenters(Array & barycenter, const ElementType & type, const GhostType & ghost_type) const { barycenter.resize(getNbElement(type, ghost_type)); for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { getBarycenter(Element{type, UInt(std::get<0>(data)), ghost_type}, std::get<1>(data)); } } /* -------------------------------------------------------------------------- */ Mesh & Mesh::initMeshFacets(const ID & id) { AKANTU_DEBUG_IN(); if (!mesh_facets) { mesh_facets = std::make_unique(spatial_dimension, this->nodes, getID() + ":" + id, getMemoryID()); mesh_facets->mesh_parent = this; mesh_facets->is_mesh_facets = true; MeshUtils::buildAllFacets(*this, *mesh_facets, 0); if (mesh.isDistributed()) { mesh_facets->is_distributed = true; mesh_facets->element_synchronizer = std::make_unique( *mesh_facets, mesh.getElementSynchronizer()); } /// transfers the the mesh physical names to the mesh facets if (not this->hasData("physical_names")) { AKANTU_DEBUG_OUT(); return *mesh_facets; } if (not mesh_facets->hasData("physical_names")) { mesh_facets->registerData("physical_names"); } auto & mesh_phys_data = this->getData("physical_names"); auto & phys_data = mesh_facets->getData("physical_names"); phys_data.initialize(*mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); ElementTypeMapArray barycenters(getID(), "temporary_barycenters"); barycenters.initialize(*mesh_facets, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); for (auto && ghost_type : ghost_types) { for (auto && type : barycenters.elementTypes( spatial_dimension - 1, ghost_type)) { mesh_facets->getBarycenters(barycenters(type, ghost_type), type, ghost_type); } } for_each_element( mesh, [&](auto && element) { Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); auto norm_barycenter = barycenter.norm(); auto tolerance = Math::getTolerance(); if (norm_barycenter > tolerance) tolerance *= norm_barycenter; const auto & element_to_facet = mesh_facets->getElementToSubelement( element.type, element.ghost_type); Vector barycenter_facet(spatial_dimension); auto range = enumerate(make_view(barycenters(element.type, element.ghost_type), spatial_dimension)); #ifndef AKANTU_NDEBUG auto min_dist = std::numeric_limits::max(); #endif // this is a spacial search coded the most inefficient way. auto facet = std::find_if(range.begin(), range.end(), [&](auto && data) { auto facet = std::get<0>(data); if (element_to_facet(facet)[1] == ElementNull) return false; auto norm_distance = barycenter.distance(std::get<1>(data)); #ifndef AKANTU_NDEBUG min_dist = std::min(min_dist, norm_distance); #endif return (norm_distance < tolerance); }); if (facet == range.end()) { AKANTU_DEBUG_INFO("The element " << element << " did not find its associated facet in the " "mesh_facets! Try to decrease math tolerance. " "The closest element was at a distance of " << min_dist); return; } // set physical name phys_data(Element{element.type, UInt(std::get<0>(*facet)), element.ghost_type}) = mesh_phys_data(element); }, _spatial_dimension = spatial_dimension - 1); mesh_facets->createGroupsFromMeshData("physical_names"); } AKANTU_DEBUG_OUT(); return *mesh_facets; } /* -------------------------------------------------------------------------- */ void Mesh::defineMeshParent(const Mesh & mesh) { AKANTU_DEBUG_IN(); this->mesh_parent = &mesh; this->is_mesh_facets = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() = default; /* -------------------------------------------------------------------------- */ void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) { MeshIO mesh_io; mesh_io.read(filename, *this, mesh_io_type); type_iterator it = this->firstType(spatial_dimension, _not_ghost, _ek_not_defined); type_iterator last = this->lastType(spatial_dimension, _not_ghost, _ek_not_defined); if (it == last) AKANTU_EXCEPTION( "The mesh contained in the file " << filename << " does not seem to be of the good dimension." << " No element of dimension " << spatial_dimension << " where read."); this->nb_global_nodes = this->nodes->size(); this->computeBoundingBox(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } } /* -------------------------------------------------------------------------- */ void Mesh::write(const std::string & filename, const MeshIOType & mesh_io_type) { MeshIO mesh_io; mesh_io.write(filename, *this, mesh_io_type); } /* -------------------------------------------------------------------------- */ void Mesh::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Mesh [" << std::endl; stream << space << " + id : " << getID() << std::endl; stream << space << " + spatial dimension : " << this->spatial_dimension << std::endl; stream << space << " + nodes [" << std::endl; nodes->printself(stream, indent + 2); stream << space << " + connectivities [" << std::endl; connectivities.printself(stream, indent + 2); stream << space << " ]" << std::endl; GroupManager::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Mesh::computeBoundingBox() { AKANTU_DEBUG_IN(); for (UInt k = 0; k < spatial_dimension; ++k) { local_lower_bounds(k) = std::numeric_limits::max(); local_upper_bounds(k) = -std::numeric_limits::max(); } for (UInt i = 0; i < nodes->size(); ++i) { // if(!isPureGhostNode(i)) for (UInt k = 0; k < spatial_dimension; ++k) { local_lower_bounds(k) = std::min(local_lower_bounds[k], (*nodes)(i, k)); local_upper_bounds(k) = std::max(local_upper_bounds[k], (*nodes)(i, k)); } } if (this->is_distributed) { Matrix reduce_bounds(spatial_dimension, 2); for (UInt k = 0; k < spatial_dimension; ++k) { reduce_bounds(k, 0) = local_lower_bounds(k); reduce_bounds(k, 1) = -local_upper_bounds(k); } communicator->allReduce(reduce_bounds, SynchronizerOperation::_min); for (UInt k = 0; k < spatial_dimension; ++k) { lower_bounds(k) = reduce_bounds(k, 0); upper_bounds(k) = -reduce_bounds(k, 1); } } else { this->lower_bounds = this->local_lower_bounds; this->upper_bounds = this->local_upper_bounds; } size = upper_bounds - lower_bounds; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initNormals() { normals.initialize(*this, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); } /* -------------------------------------------------------------------------- */ void Mesh::getGlobalConnectivity( ElementTypeMapArray & global_connectivity) { AKANTU_DEBUG_IN(); for (auto && ghost_type : ghost_types) { for (auto type : - global_connectivity.elementTypes(_ghost_type = ghost_type)) { + global_connectivity.elementTypes(_spatial_dimension = _all_dimensions, + _element_kind = _ek_not_defined, + _ghost_type = ghost_type)) { if (not connectivities.exists(type, ghost_type)) continue; auto & local_conn = connectivities(type, ghost_type); auto & g_connectivity = global_connectivity(type, ghost_type); UInt nb_nodes = local_conn.size() * local_conn.getNbComponent(); if (not nodes_global_ids && is_mesh_facets) { std::transform( local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), [& node_ids = *mesh_parent->nodes_global_ids](UInt l)->UInt { return node_ids(l); }); } else { std::transform(local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), [& node_ids = *nodes_global_ids](UInt l)->UInt { return node_ids(l); }); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Mesh::getGroupDumper(const std::string & dumper_name, const std::string & group_name) { if (group_name == "all") return this->getDumper(dumper_name); else return element_groups[group_name]->getDumper(dumper_name); } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & arrays, const ElementKind & element_kind) { ElementTypeMap nb_data_per_elem; for (auto type : elementTypes(spatial_dimension, _not_ghost, element_kind)) { UInt nb_elements = this->getNbElement(type); auto & array = arrays(type); nb_data_per_elem(type) = array.getNbComponent() * array.size(); nb_data_per_elem(type) /= nb_elements; } return nb_data_per_elem; } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER template dumper::Field * Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind) { dumper::Field * field = nullptr; ElementTypeMapArray * internal = nullptr; try { internal = &(this->getData(field_id)); } catch (...) { return nullptr; } ElementTypeMap nb_data_per_elem = this->getNbDataPerElem(*internal, element_kind); field = this->createElementalField( *internal, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); return field; } template dumper::Field * Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); template dumper::Field * Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); #endif /* -------------------------------------------------------------------------- */ void Mesh::distribute() { this->distribute(Communicator::getStaticCommunicator()); } /* -------------------------------------------------------------------------- */ void Mesh::distribute(Communicator & communicator) { AKANTU_DEBUG_ASSERT(is_distributed == false, "This mesh is already distribute"); this->communicator = &communicator; this->element_synchronizer = std::make_unique( *this, this->getID() + ":element_synchronizer", this->getMemoryID(), true); this->node_synchronizer = std::make_unique( *this, this->getID() + ":node_synchronizer", this->getMemoryID(), true); Int psize = this->communicator->getNbProc(); #ifdef AKANTU_USE_SCOTCH Int prank = this->communicator->whoAmI(); if (prank == 0) { MeshPartitionScotch partition(*this, spatial_dimension); partition.partitionate(psize); MeshUtilsDistribution::distributeMeshCentralized(*this, 0, partition); } else { MeshUtilsDistribution::distributeMeshCentralized(*this, 0); } #else if (!(psize == 1)) { AKANTU_DEBUG_ERROR("Cannot distribute a mesh without a partitioning tool"); } #endif this->is_distributed = true; this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const Array & node_list, Array & elements) { for (const auto & node : node_list) for (const auto & element : *nodes_to_elements[node]) elements.push_back(element); } /* -------------------------------------------------------------------------- */ void Mesh::fillNodesToElements() { Element e; UInt nb_nodes = nodes->size(); for (UInt n = 0; n < nb_nodes; ++n) { if (this->nodes_to_elements[n]) this->nodes_to_elements[n]->clear(); else this->nodes_to_elements[n] = std::make_unique>(); } for (auto ghost_type : ghost_types) { e.ghost_type = ghost_type; for (const auto & type : elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { e.type = type; UInt nb_element = this->getNbElement(type, ghost_type); Array::const_iterator> conn_it = connectivities(type, ghost_type) .begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; const Vector & conn = *conn_it; for (UInt n = 0; n < conn.size(); ++n) nodes_to_elements[conn(n)]->insert(e); } } } } /* -------------------------------------------------------------------------- */ void Mesh::eraseElements(const Array & elements) { ElementTypeMap last_element; RemovedElementsEvent event(*this); auto & remove_list = event.getList(); auto & new_numbering = event.getNewNumbering(); for (auto && el : elements) { if (el.ghost_type != _not_ghost) { auto & count = ghosts_counters(el); --count; if (count > 0) continue; } remove_list.push_back(el); if (not last_element.exists(el.type, el.ghost_type)) { UInt nb_element = mesh.getNbElement(el.type, el.ghost_type); last_element(nb_element - 1, el.type, el.ghost_type); auto & numbering = new_numbering.alloc(nb_element, 1, el.type, el.ghost_type); for (auto && pair : enumerate(numbering)) { std::get<1>(pair) = std::get<0>(pair); } } UInt & pos = last_element(el.type, el.ghost_type); auto & numbering = new_numbering(el.type, el.ghost_type); numbering(el.element) = UInt(-1); numbering(pos) = el.element; --pos; } this->sendEvent(event); } } // namespace akantu diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index 1cb298d0d..93f8e8090 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,2172 +1,2173 @@ /** * @file mesh_utils.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Leonardo Snozzi * @author Marco Vocialta * * @date creation: Fri Aug 20 2010 * @date last modification: Fri Jan 22 2016 * * @brief All mesh utils necessary for various tasks * * @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 "mesh_utils.hh" #include "element_synchronizer.hh" #include "fe_engine.hh" #include "mesh_accessor.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == _all_dimensions) spatial_dimension = mesh.getSpatialDimension(); /// count number of occurrence of each node UInt nb_nodes = mesh.getNbNodes(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); AKANTU_DEBUG_ASSERT( mesh.firstType(spatial_dimension) != mesh.lastType(spatial_dimension), "Some elements must be found in right dimension to compute facets!"); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt, _ek_not_defined); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_not_defined); for (; first != last; ++first) { ElementType type = *first; UInt nb_element = mesh.getNbElement(type, *gt); Array::const_iterator> conn_it = mesh.getConnectivity(type, *gt).begin( Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) for (UInt n = 0; n < conn_it->size(); ++n) ++node_to_elem.rowOffset((*conn_it)(n)); } } node_to_elem.countToCSR(); node_to_elem.resizeCols(); /// rearrange element to get the node-element list Element e; node_to_elem.beginInsertions(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt, _ek_not_defined); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_not_defined); e.ghost_type = *gt; for (; first != last; ++first) { ElementType type = *first; e.type = type; UInt nb_element = mesh.getNbElement(type, *gt); Array::const_iterator> conn_it = mesh.getConnectivity(type, *gt).begin( Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; for (UInt n = 0; n < conn_it->size(); ++n) node_to_elem.insertInRow((*conn_it)(n), e); } } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * This function should disappear in the future (used in mesh partitioning) */ // void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & // node_to_elem, // UInt spatial_dimension) { // AKANTU_DEBUG_IN(); // if (spatial_dimension == _all_dimensions) // spatial_dimension = mesh.getSpatialDimension(); // UInt nb_nodes = mesh.getNbNodes(); // const Mesh::ConnectivityTypeList & type_list = // mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator // it; // UInt nb_types = type_list.size(); // UInt nb_good_types = 0; // Vector nb_nodes_per_element(nb_types); // UInt ** conn_val = new UInt *[nb_types]; // Vector nb_element(nb_types); // for (it = type_list.begin(); it != type_list.end(); ++it) { // ElementType type = *it; // if (Mesh::getSpatialDimension(type) != spatial_dimension) // continue; // nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type); // conn_val[nb_good_types] = mesh.getConnectivity(type, // _not_ghost).storage(); nb_element[nb_good_types] = // mesh.getConnectivity(type, _not_ghost).size(); // nb_good_types++; // } // AKANTU_DEBUG_ASSERT( // nb_good_types != 0, // "Some elements must be found in right dimension to compute facets!"); // /// array for the node-element list // node_to_elem.resizeRows(nb_nodes); // node_to_elem.clearRows(); // /// count number of occurrence of each node // for (UInt t = 0; t < nb_good_types; ++t) { // for (UInt el = 0; el < nb_element[t]; ++el) { // UInt el_offset = el * nb_nodes_per_element[t]; // for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) { // ++node_to_elem.rowOffset(conn_val[t][el_offset + n]); // } // } // } // node_to_elem.countToCSR(); // node_to_elem.resizeCols(); // node_to_elem.beginInsertions(); // /// rearrange element to get the node-element list // for (UInt t = 0, linearized_el = 0; t < nb_good_types; ++t) // for (UInt el = 0; el < nb_element[t]; ++el, ++linearized_el) { // UInt el_offset = el * nb_nodes_per_element[t]; // for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) // node_to_elem.insertInRow(conn_val[t][el_offset + n], linearized_el); // } // node_to_elem.endInsertions(); // delete[] conn_val; // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsElementTypeMap(const Mesh & mesh, CSR & node_to_elem, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_elements = mesh.getConnectivity(type, ghost_type).size(); UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); /// count number of occurrence of each node for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) ++node_to_elem.rowOffset(conn_val[el_offset + n]); } /// convert the occurrence array in a csr one node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// save the element index in the node-element list for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { node_to_elem.insertInRow(conn_val[el_offset + n], el); } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacets(Mesh & mesh) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh.lastType(spatial_dimension - 1, gt_facet); for (; it != end; ++it) { mesh.getConnectivity(*it, *gt).resize(0); // \todo inform the mesh event handler } } buildFacetsDimension(mesh, mesh, true, spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt to_dimension) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); buildAllFacets(mesh, mesh_facets, spatial_dimension, to_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt from_dimension, UInt to_dimension) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT( mesh_facets.isMeshFacets(), "The mesh_facets should be initialized with initMeshFacets"); /// generate facets buildFacetsDimension(mesh, mesh_facets, false, from_dimension); /// copy nodes type MeshAccessor mesh_accessor_facets(mesh_facets); mesh_accessor_facets.getNodesType().copy(mesh.getNodesType()); /// sort facets and generate sub-facets for (UInt i = from_dimension - 1; i > to_dimension; --i) { buildFacetsDimension(mesh_facets, mesh_facets, false, i); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacetsDimension(const Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension) { AKANTU_DEBUG_IN(); // save the current parent of mesh_facets and set it temporarly to mesh since // mesh is the one containing the elements for which mesh_facets has the // sub-elements // example: if the function is called with mesh = mesh_facets const Mesh * mesh_facets_parent = nullptr; try { mesh_facets_parent = &mesh_facets.getMeshParent(); } catch (...) { } mesh_facets.defineMeshParent(mesh); MeshAccessor mesh_accessor(mesh_facets); UInt spatial_dimension = mesh.getSpatialDimension(); const Array & mesh_facets_nodes = mesh_facets.getNodes(); const Array::const_vector_iterator mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension); CSR node_to_elem; buildNode2Elements(mesh, node_to_elem, dimension); Array counter; std::vector connected_elements; // init the SubelementToElement data to improve performance for (auto && ghost_type : ghost_types) { for (auto && type : mesh.elementTypes(dimension, ghost_type)) { mesh_accessor.getSubelementToElement(type, ghost_type); auto facet_types = mesh.getAllFacetTypes(type); for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); mesh_accessor.getElementToSubelement(facet_type, ghost_type); mesh_accessor.getConnectivity(facet_type, ghost_type); } } } const ElementSynchronizer * synchronizer = nullptr; if (mesh.isDistributed()) { synchronizer = &(mesh.getElementSynchronizer()); } Element current_element; for (auto && ghost_type : ghost_types) { GhostType facet_ghost_type = ghost_type; current_element.ghost_type = ghost_type; for (auto && type : mesh.elementTypes(dimension, ghost_type)) { auto facet_types = mesh.getAllFacetTypes(type); current_element.type = type; for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); auto nb_element = mesh.getNbElement(type, ghost_type); auto element_to_subelement = &mesh_facets.getElementToSubelement(facet_type, ghost_type); auto connectivity_facets = &mesh_facets.getConnectivity(facet_type, ghost_type); auto nb_facet_per_element = mesh.getNbFacetsPerElement(type, ft); const auto & element_connectivity = mesh.getConnectivity(type, ghost_type); Matrix facet_local_connectivity( mesh.getFacetLocalConnectivity(type, ft)); auto nb_nodes_per_facet = connectivity_facets->getNbComponent(); Vector facet(nb_nodes_per_facet); for (UInt el = 0; el < nb_element; ++el) { current_element.element = el; for (UInt f = 0; f < nb_facet_per_element; ++f) { for (UInt n = 0; n < nb_nodes_per_facet; ++n) facet(n) = element_connectivity(el, facet_local_connectivity(f, n)); UInt first_node_nb_elements = node_to_elem.getNbCols(facet(0)); counter.resize(first_node_nb_elements); counter.clear(); // loop over the other nodes to search intersecting elements, // which are the elements that share another node with the // starting element after first_node UInt local_el = 0; auto first_node_elements = node_to_elem.begin(facet(0)); auto first_node_elements_end = node_to_elem.end(facet(0)); for (; first_node_elements != first_node_elements_end; ++first_node_elements, ++local_el) { for (UInt n = 1; n < nb_nodes_per_facet; ++n) { auto node_elements_begin = node_to_elem.begin(facet(n)); auto node_elements_end = node_to_elem.end(facet(n)); counter(local_el) += std::count(node_elements_begin, node_elements_end, *first_node_elements); } } // counting the number of elements connected to the facets and // taking the minimum element number, because the facet should // be inserted just once UInt nb_element_connected_to_facet = 0; Element minimum_el = ElementNull; connected_elements.clear(); for (UInt el_f = 0; el_f < first_node_nb_elements; el_f++) { Element real_el = node_to_elem(facet(0), el_f); if (not(counter(el_f) == nb_nodes_per_facet - 1)) continue; ++nb_element_connected_to_facet; minimum_el = std::min(minimum_el, real_el); connected_elements.push_back(real_el); } if (minimum_el != current_element) continue; bool full_ghost_facet = false; UInt n = 0; while (n < nb_nodes_per_facet && mesh.isPureGhostNode(facet(n))) ++n; if (n == nb_nodes_per_facet) full_ghost_facet = true; if (full_ghost_facet) continue; if (boundary_only and nb_element_connected_to_facet != 1) continue; std::vector elements; // build elements_on_facets: linearized_el must come first // in order to store the facet in the correct direction // and avoid to invert the sign in the normal computation elements.push_back(current_element); if (nb_element_connected_to_facet == 1) { /// boundary facet elements.push_back(ElementNull); } else if (nb_element_connected_to_facet == 2) { /// internal facet elements.push_back(connected_elements[1]); /// check if facet is in between ghost and normal /// elements: if it's the case, the facet is either /// ghost or not ghost. The criterion to decide this /// is arbitrary. It was chosen to check the processor /// id (prank) of the two neighboring elements. If /// prank of the ghost element is lower than prank of /// the normal one, the facet is not ghost, otherwise /// it's ghost GhostType gt[2] = {_not_ghost, _not_ghost}; for (UInt el = 0; el < connected_elements.size(); ++el) gt[el] = connected_elements[el].ghost_type; if (gt[0] != gt[1] and (gt[0] == _not_ghost or gt[1] == _not_ghost)) { UInt prank[2]; for (UInt el = 0; el < 2; ++el) { prank[el] = synchronizer->getRank(connected_elements[el]); } // ugly trick from Marco detected :P bool ghost_one = (gt[0] != _ghost); if (prank[ghost_one] > prank[!ghost_one]) facet_ghost_type = _not_ghost; else facet_ghost_type = _ghost; connectivity_facets = &mesh_facets.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_facets.getElementToSubelement( facet_type, facet_ghost_type); } } else { /// facet of facet for (UInt i = 1; i < nb_element_connected_to_facet; ++i) { elements.push_back(connected_elements[i]); } } element_to_subelement->push_back(elements); connectivity_facets->push_back(facet); /// current facet index UInt current_facet = connectivity_facets->size() - 1; /// loop on every element connected to current facet and /// insert current facet in the first free spot of the /// subelement_to_element vector for (UInt elem = 0; elem < elements.size(); ++elem) { Element loc_el = elements[elem]; if (loc_el.type == _not_defined) continue; Array & subelement_to_element = mesh_facets.getSubelementToElement(loc_el.type, loc_el.ghost_type); UInt nb_facet_per_loc_element = subelement_to_element.getNbComponent(); for (UInt f_in = 0; f_in < nb_facet_per_loc_element; ++f_in) { auto & el = subelement_to_element(loc_el.element, f_in); if (el.type != _not_defined) continue; el.type = facet_type; el.element = current_facet; el.ghost_type = facet_ghost_type; break; } } /// reset connectivity in case a facet was found in /// between ghost and normal elements if (facet_ghost_type != ghost_type) { facet_ghost_type = ghost_type; connectivity_facets = &mesh_accessor.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_accessor.getElementToSubelement( facet_type, facet_ghost_type); } } } } } } // restore the parent of mesh_facet if (mesh_facets_parent) mesh_facets.defineMeshParent(*mesh_facets_parent); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, Array & local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Array & old_nodes_numbers) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::map renumbering_map; for (UInt i = 0; i < old_nodes_numbers.size(); ++i) { renumbering_map[old_nodes_numbers(i)] = i; } /// renumber the nodes renumberNodesInConnectivity(local_connectivities, (nb_local_element + nb_ghost_element) * nb_nodes_per_element, renumbering_map); old_nodes_numbers.resize(renumbering_map.size()); for (auto & renumber_pair : renumbering_map) { old_nodes_numbers(renumber_pair.second) = renumber_pair.first; } renumbering_map.clear(); MeshAccessor mesh_accessor(mesh); /// copy the renumbered connectivity to the right place auto & local_conn = mesh_accessor.getConnectivity(type); local_conn.resize(nb_local_element); memcpy(local_conn.storage(), local_connectivities.storage(), nb_local_element * nb_nodes_per_element * sizeof(UInt)); auto & ghost_conn = mesh_accessor.getConnectivity(type, _ghost); ghost_conn.resize(nb_ghost_element); std::memcpy(ghost_conn.storage(), local_connectivities.storage() + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); auto & ghost_counter = mesh_accessor.getGhostsCounters(type, _ghost); ghost_counter.resize(nb_ghost_element, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberNodesInConnectivity( Array & list_nodes, UInt nb_nodes, std::map & renumbering_map) { AKANTU_DEBUG_IN(); UInt * connectivity = list_nodes.storage(); UInt new_node_num = renumbering_map.size(); for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) { UInt & node = *connectivity; auto it = renumbering_map.find(node); if (it == renumbering_map.end()) { UInt old_node = node; renumbering_map[old_node] = new_node_num; node = new_node_num; ++new_node_num; } else { node = it->second; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::purifyMesh(Mesh & mesh) { AKANTU_DEBUG_IN(); std::map renumbering_map; RemovedNodesEvent remove_nodes(mesh); Array & nodes_removed = remove_nodes.getList(); for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { auto ghost_type = (GhostType)gt; Mesh::type_iterator it = mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined); for (; it != end; ++it) { ElementType type(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); Array & connectivity = mesh.getConnectivity(type, ghost_type); UInt nb_element(connectivity.size()); renumberNodesInConnectivity( connectivity, nb_element * nb_nodes_per_element, renumbering_map); } } Array & new_numbering = remove_nodes.getNewNumbering(); std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1)); auto it = renumbering_map.begin(); auto end = renumbering_map.end(); for (; it != end; ++it) { new_numbering(it->first) = it->second; } for (UInt i = 0; i < new_numbering.size(); ++i) { if (new_numbering(i) == UInt(-1)) nodes_removed.push_back(i); } mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } #if defined(AKANTU_COHESIVE_ELEMENT) /* -------------------------------------------------------------------------- */ UInt MeshUtils::insertCohesiveElements( Mesh & mesh, Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion, Array & doubled_nodes, Array & new_elements, bool only_double_facets) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt elements_to_insert = updateFacetToDouble(mesh_facets, facet_insertion); if (elements_to_insert > 0) { if (spatial_dimension == 1) { doublePointFacet(mesh, mesh_facets, doubled_nodes); } else { doubleFacet(mesh, mesh_facets, spatial_dimension - 1, doubled_nodes, true); findSubfacetToDouble(mesh, mesh_facets); if (spatial_dimension == 2) { doubleSubfacet<2>(mesh, mesh_facets, doubled_nodes); } else if (spatial_dimension == 3) { doubleFacet(mesh, mesh_facets, 1, doubled_nodes, false); findSubfacetToDouble(mesh, mesh_facets); doubleSubfacet<3>(mesh, mesh_facets, doubled_nodes); } } if (!only_double_facets) updateCohesiveData(mesh, mesh_facets, new_elements); } return elements_to_insert; } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doubleNodes(Mesh & mesh, const std::vector & old_nodes, Array & doubled_nodes) { AKANTU_DEBUG_IN(); Array & position = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt old_nb_nodes = position.size(); UInt new_nb_nodes = old_nb_nodes + old_nodes.size(); UInt old_nb_doubled_nodes = doubled_nodes.size(); UInt new_nb_doubled_nodes = old_nb_doubled_nodes + old_nodes.size(); position.resize(new_nb_nodes); doubled_nodes.resize(new_nb_doubled_nodes); Array::iterator> position_begin = position.begin(spatial_dimension); for (UInt n = 0; n < old_nodes.size(); ++n) { UInt new_node = old_nb_nodes + n; /// store doubled nodes doubled_nodes(old_nb_doubled_nodes + n, 0) = old_nodes[n]; doubled_nodes(old_nb_doubled_nodes + n, 1) = new_node; /// update position std::copy(position_begin + old_nodes[n], position_begin + old_nodes[n] + 1, position_begin + new_node); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::doubleFacet(Mesh & mesh, Mesh & mesh_facets, UInt facet_dimension, Array & doubled_nodes, bool facet_mode) { AKANTU_DEBUG_IN(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(facet_dimension, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(facet_dimension, gt_facet); for (; it != end; ++it) { ElementType type_facet = *it; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); if (nb_facet_to_double == 0) continue; ElementType type_subfacet = Mesh::getFacetType(type_facet); const UInt nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(type_facet); GhostType gt_subfacet = _casper; Array> * f_to_subfacet = nullptr; Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); UInt nb_nodes_per_facet = conn_facet.getNbComponent(); UInt old_nb_facet = conn_facet.size(); UInt new_nb_facet = old_nb_facet + nb_facet_to_double; conn_facet.resize(new_nb_facet); subfacet_to_facet.resize(new_nb_facet); UInt new_facet = old_nb_facet - 1; Element new_facet_el{type_facet, 0, gt_facet}; Array::iterator> subfacet_to_facet_begin = subfacet_to_facet.begin(nb_subfacet_per_facet); Array::iterator> conn_facet_begin = conn_facet.begin(nb_nodes_per_facet); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); ++new_facet; /// adding a new facet by copying original one /// copy connectivity in new facet std::copy(conn_facet_begin + old_facet, conn_facet_begin + old_facet + 1, conn_facet_begin + new_facet); /// update subfacet_to_facet std::copy(subfacet_to_facet_begin + old_facet, subfacet_to_facet_begin + old_facet + 1, subfacet_to_facet_begin + new_facet); new_facet_el.element = new_facet; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet_per_facet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; if (gt_subfacet != subfacet.ghost_type) { gt_subfacet = subfacet.ghost_type; f_to_subfacet = &mesh_facets.getElementToSubelement( type_subfacet, subfacet.ghost_type); } /// update facet_to_subfacet array (*f_to_subfacet)(subfacet.element).push_back(new_facet_el); } } /// update facet_to_subfacet and _segment_3 facets if any if (!facet_mode) { updateSubfacetToFacet(mesh_facets, type_facet, gt_facet, true); updateFacetToSubfacet(mesh_facets, type_facet, gt_facet, true); updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } else updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt MeshUtils::updateFacetToDouble( Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); UInt nb_facets_to_double = 0.; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for (; it != end; ++it) { ElementType type_facet = *it; const Array & f_insertion = facet_insertion(type_facet, gt_facet); Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); ElementType el_type = _not_defined; GhostType el_gt = _casper; UInt nb_facet_per_element = 0; Element old_facet_el{type_facet, 0, gt_facet}; Array * facet_to_element = nullptr; for (UInt f = 0; f < f_insertion.size(); ++f) { if (f_insertion(f) == false) continue; ++nb_facets_to_double; if (element_to_facet(f)[1].type == _not_defined #if defined(AKANTU_COHESIVE_ELEMENT) || element_to_facet(f)[1].kind() == _ek_cohesive #endif ) { AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary"); continue; } f_to_double.push_back(f); UInt new_facet = mesh_facets.getNbElement(type_facet, gt_facet) + f_to_double.size() - 1; old_facet_el.element = f; /// update facet_to_element vector Element & elem_to_update = element_to_facet(f)[1]; UInt el = elem_to_update.element; if (elem_to_update.ghost_type != el_gt || elem_to_update.type != el_type) { el_type = elem_to_update.type; el_gt = elem_to_update.ghost_type; facet_to_element = &mesh_facets.getSubelementToElement(el_type, el_gt); nb_facet_per_element = facet_to_element->getNbComponent(); } Element * f_update = std::find(facet_to_element->storage() + el * nb_facet_per_element, facet_to_element->storage() + el * nb_facet_per_element + nb_facet_per_element, old_facet_el); AKANTU_DEBUG_ASSERT( facet_to_element->storage() + el * nb_facet_per_element != facet_to_element->storage() + el * nb_facet_per_element + nb_facet_per_element, "Facet not found"); f_update->element = new_facet; /// update elements connected to facet std::vector first_facet_list = element_to_facet(f); element_to_facet.push_back(first_facet_list); /// set new and original facets as boundary facets element_to_facet(new_facet)[0] = element_to_facet(new_facet)[1]; element_to_facet(f)[1] = ElementNull; element_to_facet(new_facet)[1] = ElementNull; } } } AKANTU_DEBUG_OUT(); return nb_facets_to_double; } /* -------------------------------------------------------------------------- */ void MeshUtils::resetFacetToDouble(Mesh & mesh_facets) { AKANTU_DEBUG_IN(); for (UInt g = _not_ghost; g <= _ghost; ++g) { auto gt = (GhostType)g; Mesh::type_iterator it = mesh_facets.firstType(_all_dimensions, gt); Mesh::type_iterator end = mesh_facets.lastType(_all_dimensions, gt); for (; it != end; ++it) { ElementType type = *it; mesh_facets.getDataPointer("facet_to_double", type, gt, 1, false); mesh_facets.getDataPointer>( "facets_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer>( "elements_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer>( "subfacets_to_subsubfacet_double", type, gt, 1, false); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::findSubfacetToDouble(Mesh & mesh, Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); if (spatial_dimension == 1) return; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for (; it != end; ++it) { ElementType type_facet = *it; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); if (nb_facet_to_double == 0) continue; ElementType type_subfacet = Mesh::getFacetType(type_facet); GhostType gt_subfacet = _casper; ElementType type_subsubfacet = Mesh::getFacetType(type_subfacet); GhostType gt_subsubfacet = _casper; Array * conn_subfacet = nullptr; Array * sf_to_double = nullptr; Array> * sf_to_subfacet_double = nullptr; Array> * f_to_subfacet_double = nullptr; Array> * el_to_subfacet_double = nullptr; UInt nb_subfacet = Mesh::getNbFacetsPerElement(type_facet); UInt nb_subsubfacet; UInt nb_nodes_per_sf_el; if (subsubfacet_mode) { nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subsubfacet); nb_subsubfacet = Mesh::getNbFacetsPerElement(type_subfacet); } else nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subfacet); Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); Array * subsubfacet_to_subfacet = nullptr; UInt old_nb_facet = subfacet_to_facet.size() - nb_facet_to_double; Element current_facet{type_facet, 0, gt_facet}; std::vector element_list; std::vector facet_list; std::vector * subfacet_list; if (subsubfacet_mode) subfacet_list = new std::vector; /// map to filter subfacets Array> * facet_to_subfacet = nullptr; /// this is used to make sure that both new and old facets are /// checked UInt facets[2]; /// loop on every facet for (UInt f_index = 0; f_index < 2; ++f_index) { for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { facets[bool(f_index)] = f_to_double(facet); facets[!bool(f_index)] = old_nb_facet + facet; UInt old_facet = facets[0]; UInt new_facet = facets[1]; Element & starting_element = element_to_facet(new_facet)[0]; current_facet.element = old_facet; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; if (gt_subfacet != subfacet.ghost_type) { gt_subfacet = subfacet.ghost_type; if (subsubfacet_mode) { subsubfacet_to_subfacet = &mesh_facets.getSubelementToElement( type_subfacet, gt_subfacet); } else { conn_subfacet = &mesh_facets.getConnectivity(type_subfacet, gt_subfacet); sf_to_double = &mesh_facets.getData( "facet_to_double", type_subfacet, gt_subfacet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_subfacet, gt_subfacet); facet_to_subfacet = &mesh_facets.getElementToSubelement( type_subfacet, gt_subfacet); } } if (subsubfacet_mode) { /// loop on every subsubfacet for (UInt ssf = 0; ssf < nb_subsubfacet; ++ssf) { Element & subsubfacet = (*subsubfacet_to_subfacet)(subfacet.element, ssf); if (subsubfacet == ElementNull) continue; if (gt_subsubfacet != subsubfacet.ghost_type) { gt_subsubfacet = subsubfacet.ghost_type; conn_subfacet = &mesh_facets.getConnectivity(type_subsubfacet, gt_subsubfacet); sf_to_double = &mesh_facets.getData( "facet_to_double", type_subsubfacet, gt_subsubfacet); sf_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subsubfacet, gt_subsubfacet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subsubfacet, gt_subsubfacet); el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_subsubfacet, gt_subsubfacet); facet_to_subfacet = &mesh_facets.getElementToSubelement( type_subsubfacet, gt_subsubfacet); } UInt global_ssf = subsubfacet.element; Vector subsubfacet_connectivity( conn_subfacet->storage() + global_ssf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subsubfacet is to be doubled if (findElementsAroundSubfacet( mesh, mesh_facets, starting_element, current_facet, subsubfacet_connectivity, element_list, facet_list, subfacet_list) == false && removeElementsInVector(*subfacet_list, (*facet_to_subfacet)(global_ssf)) == false) { sf_to_double->push_back(global_ssf); sf_to_subfacet_double->push_back(*subfacet_list); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } else { const UInt global_sf = subfacet.element; Vector subfacet_connectivity( conn_subfacet->storage() + global_sf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subfacet is to be doubled if (findElementsAroundSubfacet( mesh, mesh_facets, starting_element, current_facet, subfacet_connectivity, element_list, facet_list) == false && removeElementsInVector( facet_list, (*facet_to_subfacet)(global_sf)) == false) { sf_to_double->push_back(global_sf); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } } } if (subsubfacet_mode) delete subfacet_list; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) void MeshUtils::updateCohesiveData(Mesh & mesh, Mesh & mesh_facets, Array & new_elements) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); bool third_dimension = spatial_dimension == 3; MeshAccessor mesh_facets_accessor(mesh_facets); for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); if (nb_facet_to_double == 0) continue; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); auto & facet_to_coh_element = mesh_facets_accessor.getSubelementToElement(type_cohesive, gt_facet); auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto & conn_cohesive = mesh.getConnectivity(type_cohesive, gt_facet); UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(type_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); UInt old_nb_cohesive_elements = conn_cohesive.size(); UInt new_nb_cohesive_elements = conn_cohesive.size() + nb_facet_to_double; UInt old_nb_facet = element_to_facet.size() - nb_facet_to_double; facet_to_coh_element.resize(new_nb_cohesive_elements); conn_cohesive.resize(new_nb_cohesive_elements); UInt new_elements_old_size = new_elements.size(); new_elements.resize(new_elements_old_size + nb_facet_to_double); Element c_element{type_cohesive, 0, gt_facet}; Element f_element{type_facet, 0, gt_facet}; UInt facets[2]; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { /// (in 3D cohesive elements connectivity is inverted) facets[third_dimension] = f_to_double(facet); facets[!third_dimension] = old_nb_facet + facet; UInt cohesive_element = old_nb_cohesive_elements + facet; /// store doubled facets f_element.element = facets[0]; facet_to_coh_element(cohesive_element, 0) = f_element; f_element.element = facets[1]; facet_to_coh_element(cohesive_element, 1) = f_element; /// modify cohesive elements' connectivity for (UInt n = 0; n < nb_nodes_per_facet; ++n) { conn_cohesive(cohesive_element, n) = conn_facet(facets[0], n); conn_cohesive(cohesive_element, n + nb_nodes_per_facet) = conn_facet(facets[1], n); } /// update element_to_facet vectors c_element.element = cohesive_element; element_to_facet(facets[0])[1] = c_element; element_to_facet(facets[1])[1] = c_element; /// add cohesive element to the element event list new_elements(new_elements_old_size + facet) = c_element; } } } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doublePointFacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); if (spatial_dimension != 1) return; Array & position = mesh.getNodes(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for (; it != end; ++it) { ElementType type_facet = *it; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); const Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); UInt old_nb_facet = element_to_facet.size() - nb_facet_to_double; UInt new_nb_facet = element_to_facet.size(); UInt old_nb_nodes = position.size(); UInt new_nb_nodes = old_nb_nodes + nb_facet_to_double; position.resize(new_nb_nodes); conn_facet.resize(new_nb_facet); UInt old_nb_doubled_nodes = doubled_nodes.size(); doubled_nodes.resize(old_nb_doubled_nodes + nb_facet_to_double); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt new_facet = old_nb_facet + facet; ElementType type = element_to_facet(new_facet)[0].type; UInt el = element_to_facet(new_facet)[0].element; GhostType gt = element_to_facet(new_facet)[0].ghost_type; UInt old_node = conn_facet(old_facet); UInt new_node = old_nb_nodes + facet; /// update position position(new_node) = position(old_node); conn_facet(new_facet) = new_node; Array & conn_segment = mesh.getConnectivity(type, gt); UInt nb_nodes_per_segment = conn_segment.getNbComponent(); /// update facet connectivity UInt i; for (i = 0; conn_segment(el, i) != old_node && i <= nb_nodes_per_segment; ++i) ; conn_segment(el, i) = new_node; doubled_nodes(old_nb_doubled_nodes + facet, 0) = old_node; doubled_nodes(old_nb_doubled_nodes + facet, 1) = new_node; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::updateQuadraticSegments(Mesh & mesh, Mesh & mesh_facets, ElementType type_facet, GhostType gt_facet, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (type_facet != _segment_3) return; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); UInt old_nb_facet = mesh_facets.getNbElement(type_facet, gt_facet) - nb_facet_to_double; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); /// this ones matter only for segments in 3D Array> * el_to_subfacet_double = nullptr; Array> * f_to_subfacet_double = nullptr; if (third_dim_segments) { el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_facet, gt_facet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_facet, gt_facet); } std::vector middle_nodes; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt node = conn_facet(old_facet, 2); if (!mesh.isPureGhostNode(node)) middle_nodes.push_back(node); } UInt n = doubled_nodes.size(); doubleNodes(mesh, middle_nodes, doubled_nodes); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt old_node = conn_facet(old_facet, 2); if (mesh.isPureGhostNode(old_node)) continue; UInt new_node = doubled_nodes(n, 1); UInt new_facet = old_nb_facet + facet; conn_facet(new_facet, 2) = new_node; if (third_dim_segments) { updateElementalConnectivity(mesh_facets, old_node, new_node, element_to_facet(new_facet)); updateElementalConnectivity(mesh, old_node, new_node, (*el_to_subfacet_double)(facet), &(*f_to_subfacet_double)(facet)); } else { updateElementalConnectivity(mesh, old_node, new_node, element_to_facet(new_facet)); } ++n; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateSubfacetToFacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); /// update subfacet_to_facet vector ElementType type_facet = _not_defined; GhostType gt_facet = _casper; Array * subfacet_to_facet = nullptr; UInt nb_subfacet_per_facet = 0; UInt old_nb_subfacet = mesh_facets.getNbElement(type_subfacet, gt_subfacet) - nb_subfacet_to_double; Array> * facet_list = nullptr; if (facet_mode) facet_list = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); else facet_list = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); Element old_subfacet_el{type_subfacet, 0, gt_subfacet}; Element new_subfacet_el{type_subfacet, 0, gt_subfacet}; for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { old_subfacet_el.element = sf_to_double(sf); new_subfacet_el.element = old_nb_subfacet + sf; for (UInt f = 0; f < (*facet_list)(sf).size(); ++f) { Element & facet = (*facet_list)(sf)[f]; if (facet.type != type_facet || facet.ghost_type != gt_facet) { type_facet = facet.type; gt_facet = facet.ghost_type; subfacet_to_facet = &mesh_facets.getSubelementToElement(type_facet, gt_facet); nb_subfacet_per_facet = subfacet_to_facet->getNbComponent(); } Element * sf_update = std::find( subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet, subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, old_subfacet_el); AKANTU_DEBUG_ASSERT(subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet != subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, "Subfacet not found"); *sf_update = new_subfacet_el; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateFacetToSubfacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); Array> & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet); Array> * facet_to_subfacet_double = nullptr; if (facet_mode) { facet_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); } else { facet_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); } UInt old_nb_subfacet = facet_to_subfacet.size(); facet_to_subfacet.resize(old_nb_subfacet + nb_subfacet_to_double); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) facet_to_subfacet(old_nb_subfacet + sf) = (*facet_to_subfacet_double)(sf); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::doubleSubfacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (spatial_dimension == 1) return; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_subfacet = *gt; Mesh::type_iterator it = mesh_facets.firstType(0, gt_subfacet); Mesh::type_iterator end = mesh_facets.lastType(0, gt_subfacet); for (; it != end; ++it) { ElementType type_subfacet = *it; Array & sf_to_double = mesh_facets.getData( "facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); if (nb_subfacet_to_double == 0) continue; AKANTU_DEBUG_ASSERT( type_subfacet == _point_1, "Only _point_1 subfacet doubling is supported at the moment"); Array & conn_subfacet = mesh_facets.getConnectivity(type_subfacet, gt_subfacet); UInt old_nb_subfacet = conn_subfacet.size(); UInt new_nb_subfacet = old_nb_subfacet + nb_subfacet_to_double; conn_subfacet.resize(new_nb_subfacet); std::vector nodes_to_double; UInt old_nb_doubled_nodes = doubled_nodes.size(); /// double nodes for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_subfacet = sf_to_double(sf); nodes_to_double.push_back(conn_subfacet(old_subfacet)); } doubleNodes(mesh, nodes_to_double, doubled_nodes); /// add new nodes in connectivity for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt new_subfacet = old_nb_subfacet + sf; UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); conn_subfacet(new_subfacet) = new_node; } /// update facet and element connectivity Array> & f_to_subfacet_double = mesh_facets.getData>("facets_to_subfacet_double", type_subfacet, gt_subfacet); Array> & el_to_subfacet_double = mesh_facets.getData>( "elements_to_subfacet_double", type_subfacet, gt_subfacet); Array> * sf_to_subfacet_double = nullptr; if (spatial_dimension == 3) sf_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_node = doubled_nodes(old_nb_doubled_nodes + sf, 0); UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); updateElementalConnectivity(mesh, old_node, new_node, el_to_subfacet_double(sf), &f_to_subfacet_double(sf)); updateElementalConnectivity(mesh_facets, old_node, new_node, f_to_subfacet_double(sf)); if (spatial_dimension == 3) updateElementalConnectivity(mesh_facets, old_node, new_node, (*sf_to_subfacet_double)(sf)); } if (spatial_dimension == 2) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, true); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, true); } else if (spatial_dimension == 3) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, false); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::flipFacets( Mesh & mesh_facets, const ElementTypeMapArray & global_connectivity, GhostType gt_facet) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); /// get global connectivity for local mesh ElementTypeMapArray global_connectivity_tmp("global_connectivity_tmp", mesh_facets.getID(), mesh_facets.getMemoryID()); global_connectivity_tmp.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, + _ghost_type = gt_facet, _with_nb_nodes_per_element = true, _with_nb_element = true); mesh_facets.getGlobalConnectivity(global_connectivity_tmp); /// loop on every facet for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & connectivity = mesh_facets.getConnectivity(type_facet, gt_facet); const auto & g_connectivity = global_connectivity(type_facet, gt_facet); auto & el_to_f = mesh_facets.getElementToSubelement(type_facet, gt_facet); auto & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); UInt nb_subfacet_per_facet = subfacet_to_facet.getNbComponent(); UInt nb_nodes_per_facet = connectivity.getNbComponent(); UInt nb_facet = connectivity.size(); UInt nb_nodes_per_P1_facet = Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type_facet)); auto & global_conn_tmp = global_connectivity_tmp(type_facet, gt_facet); auto conn_it = connectivity.begin(nb_nodes_per_facet); auto gconn_tmp_it = global_conn_tmp.begin(nb_nodes_per_facet); auto conn_glob_it = g_connectivity.begin(nb_nodes_per_facet); auto subf_to_f = subfacet_to_facet.begin(nb_subfacet_per_facet); auto * conn_tmp_pointer = new UInt[nb_nodes_per_facet]; Vector conn_tmp(conn_tmp_pointer, nb_nodes_per_facet); Element el_tmp; auto * subf_tmp_pointer = new Element[nb_subfacet_per_facet]; Vector subf_tmp(subf_tmp_pointer, nb_subfacet_per_facet); for (UInt f = 0; f < nb_facet; ++f, ++conn_it, ++gconn_tmp_it, ++subf_to_f, ++conn_glob_it) { Vector & gconn_tmp = *gconn_tmp_it; const Vector & conn_glob = *conn_glob_it; Vector & conn_local = *conn_it; /// skip facet if connectivities are the same if (gconn_tmp == conn_glob) continue; /// re-arrange connectivity conn_tmp = conn_local; UInt * begin = conn_glob.storage(); UInt * end = conn_glob.storage() + nb_nodes_per_facet; for (UInt n = 0; n < nb_nodes_per_facet; ++n) { UInt * new_node = std::find(begin, end, gconn_tmp(n)); AKANTU_DEBUG_ASSERT(new_node != end, "Node not found"); UInt new_position = new_node - begin; conn_local(new_position) = conn_tmp(n); } /// if 3D, check if facets are just rotated if (spatial_dimension == 3) { /// find first node UInt * new_node = std::find(begin, end, gconn_tmp(0)); AKANTU_DEBUG_ASSERT(new_node != end, "Node not found"); UInt new_position = new_node - begin; UInt n = 1; /// count how many nodes in the received connectivity follow /// the same order of those in the local connectivity for (; n < nb_nodes_per_P1_facet && gconn_tmp(n) == conn_glob((new_position + n) % nb_nodes_per_P1_facet); ++n) ; /// skip the facet inversion if facet is just rotated if (n == nb_nodes_per_P1_facet) continue; } /// update data to invert facet el_tmp = el_to_f(f)[0]; el_to_f(f)[0] = el_to_f(f)[1]; el_to_f(f)[1] = el_tmp; subf_tmp = (*subf_to_f); (*subf_to_f)(0) = subf_tmp(1); (*subf_to_f)(1) = subf_tmp(0); } delete[] subf_tmp_pointer; delete[] conn_tmp_pointer; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::fillElementToSubElementsData(Mesh & mesh) { AKANTU_DEBUG_IN(); if (mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) { AKANTU_DEBUG_INFO("There are not facets, add them in the mesh file or call " "the buildFacet method."); return; } UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray barycenters("barycenter_tmp", mesh.getID(), mesh.getMemoryID()); barycenters.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = _all_dimensions); // mesh.initElementTypeMapArray(barycenters, spatial_dimension, // _all_dimensions); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (auto & type : mesh.elementTypes(_all_dimensions, ghost_type)) { element.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & barycenters_arr = barycenters(type, ghost_type); barycenters_arr.resize(nb_element); auto bary = barycenters_arr.begin(spatial_dimension); auto bary_end = barycenters_arr.end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { element.element = el; mesh.getBarycenter(element, *bary); } } } MeshAccessor mesh_accessor(mesh); for (Int sp(spatial_dimension); sp >= 1; --sp) { if (mesh.getNbElement(sp) == 0) continue; for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(sp, ghost_type)) { mesh_accessor.getSubelementToElement(type, ghost_type) .resize(mesh.getNbElement(type, ghost_type)); mesh_accessor.getSubelementToElement(type, ghost_type).clear(); } for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) { mesh_accessor.getElementToSubelement(type, ghost_type) .resize(mesh.getNbElement(type, ghost_type)); mesh.getElementToSubelement(type, ghost_type).clear(); } } CSR nodes_to_elements; buildNode2Elements(mesh, nodes_to_elements, sp); Element facet_element; for (auto ghost_type : ghost_types) { facet_element.ghost_type = ghost_type; for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) { facet_element.type = type; auto & element_to_subelement = mesh.getElementToSubelement(type, ghost_type); const auto & connectivity = mesh.getConnectivity(type, ghost_type); auto fit = connectivity.begin(mesh.getNbNodesPerElement(type)); auto fend = connectivity.end(mesh.getNbNodesPerElement(type)); UInt fid = 0; for (; fit != fend; ++fit, ++fid) { const Vector & facet = *fit; facet_element.element = fid; std::map element_seen_counter; UInt nb_nodes_per_facet = mesh.getNbNodesPerElement(Mesh::getP1ElementType(type)); for (UInt n(0); n < nb_nodes_per_facet; ++n) { auto eit = nodes_to_elements.begin(facet(n)); auto eend = nodes_to_elements.end(facet(n)); for (; eit != eend; ++eit) { auto & elem = *eit; auto cit = element_seen_counter.find(elem); if (cit != element_seen_counter.end()) { cit->second++; } else { element_seen_counter[elem] = 1; } } } std::vector connected_elements; auto cit = element_seen_counter.begin(); auto cend = element_seen_counter.end(); for (; cit != cend; ++cit) { if (cit->second == nb_nodes_per_facet) connected_elements.push_back(cit->first); } auto ceit = connected_elements.begin(); auto ceend = connected_elements.end(); for (; ceit != ceend; ++ceit) element_to_subelement(fid).push_back(*ceit); for (UInt ce = 0; ce < connected_elements.size(); ++ce) { Element & elem = connected_elements[ce]; Array & subelement_to_element = mesh_accessor.getSubelementToElement(elem.type, elem.ghost_type); UInt f(0); for (; f < mesh.getNbFacetsPerElement(elem.type) && subelement_to_element(elem.element, f) != ElementNull; ++f) ; AKANTU_DEBUG_ASSERT( f < mesh.getNbFacetsPerElement(elem.type), "The element " << elem << " seems to have too many facets!! (" << f << " < " << mesh.getNbFacetsPerElement(elem.type) << ")"); subelement_to_element(elem.element, f) = facet_element; } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MeshUtils::findElementsAroundSubfacet( const Mesh & mesh, const Mesh & mesh_facets, const Element & starting_element, const Element & end_facet, const Vector & subfacet_connectivity, std::vector & elem_list, std::vector & facet_list, std::vector * subfacet_list) { AKANTU_DEBUG_IN(); /// preallocated stuff before starting bool facet_matched = false; elem_list.clear(); facet_list.clear(); if (third_dim_points) subfacet_list->clear(); elem_list.push_back(starting_element); const Array * facet_connectivity = nullptr; const Array * sf_connectivity = nullptr; const Array * facet_to_element = nullptr; const Array * subfacet_to_facet = nullptr; ElementType current_type = _not_defined; GhostType current_ghost_type = _casper; ElementType current_facet_type = _not_defined; GhostType current_facet_ghost_type = _casper; ElementType current_subfacet_type = _not_defined; GhostType current_subfacet_ghost_type = _casper; const Array> * element_to_facet = nullptr; const Element * opposing_el = nullptr; std::queue elements_to_check; elements_to_check.push(starting_element); /// keep going until there are elements to check while (!elements_to_check.empty()) { /// check current element Element & current_el = elements_to_check.front(); if (current_el.type != current_type || current_el.ghost_type != current_ghost_type) { current_type = current_el.type; current_ghost_type = current_el.ghost_type; facet_to_element = &mesh_facets.getSubelementToElement(current_type, current_ghost_type); } /// loop over each facet of the element for (UInt f = 0; f < facet_to_element->getNbComponent(); ++f) { const Element & current_facet = (*facet_to_element)(current_el.element, f); if (current_facet == ElementNull) continue; if (current_facet_type != current_facet.type || current_facet_ghost_type != current_facet.ghost_type) { current_facet_type = current_facet.type; current_facet_ghost_type = current_facet.ghost_type; element_to_facet = &mesh_facets.getElementToSubelement( current_facet_type, current_facet_ghost_type); facet_connectivity = &mesh_facets.getConnectivity( current_facet_type, current_facet_ghost_type); if (third_dim_points) subfacet_to_facet = &mesh_facets.getSubelementToElement( current_facet_type, current_facet_ghost_type); } /// check if end facet is reached if (current_facet == end_facet) facet_matched = true; /// add this facet if not already passed if (std::find(facet_list.begin(), facet_list.end(), current_facet) == facet_list.end() && hasElement(*facet_connectivity, current_facet, subfacet_connectivity)) { facet_list.push_back(current_facet); if (third_dim_points) { /// check subfacets for (UInt sf = 0; sf < subfacet_to_facet->getNbComponent(); ++sf) { const Element & current_subfacet = (*subfacet_to_facet)(current_facet.element, sf); if (current_subfacet == ElementNull) continue; if (current_subfacet_type != current_subfacet.type || current_subfacet_ghost_type != current_subfacet.ghost_type) { current_subfacet_type = current_subfacet.type; current_subfacet_ghost_type = current_subfacet.ghost_type; sf_connectivity = &mesh_facets.getConnectivity( current_subfacet_type, current_subfacet_ghost_type); } if (std::find(subfacet_list->begin(), subfacet_list->end(), current_subfacet) == subfacet_list->end() && hasElement(*sf_connectivity, current_subfacet, subfacet_connectivity)) subfacet_list->push_back(current_subfacet); } } } else continue; /// consider opposing element if ((*element_to_facet)(current_facet.element)[0] == current_el) opposing_el = &(*element_to_facet)(current_facet.element)[1]; else opposing_el = &(*element_to_facet)(current_facet.element)[0]; /// skip null elements since they are on a boundary if (*opposing_el == ElementNull) continue; /// skip this element if already added if (std::find(elem_list.begin(), elem_list.end(), *opposing_el) != elem_list.end()) continue; /// only regular elements have to be checked if (opposing_el->kind() == _ek_regular) elements_to_check.push(*opposing_el); elem_list.push_back(*opposing_el); #ifndef AKANTU_NDEBUG const Array & conn_elem = mesh.getConnectivity(opposing_el->type, opposing_el->ghost_type); AKANTU_DEBUG_ASSERT( hasElement(conn_elem, *opposing_el, subfacet_connectivity), "Subfacet doesn't belong to this element"); #endif } /// erased checked element from the list elements_to_check.pop(); } AKANTU_DEBUG_OUT(); return facet_matched; } /* -------------------------------------------------------------------------- */ UInt MeshUtils::updateLocalMasterGlobalConnectivity(Mesh & mesh, UInt local_nb_new_nodes) { const auto & comm = mesh.getCommunicator(); Int rank = comm.whoAmI(); Int nb_proc = comm.getNbProc(); if (nb_proc == 1) return local_nb_new_nodes; /// resize global ids array Array & nodes_global_ids = mesh.getGlobalNodesIds(); UInt old_nb_nodes = mesh.getNbNodes() - local_nb_new_nodes; nodes_global_ids.resize(mesh.getNbNodes()); /// compute the number of global nodes based on the number of old nodes Vector old_local_master_nodes(nb_proc); for (UInt n = 0; n < old_nb_nodes; ++n) if (mesh.isLocalOrMasterNode(n)) ++old_local_master_nodes(rank); comm.allGather(old_local_master_nodes); UInt old_global_nodes = std::accumulate(old_local_master_nodes.storage(), old_local_master_nodes.storage() + nb_proc, 0); /// compute amount of local or master doubled nodes Vector local_master_nodes(nb_proc); for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) if (mesh.isLocalOrMasterNode(n)) ++local_master_nodes(rank); comm.allGather(local_master_nodes); /// update global number of nodes UInt total_nb_new_nodes = std::accumulate( local_master_nodes.storage(), local_master_nodes.storage() + nb_proc, 0); if (total_nb_new_nodes == 0) return 0; /// set global ids of local and master nodes UInt starting_index = std::accumulate(local_master_nodes.storage(), local_master_nodes.storage() + rank, old_global_nodes); for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) { if (mesh.isLocalOrMasterNode(n)) { nodes_global_ids(n) = starting_index; ++starting_index; } } MeshAccessor mesh_accessor(mesh); mesh_accessor.setNbGlobalNodes(old_global_nodes + total_nb_new_nodes); return total_nb_new_nodes; } /* -------------------------------------------------------------------------- */ void MeshUtils::updateElementalConnectivity( Mesh & mesh, UInt old_node, UInt new_node, const std::vector & element_list, __attribute__((unused)) const std::vector * facet_list) { AKANTU_DEBUG_IN(); ElementType el_type = _not_defined; GhostType gt_type = _casper; Array * conn_elem = nullptr; #if defined(AKANTU_COHESIVE_ELEMENT) const Array * cohesive_facets = nullptr; #endif UInt nb_nodes_per_element = 0; UInt * n_update = nullptr; for (UInt el = 0; el < element_list.size(); ++el) { const Element & elem = element_list[el]; if (elem.type == _not_defined) continue; if (elem.type != el_type || elem.ghost_type != gt_type) { el_type = elem.type; gt_type = elem.ghost_type; conn_elem = &mesh.getConnectivity(el_type, gt_type); nb_nodes_per_element = conn_elem->getNbComponent(); #if defined(AKANTU_COHESIVE_ELEMENT) if (elem.kind() == _ek_cohesive) cohesive_facets = &mesh.getMeshFacets().getSubelementToElement(el_type, gt_type); #endif } #if defined(AKANTU_COHESIVE_ELEMENT) if (elem.kind() == _ek_cohesive) { AKANTU_DEBUG_ASSERT( facet_list != nullptr, "Provide a facet list in order to update cohesive elements"); /// loop over cohesive element's facets for (UInt f = 0, n = 0; f < 2; ++f, n += nb_nodes_per_element / 2) { const Element & facet = (*cohesive_facets)(elem.element, f); /// skip facets if not present in the list if (std::find(facet_list->begin(), facet_list->end(), facet) == facet_list->end()) continue; n_update = std::find( conn_elem->storage() + elem.element * nb_nodes_per_element + n, conn_elem->storage() + elem.element * nb_nodes_per_element + n + nb_nodes_per_element / 2, old_node); AKANTU_DEBUG_ASSERT(n_update != conn_elem->storage() + elem.element * nb_nodes_per_element + n + nb_nodes_per_element / 2, "Node not found in current element"); /// update connectivity *n_update = new_node; } } else { #endif n_update = std::find(conn_elem->storage() + elem.element * nb_nodes_per_element, conn_elem->storage() + elem.element * nb_nodes_per_element + nb_nodes_per_element, old_node); AKANTU_DEBUG_ASSERT(n_update != conn_elem->storage() + elem.element * nb_nodes_per_element + nb_nodes_per_element, "Node not found in current element"); /// update connectivity *n_update = new_node; #if defined(AKANTU_COHESIVE_ELEMENT) } #endif } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index 51cbada6c..16058656f 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1553 +1,1554 @@ /** * @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(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("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/material_inline_impl.cc b/src/model/solid_mechanics/material_inline_impl.cc index f6141480e..4b06b4d13 100644 --- a/src/model/solid_mechanics/material_inline_impl.cc +++ b/src/model/solid_mechanics/material_inline_impl.cc @@ -1,463 +1,464 @@ /** * @file material_inline_impl.cc * * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Nov 25 2015 * * @brief Implementation of the inline functions of the class 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 . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_INLINE_IMPL_CC__ #define __AKANTU_MATERIAL_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const ElementType & type, UInt element, const GhostType & ghost_type) { Array & el_filter = this->element_filter(type, ghost_type); el_filter.push_back(element); return el_filter.size() - 1; } /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const Element & element) { return this->addElement(element.type, element.element, element.ghost_type); } /* -------------------------------------------------------------------------- */ inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const { return (dim * (dim - 1) / 2 + dim); } /* -------------------------------------------------------------------------- */ inline UInt Material::getCauchyStressMatrixSize(UInt dim) const { return (dim * dim); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToF(const Matrix & grad_u, Matrix & F) { AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim * dim, "The dimension of the tensor F should be greater or " "equal to the dimension of the tensor grad_u."); F.eye(); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) F(i, j) += grad_u(i, j); } /* -------------------------------------------------------------------------- */ template inline void Material::computeCauchyStressOnQuad(const Matrix & F, const Matrix & piola, Matrix & sigma, const Real & C33) const { Real J = F.det() * sqrt(C33); Matrix F_S(dim, dim); F_S.mul(F, piola); Real constant = J ? 1. / J : 0; sigma.mul(F_S, F, constant); } /* -------------------------------------------------------------------------- */ inline void Material::rightCauchy(const Matrix & F, Matrix & C) { C.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::leftCauchy(const Matrix & F, Matrix & B) { B.mul(F, F); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i)); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToGreenStrain(const Matrix & grad_u, Matrix & epsilon) { epsilon.mul(grad_u, grad_u, .5); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) epsilon(i, j) += 0.5 * (grad_u(i, j) + grad_u(j, i)); } /* -------------------------------------------------------------------------- */ inline Real Material::stressToVonMises(const Matrix & stress) { // compute deviatoric stress UInt dim = stress.cols(); Matrix deviatoric_stress = Matrix::eye(dim, -1. * stress.trace() / 3.); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) deviatoric_stress(i, j) += stress(i, j); // return Von Mises stress return std::sqrt(3. * deviatoric_stress.doubleDot(deviatoric_stress) / 2.); } /* ---------------------------------------------------------------------------*/ template inline void Material::setCauchyStressArray(const Matrix & S_t, Matrix & sigma_voight) { AKANTU_DEBUG_IN(); sigma_voight.clear(); // see Finite element formulations for large deformation dynamic analysis, // Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau /* * 1d: [ s11 ]' * 2d: [ s11 s22 s12 ]' * 3d: [ s11 s22 s33 s23 s13 s12 ] */ for (UInt i = 0; i < dim; ++i) // diagonal terms sigma_voight(i, 0) = S_t(i, i); for (UInt i = 1; i < dim; ++i) // term s12 in 2D and terms s23 s13 in 3D sigma_voight(dim + i - 1, 0) = S_t(dim - i - 1, dim - 1); for (UInt i = 2; i < dim; ++i) // term s13 in 3D sigma_voight(dim + i, 0) = S_t(0, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void Material::setCauchyStressMatrix(const Matrix & S_t, Matrix & sigma) { AKANTU_DEBUG_IN(); sigma.clear(); /// see Finite ekement formulations for large deformation dynamic analysis, /// Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau for (UInt i = 0; i < dim; ++i) { for (UInt m = 0; m < dim; ++m) { for (UInt n = 0; n < dim; ++n) { sigma(i * dim + m, i * dim + n) = S_t(m, n); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline Element Material::convertToLocalElement(const Element & global_element) const { UInt ge = global_element.element; #ifndef AKANTU_NDEBUG UInt model_mat_index = this->model.getMaterialByElement( global_element.type, global_element.ghost_type)(ge); UInt mat_index = this->model.getMaterialIndex(this->name); AKANTU_DEBUG_ASSERT(model_mat_index == mat_index, "Conversion of a global element in a local element for " "the wrong material " << this->name << std::endl); #endif UInt le = this->model.getMaterialLocalNumbering( global_element.type, global_element.ghost_type)(ge); Element tmp_quad{global_element.type, le, global_element.ghost_type}; return tmp_quad; } /* -------------------------------------------------------------------------- */ inline Element Material::convertToGlobalElement(const Element & local_element) const { UInt le = local_element.element; UInt ge = this->element_filter(local_element.type, local_element.ghost_type)(le); Element tmp_quad{local_element.type, ge, local_element.ghost_type}; return tmp_quad; } /* -------------------------------------------------------------------------- */ inline IntegrationPoint Material::convertToLocalPoint(const IntegrationPoint & global_point) const { const FEEngine & fem = this->model.getFEEngine(); UInt nb_quad = fem.getNbIntegrationPoints(global_point.type); Element el = this->convertToLocalElement(static_cast(global_point)); IntegrationPoint tmp_quad(el, global_point.num_point, nb_quad); return tmp_quad; } /* -------------------------------------------------------------------------- */ inline IntegrationPoint Material::convertToGlobalPoint(const IntegrationPoint & local_point) const { const FEEngine & fem = this->model.getFEEngine(); UInt nb_quad = fem.getNbIntegrationPoints(local_point.type); Element el = this->convertToGlobalElement(static_cast(local_point)); IntegrationPoint tmp_quad(el, local_point.num_point, nb_quad); return tmp_quad; } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbData(const Array & elements, const SynchronizationTag & tag) const { if (tag == _gst_smm_stress) { return (this->isFiniteDeformation() ? 3 : 1) * spatial_dimension * spatial_dimension * sizeof(Real) * this->getModel().getNbIntegrationPoints(elements); } return 0; } /* -------------------------------------------------------------------------- */ inline void Material::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { if (tag == _gst_smm_stress) { if (this->isFiniteDeformation()) { packElementDataHelper(piola_kirchhoff_2, buffer, elements); packElementDataHelper(gradu, buffer, elements); } packElementDataHelper(stress, buffer, elements); } } /* -------------------------------------------------------------------------- */ inline void Material::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { if (tag == _gst_smm_stress) { if (this->isFiniteDeformation()) { unpackElementDataHelper(piola_kirchhoff_2, buffer, elements); unpackElementDataHelper(gradu, buffer, elements); } unpackElementDataHelper(stress, buffer, elements); } } /* -------------------------------------------------------------------------- */ inline const Parameter & Material::getParam(const ID & param) const { try { return get(param); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID()); } } /* -------------------------------------------------------------------------- */ template inline void Material::setParam(const ID & param, T value) { try { set(param, value); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID()); } updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template inline void Material::packElementDataHelper( const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) const { DataAccessor::packElementalDataHelper(data_to_pack, buffer, elements, true, model.getFEEngine(fem_id)); } /* -------------------------------------------------------------------------- */ template inline void Material::unpackElementDataHelper( ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) { DataAccessor::unpackElementalDataHelper(data_to_unpack, buffer, elements, true, model.getFEEngine(fem_id)); } /* -------------------------------------------------------------------------- */ template <> inline void Material::registerInternal(InternalField & vect) { internal_vectors_real[vect.getID()] = &vect; } template <> inline void Material::registerInternal(InternalField & vect) { internal_vectors_uint[vect.getID()] = &vect; } template <> inline void Material::registerInternal(InternalField & vect) { internal_vectors_bool[vect.getID()] = &vect; } /* -------------------------------------------------------------------------- */ template <> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_real.erase(vect.getID()); } template <> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_uint.erase(vect.getID()); } template <> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_bool.erase(vect.getID()); } /* -------------------------------------------------------------------------- */ template inline bool Material::isInternal(__attribute__((unused)) const ID & id, - __attribute__((unused)) const ElementKind & element_kind) const { + __attribute__((unused)) + const ElementKind & element_kind) const { AKANTU_DEBUG_TO_IMPLEMENT(); } template <> inline bool Material::isInternal(const ID & id, const ElementKind & element_kind) const { auto internal_array = internal_vectors_real.find(this->getID() + ":" + id); if (internal_array == internal_vectors_real.end() || internal_array->second->getElementKind() != element_kind) return false; return true; } /* -------------------------------------------------------------------------- */ template inline ElementTypeMap Material::getInternalDataPerElem(const ID & field_id, const ElementKind & element_kind) const { if (!this->template isInternal(field_id, element_kind)) AKANTU_EXCEPTION("Cannot find internal field " << id << " in material " << this->name); const InternalField & internal_field = this->template getInternal(field_id); const FEEngine & fe_engine = internal_field.getFEEngine(); UInt nb_data_per_quad = internal_field.getNbComponent(); ElementTypeMap res; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { using type_iterator = typename InternalField::type_iterator; type_iterator tit = internal_field.firstType(*gt); type_iterator tend = internal_field.lastType(*gt); for (; tit != tend; ++tit) { UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(*tit, *gt); res(*tit, *gt) = nb_data_per_quad * nb_quadrature_points; } } return res; } /* -------------------------------------------------------------------------- */ template void Material::flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type, ElementKind element_kind) const { if (!this->template isInternal(field_id, element_kind)) AKANTU_EXCEPTION("Cannot find internal field " << id << " in material " << this->name); const InternalField & internal_field = this->template getInternal(field_id); const FEEngine & fe_engine = internal_field.getFEEngine(); const Mesh & mesh = fe_engine.getMesh(); using type_iterator = typename InternalField::filter_type_iterator; type_iterator tit = internal_field.filterFirstType(ghost_type); type_iterator tend = internal_field.filterLastType(ghost_type); for (; tit != tend; ++tit) { ElementType type = *tit; const Array & src_vect = internal_field(type, ghost_type); const Array & filter = internal_field.getFilter(type, ghost_type); // total number of elements in the corresponding mesh UInt nb_element_dst = mesh.getNbElement(type, ghost_type); // number of element in the internal field UInt nb_element_src = filter.size(); // number of quadrature points per elem UInt nb_quad_per_elem = fe_engine.getNbIntegrationPoints(type); // number of data per quadrature point UInt nb_data_per_quad = internal_field.getNbComponent(); if (!internal_flat.exists(type, ghost_type)) { internal_flat.alloc(nb_element_dst * nb_quad_per_elem, nb_data_per_quad, type, ghost_type); } if (nb_element_src == 0) continue; // number of data per element UInt nb_data = nb_quad_per_elem * nb_data_per_quad; Array & dst_vect = internal_flat(type, ghost_type); dst_vect.resize(nb_element_dst * nb_quad_per_elem); Array::const_scalar_iterator it = filter.begin(); Array::const_scalar_iterator end = filter.end(); Array::const_vector_iterator it_src = src_vect.begin_reinterpret(nb_data, nb_element_src); Array::vector_iterator it_dst = dst_vect.begin_reinterpret(nb_data, nb_element_dst); for (; it != end; ++it, ++it_src) { it_dst[*it] = *it_src; } } } } // akantu #endif /* __AKANTU_MATERIAL_INLINE_IMPL_CC__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc index 9369e9df5..ab54a0dc6 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc @@ -1,742 +1,713 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Jan 13 2016 * * @brief Solid mechanics model for cohesive 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 "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "cohesive_element_inserter.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "fe_engine_template.hh" #include "integrator_gauss.hh" #include "material_cohesive.hh" #include "parser.hh" #include "shape_cohesive.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { -const SolidMechanicsModelCohesiveOptions - default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass, - false); - /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id, ModelType::_solid_mechanics_model_cohesive), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->material_selector); this->material_selector = tmp_material_selector; #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, _gst_material_id); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_stress); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_boundary); } this->inserter = std::make_unique( this->mesh, id + ":cohesive_element_inserter"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step, const ID & solver_id) { SolidMechanicsModel::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); const auto & smmc_options = dynamic_cast(options); this->is_extrinsic = smmc_options.extrinsic; inserter->setIsExtrinsic(is_extrinsic); if (mesh.isDistributed()) { auto & mesh_facets = inserter->getMeshFacets(); auto & synchronizer = dynamic_cast(mesh_facets.getElementSynchronizer()); this->registerSynchronizer(synchronizer, _gst_smmc_facets); this->registerSynchronizer(synchronizer, _gst_smmc_facets_conn); synchronizeGhostFacetsConnectivity(); /// create the facet synchronizer for extrinsic simulations if (is_extrinsic) { - facet_stress_synchronizer = - std::make_unique(mesh_facets); - - facet_stress_synchronizer->updateSchemes( - [&](auto & scheme, auto & proc, auto & direction) { - scheme.copy(const_cast(synchronizer) - .getCommunications() - .getScheme(proc, direction)); - }); + facet_stress_synchronizer = std::make_unique( + synchronizer, id + ":facet_stress_synchronizer"); + this->registerSynchronizer(*facet_stress_synchronizer, _gst_smmc_facets_stress); } inserter->initParallel(*cohesive_synchronizer); } ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { auto inserter_section = section.getSubSections(ParserType::_cohesive_inserter); if (inserter_section.begin() != inserter_section.end()) { inserter->parseSection(*inserter_section.begin()); } } SolidMechanicsModel::initFullImpl(options); AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); /// find the first cohesive material UInt cohesive_index = UInt(-1); for (auto && material : enumerate(materials)) { if (dynamic_cast(std::get<1>(material).get())) { cohesive_index = std::get<0>(material); break; } } if (cohesive_index == UInt(-1)) AKANTU_EXCEPTION("No cohesive materials in the material input file"); material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion // to know what material to call for stress checks if (is_extrinsic) { this->initExtrinsicMaterials(); } else { this->initIntrinsicMaterials(); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initExtrinsicMaterials() { const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = material_selector->getFallbackValue()); - for_each_element(mesh_facets, - [&](auto && element) { - auto mat_index = (*material_selector)(element); - auto & mat = dynamic_cast( - *materials[mat_index]); - facet_material(element) = mat_index; - mat.addFacet(element); - }, - _spatial_dimension = spatial_dimension - 1); + for_each_element( + mesh_facets, + [&](auto && element) { + auto mat_index = (*material_selector)(element); + auto & mat = dynamic_cast(*materials[mat_index]); + facet_material(element) = mat_index; + mat.addFacet(element); + }, + _spatial_dimension = spatial_dimension - 1, _ghost_type = _not_ghost); SolidMechanicsModel::initMaterials(); this->initAutomaticInsertion(); } /* -------------------------------------------------------------------------- */ #if 0 void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials( const std::string & cohesive_surfaces) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != nullptr) inserter->initParallel(facet_synchronizer, cohesive_element_synchronizer); // inserter->initParallel(facet_synchronizer, synch_parallel); #endif std::istringstream split(cohesive_surfaces); std::string physname; while (std::getline(split, physname, ',')) { AKANTU_DEBUG_INFO( "Pre-inserting cohesive elements along facets from physical surface: " << physname); insertElementsFromMeshData(physname); } synchronizeInsertionData(); SolidMechanicsModel::initMaterials(); auto && tmp = std::make_shared(*this); tmp->setFallback(material_selector->getFallbackValue()); tmp->setFallback(material_selector->getFallbackSelector()); material_selector = tmp; // UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::synchronizeInsertionData() { if (inserter->getMeshFacets().isDistributed()) { inserter->getMeshFacets().getElementSynchronizer().synchronizeOnce( *inserter, _gst_ce_groups); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicMaterials() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initMaterials(); this->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); /// add cohesive type connectivity ElementType type = _not_defined; for (auto && type_ghost : ghost_types) { for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost); if (connectivity.size() == 0) continue; type = tmp_type; auto type_facet = Mesh::getFacetType(type); auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModelCohesive::updateFacetStressSynchronizer() { - if (facet_stress_synchronizer != nullptr) { - const auto & rank_to_element = - mesh.getElementSynchronizer().getElementToRank(); - const auto & facet_checks = inserter->getCheckFacets(); - const auto & element_to_facet = mesh.getElementToSubelement(); - UInt rank = mesh.getCommunicator().whoAmI(); - - facet_stress_synchronizer->updateSchemes( - [&](auto & scheme, auto & proc, auto & /*direction*/) { - UInt el = 0; - for (auto && element : scheme) { - if (not facet_checks(element)) - continue; - - const auto & next_el = element_to_facet(element); - UInt rank_left = rank_to_element(next_el[0]); - UInt rank_right = rank_to_element(next_el[1]); - - if ((rank_left == rank and rank_right == proc) or - (rank_left == proc and rank_right == rank)) { - scheme[el] = element; - ++el; - } - } - scheme.resize(el); - }); - } -} - /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); - this->updateFacetStressSynchronizer(); - this->resizeFacetStress(); /// compute normals on facets this->computeNormals(); this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension, // Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); // mesh.initElementTypeMapArray(elements_quad_facets, // Model::spatial_dimension, // Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) continue; /// compute elements' quadrature points and list of facet /// quadrature points positions by element - Array & facet_to_element = + const auto & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); - UInt nb_facet_per_elem = facet_to_element.getNbComponent(); - Array & el_q_facet = elements_quad_facets(type, elem_gt); - ElementType facet_type = Mesh::getFacetType(type); - UInt nb_quad_per_facet = - getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); + auto & el_q_facet = elements_quad_facets(type, elem_gt); - el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet); - - for (UInt el = 0; el < nb_element; ++el) { - for (UInt f = 0; f < nb_facet_per_elem; ++f) { - Element global_facet_elem = facet_to_element(el, f); - UInt global_facet = global_facet_elem.element; - GhostType facet_gt = global_facet_elem.ghost_type; - const Array & quad_f = quad_facets(facet_type, facet_gt); - - for (UInt q = 0; q < nb_quad_per_facet; ++q) { - for (UInt s = 0; s < Model::spatial_dimension; ++s) { - el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + - f * nb_quad_per_facet + q, - s) = quad_f(global_facet * nb_quad_per_facet + q, s); - } - } - } + auto facet_type = Mesh::getFacetType(type); + auto nb_quad_per_facet = + getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); + auto nb_facet_per_elem = facet_to_element.getNbComponent(); + + // small hack in the loop to skip boundary elements, they are silently + // initialized to NaN to see if this causes problems + el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet, + std::numeric_limits::quiet_NaN()); + + for (auto && data : + zip(make_view(facet_to_element), + make_view(el_q_facet, spatial_dimension, nb_quad_per_facet))) { + const auto & global_facet = std::get<0>(data); + auto & el_q = std::get<1>(data); + + if (global_facet == ElementNull) + continue; + + Matrix quad_f = + make_view(quad_facets(global_facet.type, global_facet.ghost_type), + spatial_dimension, nb_quad_per_facet) + .begin()[global_facet.element]; + + el_q = quad_f; + + // for (UInt q = 0; q < nb_quad_per_facet; ++q) { + // for (UInt s = 0; s < Model::spatial_dimension; ++s) { + // el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + + // f * nb_quad_per_facet + q, + // s) = quad_f(global_facet * nb_quad_per_facet + q, + // s); + // } + // } + //} } } } /// loop over non cohesive materials for (auto && material : materials) { if (dynamic_cast(material.get())) continue; /// initialize the interpolation function material->initElementalFieldInterpolation(elements_quad_facets); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe for (auto & material : this->materials) { try { auto & mat = dynamic_cast(*material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(tangents, tangent_components, // Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); for (auto & material : materials) { - try { - auto & mat __attribute__((unused)) = - dynamic_cast(*material); - } catch (std::bad_cast &) { + auto * mat = dynamic_cast(material.get()); + if (mat == nullptr) /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); - } } #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData( debug::_dm_model_cohesive, "Interpolated stresses before", facet_stress); #endif this->synchronize(_gst_smmc_facets_stress); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_model_cohesive, "Interpolated stresses", facet_stress); #endif } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { + AKANTU_DEBUG_IN(); + interpolateStress(); for (auto & mat : materials) { - try { - auto & mat_cohesive = dynamic_cast(*mat); + auto * mat_cohesive = dynamic_cast(mat.get()); + if(mat_cohesive) { /// check which not ghost cohesive elements are to be created - mat_cohesive.checkInsertion(); - } catch (std::bad_cast &) { + mat_cohesive->checkInsertion(); } } /// communicate data among processors this->synchronize(_gst_smmc_facets); /// insert cohesive elements UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } + AKANTU_DEBUG_OUT(); + return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); updateCohesiveSynchronizers(); SolidMechanicsModel::onElementsAdded(element_list, event); if (cohesive_synchronizer != nullptr) cohesive_synchronizer->computeAllBufferSizes(*this); if (is_extrinsic) resizeFacetStress(); /// if (method != _explicit_lumped_mass) { /// this->initSolver(); /// } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); // Array nodes_list(nb_new_nodes); // for (UInt n = 0; n < nb_new_nodes; ++n) // nodes_list(n) = doubled_nodes(n, 1); SolidMechanicsModel::onNodesAdded(new_nodes, event); UInt new_node, old_node; try { const auto & cohesive_event = dynamic_cast(event); const auto & old_nodes = cohesive_event.getOldNodesList(); auto copy = [this, &new_node, &old_node](auto & arr) { for (UInt s = 0; s < spatial_dimension; ++s) { arr(new_node, s) = arr(old_node, s); } }; for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; copy(*displacement); copy(*blocked_dofs); if (velocity) copy(*velocity); if (acceleration) copy(*acceleration); if (current_position) copy(*current_position); if (previous_displacement) copy(*previous_displacement); } // if (this->getDOFManager().hasMatrix("M")) { // this->assembleMass(old_nodes); // } // if (this->getDOFManager().hasLumpedMatrix("M")) { // this->assembleMassLumped(old_nodes); // } } catch (std::bad_cast &) { } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onEndSolveStep(const AnalysisMethod &) { AKANTU_DEBUG_IN(); /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ for (auto & mat : materials) { if (mat->isFiniteDeformation()) mat->computeAllCauchyStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "SolidMechanicsModelCohesive [" << std::endl; SolidMechanicsModel::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); this->facet_stress.initialize(getFEEngine("FacetsFEEngine"), _nb_component = 2 * spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension - 1); // for (auto && ghost_type : ghost_types) { // for (const auto & type : // mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { // UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); // UInt nb_quadrature_points = getFEEngine("FacetsFEEngine") // .getNbIntegrationPoints(type, // ghost_type); // UInt new_size = nb_facet * nb_quadrature_points; // facet_stress(type, ghost_type).resize(new_size); // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc index c9deec2a9..ba6037031 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc @@ -1,399 +1,484 @@ /** * @file solid_mechanics_model_cohesive_parallel.cc * * @author Marco Vocialta * * * @brief Functions for parallel cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ +#include "communicator.hh" #include "element_synchronizer.hh" +#include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "solid_mechanics_model_tmpl.hh" -#include "communicator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::synchronizeGhostFacetsConnectivity() { AKANTU_DEBUG_IN(); const Communicator & comm = mesh.getCommunicator(); Int psize = comm.getNbProc(); if (psize > 1) { /// get global connectivity for not ghost facets global_connectivity = new ElementTypeMapArray("global_connectivity", id); auto & mesh_facets = inserter->getMeshFacets(); global_connectivity->initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _with_nb_nodes_per_element = true, - _element_kind = _ek_regular, - _ghost_type = _not_ghost); + _element_kind = _ek_regular); mesh_facets.getGlobalConnectivity(*global_connectivity); /// communicate synchronize(_gst_smmc_facets_conn); /// flip facets MeshUtils::flipFacets(mesh_facets, *global_connectivity, _ghost); delete global_connectivity; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateCohesiveSynchronizers() { /// update synchronizers if needed if (not mesh.isDistributed()) return; auto & mesh_facets = inserter->getMeshFacets(); auto & facet_synchronizer = mesh_facets.getElementSynchronizer(); const auto & cfacet_synchronizer = facet_synchronizer; // update the cohesive element synchronizer - cohesive_synchronizer->updateSchemes( - [&](auto && scheme, auto && proc, auto && direction) { - auto & facet_scheme = - cfacet_synchronizer.getCommunications().getScheme(proc, direction); - - for (auto && facet : facet_scheme) { - const auto & connected_elements = mesh_facets.getElementToSubelement( - facet.type, - facet.ghost_type)(facet.element); // slow access here - const auto & cohesive_element = connected_elements[1]; - - auto && cohesive_type = FEEngine::getCohesiveElementType(facet.type); - auto old_nb_cohesive_elements = - mesh.getNbElement(cohesive_type, facet.ghost_type); - old_nb_cohesive_elements -= - mesh.getData("facet_to_double", facet.type, - facet.ghost_type).size(); - - if (cohesive_element.kind() == _ek_cohesive and - cohesive_element.element >= old_nb_cohesive_elements) { - scheme.push_back(cohesive_element); - } - } - }); + cohesive_synchronizer->updateSchemes([&](auto && scheme, auto && proc, + auto && direction) { + auto & facet_scheme = + cfacet_synchronizer.getCommunications().getScheme(proc, direction); + + for (auto && facet : facet_scheme) { + const auto & connected_elements = mesh_facets.getElementToSubelement( + facet.type, + facet.ghost_type)(facet.element); // slow access here + const auto & cohesive_element = connected_elements[1]; + + auto && cohesive_type = FEEngine::getCohesiveElementType(facet.type); + auto old_nb_cohesive_elements = + mesh.getNbElement(cohesive_type, facet.ghost_type); + old_nb_cohesive_elements -= + mesh.getData("facet_to_double", facet.type, facet.ghost_type) + .size(); + + if (cohesive_element.kind() == _ek_cohesive and + cohesive_element.element >= old_nb_cohesive_elements) { + scheme.push_back(cohesive_element); + } + } + }); const auto & comm = mesh.getCommunicator(); auto && my_rank = comm.whoAmI(); // update the facet stress synchronizer - facet_stress_synchronizer->updateSchemes([&](auto && scheme, auto && proc, - auto && /*direction*/) { - auto it_element = scheme.begin(); - for (auto && element : scheme) { - auto && facet_check = inserter->getCheckFacets( - element.type, element.ghost_type)(element.element); // slow access - // here - - if (facet_check) { - auto && connected_elements = mesh_facets.getElementToSubelement( + if (facet_stress_synchronizer) + facet_stress_synchronizer->updateSchemes([&](auto && scheme, auto && proc, + auto && /*direction*/) { + auto it_element = scheme.begin(); + for (auto && element : scheme) { + auto && facet_check = inserter->getCheckFacets( element.type, element.ghost_type)(element.element); // slow access // here - auto && rank_left = facet_synchronizer.getRank(connected_elements[0]); - auto && rank_right = facet_synchronizer.getRank(connected_elements[1]); - - // keep element if the element is still a boundary element between two - // processors - if ((rank_left == Int(proc) and rank_right == my_rank) or - (rank_left == my_rank and rank_right == Int(proc))) { - *it_element = element; - ++it_element; + + if (facet_check) { + auto && connected_elements = mesh_facets.getElementToSubelement( + element.type, element.ghost_type)(element.element); // slow access + // here + auto && rank_left = facet_synchronizer.getRank(connected_elements[0]); + auto && rank_right = + facet_synchronizer.getRank(connected_elements[1]); + + // keep element if the element is still a boundary element between two + // processors + if ((rank_left == Int(proc) and rank_right == my_rank) or + (rank_left == my_rank and rank_right == Int(proc))) { + *it_element = element; + ++it_element; + } } } - } - scheme.resize(it_element - scheme.begin()); - }); + scheme.resize(it_element - scheme.begin()); + }); +} + +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModelCohesive::updateFacetStressSynchronizer() { + if (facet_stress_synchronizer != nullptr) { + const auto & rank_to_element = + mesh.getElementSynchronizer().getElementToRank(); + const auto & facet_checks = inserter->getCheckFacets(); + const auto & mesh_facets = inserter->getMeshFacets(); + const auto & element_to_facet = mesh_facets.getElementToSubelement(); + UInt rank = mesh.getCommunicator().whoAmI(); + + facet_stress_synchronizer->updateSchemes( + [&](auto & scheme, auto & proc, auto & /*direction*/) { + UInt el = 0; + for (auto && element : scheme) { + if (not facet_checks(element)) + continue; + + const auto & next_el = element_to_facet(element); + UInt rank_left = rank_to_element(next_el[0]); + UInt rank_right = rank_to_element(next_el[1]); + + if ((rank_left == rank and rank_right == proc) or + (rank_left == proc and rank_right == rank)) { + scheme[el] = element; + ++el; + } + } + scheme.resize(el); + }); + } } /* -------------------------------------------------------------------------- */ template void SolidMechanicsModelCohesive::packFacetStressDataHelper( const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements) const { packUnpackFacetStressDataHelper( const_cast &>(data_to_pack), buffer, elements); } /* -------------------------------------------------------------------------- */ template void SolidMechanicsModelCohesive::unpackFacetStressDataHelper( ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements) const { packUnpackFacetStressDataHelper(data_to_unpack, buffer, elements); } /* -------------------------------------------------------------------------- */ template void SolidMechanicsModelCohesive::packUnpackFacetStressDataHelper( ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements) const { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt nb_quad_per_elem = 0; UInt sp2 = spatial_dimension * spatial_dimension; UInt nb_component = sp2 * 2; bool element_rank = false; Mesh & mesh_facets = inserter->getMeshFacets(); Array * vect = nullptr; Array> * element_to_facet = nullptr; auto & fe_engine = this->getFEEngine("FacetsFEEngine"); for (auto && el : elements) { if (el.type == _not_defined) - AKANTU_EXCEPTION("packUnpackFacetStressDataHelper called with wrong inputs"); + AKANTU_EXCEPTION( + "packUnpackFacetStressDataHelper called with wrong inputs"); if (el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; vect = &data_to_pack(el.type, el.ghost_type); element_to_facet = &(mesh_facets.getElementToSubelement(el.type, el.ghost_type)); nb_quad_per_elem = fe_engine.getNbIntegrationPoints(el.type, el.ghost_type); } if (pack_helper) element_rank = (*element_to_facet)(el.element)[0].ghost_type != _not_ghost; else element_rank = (*element_to_facet)(el.element)[0].ghost_type == _not_ghost; for (UInt q = 0; q < nb_quad_per_elem; ++q) { Vector data(vect->storage() + (el.element * nb_quad_per_elem + q) * nb_component + element_rank * sp2, sp2); if (pack_helper) buffer << data; else buffer >> data; } } } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::getNbQuadsForFacetCheck( const Array & elements) const { UInt nb_quads = 0; UInt nb_quad_per_facet = 0; ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; auto & fe_engine = this->getFEEngine("FacetsFEEngine"); - for(auto & el : elements) { + for (auto & el : elements) { if (el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; - nb_quad_per_facet = fe_engine - .getNbIntegrationPoints(el.type, el.ghost_type); + nb_quad_per_facet = + fe_engine.getNbIntegrationPoints(el.type, el.ghost_type); } nb_quads += nb_quad_per_facet; } return nb_quads; } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::getNbData( const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; if (elements.size() == 0) return 0; /// regular element case if (elements(0).kind() == _ek_regular) { switch (tag) { case _gst_smmc_facets: { size += elements.size() * sizeof(bool); break; } case _gst_smmc_facets_conn: { UInt nb_nodes = Mesh::getNbNodesPerElementList(elements); size += nb_nodes * sizeof(UInt); break; } case _gst_smmc_facets_stress: { UInt nb_quads = getNbQuadsForFacetCheck(elements); size += nb_quads * spatial_dimension * spatial_dimension * sizeof(Real); break; } + case _gst_material_id: { + for (auto && element : elements) { + if (Mesh::getSpatialDimension(element.type) == (spatial_dimension - 1)) + size += sizeof(UInt); + } + + size += SolidMechanicsModel::getNbData(elements, tag); + break; + } + default: { size += SolidMechanicsModel::getNbData(elements, tag); } } } /// cohesive element case else if (elements(0).kind() == _ek_cohesive) { switch (tag) { case _gst_material_id: { size += elements.size() * sizeof(UInt); break; } case _gst_smm_boundary: { UInt nb_nodes_per_element = 0; for (auto && el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } // force, displacement, boundary size += nb_nodes_per_element * spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } default: break; } if (tag != _gst_material_id && tag != _gst_smmc_facets) { splitByMaterial(elements, [&](auto && mat, auto && elements) { size += mat.getNbData(elements, tag); }); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::packData( CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); if (elements.size() == 0) return; if (elements(0).kind() == _ek_regular) { switch (tag) { case _gst_smmc_facets: { packElementalDataHelper(inserter->getInsertionFacetsByElement(), buffer, elements, false, getFEEngine()); break; } case _gst_smmc_facets_conn: { packElementalDataHelper(*global_connectivity, buffer, elements, false, getFEEngine()); break; } case _gst_smmc_facets_stress: { packFacetStressDataHelper(facet_stress, buffer, elements); break; } + case _gst_material_id: { + for (auto && element : elements) { + if (Mesh::getSpatialDimension(element.type) != (spatial_dimension - 1)) + continue; + buffer << material_index(element); + } + + SolidMechanicsModel::packData(buffer, elements, tag); + break; + } default: { SolidMechanicsModel::packData(buffer, elements, tag); } } - } else if (elements(0).kind() == _ek_cohesive) { + + AKANTU_DEBUG_OUT(); + return; + } + + if (elements(0).kind() == _ek_cohesive) { switch (tag) { case _gst_material_id: { packElementalDataHelper(material_index, buffer, elements, false, getFEEngine("CohesiveFEEngine")); break; } case _gst_smm_boundary: { packNodalDataHelper(*internal_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != _gst_material_id && tag != _gst_smmc_facets) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.packData(buffer, elements, tag); }); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); if (elements.size() == 0) return; if (elements(0).kind() == _ek_regular) { switch (tag) { case _gst_smmc_facets: { unpackElementalDataHelper(inserter->getInsertionFacetsByElement(), buffer, elements, false, getFEEngine()); break; } case _gst_smmc_facets_conn: { unpackElementalDataHelper(*global_connectivity, buffer, elements, false, getFEEngine()); break; } case _gst_smmc_facets_stress: { unpackFacetStressDataHelper(facet_stress, buffer, elements); break; } + case _gst_material_id: { + for (auto && element : elements) { + if (Mesh::getSpatialDimension(element.type) != (spatial_dimension - 1)) + continue; + + UInt recv_mat_index; + buffer >> recv_mat_index; + UInt & mat_index = material_index(element); + if (mat_index != UInt(-1)) + continue; + + // add ghosts element to the correct material + mat_index = recv_mat_index; + auto & mat = dynamic_cast(*materials[mat_index]); + mat.addFacet(element); + facet_material(element) = recv_mat_index; + } + SolidMechanicsModel::unpackData(buffer, elements, tag); + break; + } default: { SolidMechanicsModel::unpackData(buffer, elements, tag); } } - } else if (elements(0).kind() == _ek_cohesive) { + + AKANTU_DEBUG_OUT(); + return; + } + + if (elements(0).kind() == _ek_cohesive) { switch (tag) { case _gst_material_id: { unpackElementalDataHelper(material_index, buffer, elements, false, getFEEngine("CohesiveFEEngine")); break; } case _gst_smm_boundary: { unpackNodalDataHelper(*internal_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != _gst_material_id && tag != _gst_smmc_facets) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.unpackData(buffer, elements, tag); }); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/communications.hh b/src/synchronizer/communications.hh index 528a8274f..58f56ca01 100644 --- a/src/synchronizer/communications.hh +++ b/src/synchronizer/communications.hh @@ -1,272 +1,273 @@ /** * @file communications.hh * * @author Nicolas Richart * * @date Wed Aug 24 13:56:14 2016 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 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 "communication_descriptor.hh" #include "communicator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMUNICATIONS_HH__ #define __AKANTU_COMMUNICATIONS_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ template class Communications { public: using Scheme = Array; protected: using CommunicationPerProcs = std::map; using CommunicationsPerTags = std::map; using CommunicationSchemes = std::map; using Request = std::map>; friend class CommunicationDescriptor; public: using scheme_iterator = typename CommunicationSchemes::iterator; using const_scheme_iterator = typename CommunicationSchemes::const_iterator; /* ------------------------------------------------------------------------ */ class iterator; class tag_iterator; /* ------------------------------------------------------------------------ */ public: CommunicationPerProcs & getCommunications(const SynchronizationTag & tag, const CommunicationSendRecv & sr); /* ------------------------------------------------------------------------ */ bool hasPending(const SynchronizationTag & tag, const CommunicationSendRecv & sr) const; UInt getPending(const SynchronizationTag & tag, const CommunicationSendRecv & sr) const; /* ------------------------------------------------------------------------ */ iterator begin(const SynchronizationTag & tag, const CommunicationSendRecv & sr); iterator end(const SynchronizationTag & tag, const CommunicationSendRecv & sr); /* ------------------------------------------------------------------------ */ iterator waitAny(const SynchronizationTag & tag, const CommunicationSendRecv & sr); /* ------------------------------------------------------------------------ */ void waitAll(const SynchronizationTag & tag, const CommunicationSendRecv & sr); void incrementPending(const SynchronizationTag & tag, const CommunicationSendRecv & sr); void decrementPending(const SynchronizationTag & tag, const CommunicationSendRecv & sr); void freeRequests(const SynchronizationTag & tag, const CommunicationSendRecv & sr); /* ------------------------------------------------------------------------ */ Scheme & createScheme(UInt proc, const CommunicationSendRecv & sr); void resetSchemes(const CommunicationSendRecv & sr); /* ------------------------------------------------------------------------ */ void setCommunicationSize(const SynchronizationTag & tag, UInt proc, UInt size, const CommunicationSendRecv & sr); public: explicit Communications(const Communicator & communicator); + explicit Communications(const Communications & communications); /* ------------------------------------------------------------------------ */ class IterableCommunicationDesc { public: IterableCommunicationDesc(Communications & communications, SynchronizationTag tag, CommunicationSendRecv sr) : communications(communications), tag(std::move(tag)), sr(std::move(sr)) {} auto begin() { return communications.begin(tag, sr); } auto end() { return communications.end(tag, sr); } private: Communications & communications; SynchronizationTag tag; CommunicationSendRecv sr; }; auto iterateRecv(const SynchronizationTag & tag) { return IterableCommunicationDesc(*this, tag, _recv); } auto iterateSend(const SynchronizationTag & tag) { return IterableCommunicationDesc(*this, tag, _send); } /* ------------------------------------------------------------------------ */ // iterator begin_send(const SynchronizationTag & tag); // iterator end_send(const SynchronizationTag & tag); /* ------------------------------------------------------------------------ */ // iterator begin_recv(const SynchronizationTag & tag); // iterator end_recv(const SynchronizationTag & tag); /* ------------------------------------------------------------------------ */ class IterableTags { public: explicit IterableTags(Communications & communications) : communications(communications) {} decltype(auto) begin() { return communications.begin_tag(); } decltype(auto) end() { return communications.end_tag(); } private: Communications & communications; }; decltype(auto) iterateTags() { return IterableTags(*this); } tag_iterator begin_tag(); tag_iterator end_tag(); /* ------------------------------------------------------------------------ */ bool hasCommunication(const SynchronizationTag & tag) const; void incrementCounter(const SynchronizationTag & tag); UInt getCounter(const SynchronizationTag & tag) const; bool hasCommunicationSize(const SynchronizationTag & tag) const; void invalidateSizes(); /* ------------------------------------------------------------------------ */ bool hasPendingRecv(const SynchronizationTag & tag) const; bool hasPendingSend(const SynchronizationTag & tag) const; const auto & getCommunicator() const; /* ------------------------------------------------------------------------ */ iterator waitAnyRecv(const SynchronizationTag & tag); iterator waitAnySend(const SynchronizationTag & tag); void waitAllRecv(const SynchronizationTag & tag); void waitAllSend(const SynchronizationTag & tag); void freeSendRequests(const SynchronizationTag & tag); void freeRecvRequests(const SynchronizationTag & tag); /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ class IterableSchemes { public: IterableSchemes(Communications & communications, CommunicationSendRecv sr) : communications(communications), sr(std::move(sr)) {} decltype(auto) begin() { return communications.begin_scheme(sr); } decltype(auto) end() { return communications.end_scheme(sr); } private: Communications & communications; CommunicationSendRecv sr; }; class ConstIterableSchemes { public: ConstIterableSchemes(const Communications & communications, CommunicationSendRecv sr) : communications(communications), sr(std::move(sr)) {} decltype(auto) begin() const { return communications.begin_scheme(sr); } decltype(auto) end() const { return communications.end_scheme(sr); } private: const Communications & communications; CommunicationSendRecv sr; }; decltype(auto) iterateSchemes(const CommunicationSendRecv & sr) { return IterableSchemes(*this, sr); } decltype(auto) iterateSchemes(const CommunicationSendRecv & sr) const { return ConstIterableSchemes(*this, sr); } decltype(auto) iterateSendSchemes() { return IterableSchemes(*this, _send); } decltype(auto) iterateSendSchemes() const { return ConstIterableSchemes(*this, _send); } decltype(auto) iterateRecvSchemes() { return IterableSchemes(*this, _recv); } decltype(auto) iterateRecvSchemes() const { return ConstIterableSchemes(*this, _recv); } scheme_iterator begin_scheme(const CommunicationSendRecv & sr); scheme_iterator end_scheme(const CommunicationSendRecv & sr); const_scheme_iterator begin_scheme(const CommunicationSendRecv & sr) const; const_scheme_iterator end_scheme(const CommunicationSendRecv & sr) const; /* ------------------------------------------------------------------------ */ scheme_iterator begin_send_scheme(); scheme_iterator end_send_scheme(); const_scheme_iterator begin_send_scheme() const; const_scheme_iterator end_send_scheme() const; /* ------------------------------------------------------------------------ */ scheme_iterator begin_recv_scheme(); scheme_iterator end_recv_scheme(); const_scheme_iterator begin_recv_scheme() const; const_scheme_iterator end_recv_scheme() const; /* ------------------------------------------------------------------------ */ Scheme & createSendScheme(UInt proc); Scheme & createRecvScheme(UInt proc); /* ------------------------------------------------------------------------ */ Scheme & getScheme(UInt proc, const CommunicationSendRecv & sr); const Scheme & getScheme(UInt proc, const CommunicationSendRecv & sr) const; /* ------------------------------------------------------------------------ */ void resetSchemes(); /* ------------------------------------------------------------------------ */ void setSendCommunicationSize(const SynchronizationTag & tag, UInt proc, UInt size); void setRecvCommunicationSize(const SynchronizationTag & tag, UInt proc, UInt size); void initializeCommunications(const SynchronizationTag & tag); protected: CommunicationSchemes schemes[2]; CommunicationsPerTags communications[2]; std::map comm_counter; std::map pending_communications[2]; std::map comm_size_computed; const Communicator & communicator; }; } // namespace akantu #include "communications_tmpl.hh" #endif /* __AKANTU_COMMUNICATIONS_HH__ */ diff --git a/src/synchronizer/communications_tmpl.hh b/src/synchronizer/communications_tmpl.hh index 5bef47632..f87ced61a 100644 --- a/src/synchronizer/communications_tmpl.hh +++ b/src/synchronizer/communications_tmpl.hh @@ -1,534 +1,550 @@ /** * @file communications_tmpl.hh * * @author Nicolas Richart * * @date Tue Sep 6 17:14:06 2016 * * @brief Implementation of Communications * * @section LICENSE * * Copyright (©) 2010-2011 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 "communications.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMUNICATIONS_TMPL_HH__ #define __AKANTU_COMMUNICATIONS_TMPL_HH__ namespace akantu { +/* -------------------------------------------------------------------------- */ +template +Communications::Communications(const Communicator & communicator) + : communicator(communicator) {} + +/* -------------------------------------------------------------------------- */ +template +Communications::Communications(const Communications & other) + : communicator(other.communicator) { + for(auto && sr : iterate_send_recv) { + for (const auto & scheme_pair : other.iterateSchemes(sr)) { + auto proc = scheme_pair.first; + auto & other_scheme = scheme_pair.second; + auto & scheme = this->createScheme(proc, sr); + scheme.copy(other_scheme); + } + } + + this->invalidateSizes(); +} + /* -------------------------------------------------------------------------- */ template class Communications::iterator { using communication_iterator = typename std::map::iterator; public: iterator() : communications(nullptr) {} iterator(scheme_iterator scheme_it, communication_iterator comm_it, Communications & communications, const SynchronizationTag & tag) : scheme_it(scheme_it), comm_it(comm_it), communications(&communications), tag(tag) {} iterator & operator=(const iterator & other) { if (this != &other) { this->scheme_it = other.scheme_it; this->comm_it = other.comm_it; this->communications = other.communications; this->tag = other.tag; } return *this; } iterator & operator++() { ++scheme_it; ++comm_it; return *this; } CommunicationDescriptor operator*() { AKANTU_DEBUG_ASSERT( scheme_it->first == comm_it->first, "The two iterators are not in phase, something wrong" << " happened, time to take out your favorite debugger (" << scheme_it->first << " != " << comm_it->first << ")"); return CommunicationDescriptor(comm_it->second, scheme_it->second, *communications, tag, scheme_it->first); } bool operator==(const iterator & other) const { return scheme_it == other.scheme_it && comm_it == other.comm_it; } bool operator!=(const iterator & other) const { return scheme_it != other.scheme_it || comm_it != other.comm_it; } private: scheme_iterator scheme_it; communication_iterator comm_it; Communications * communications; SynchronizationTag tag; }; /* -------------------------------------------------------------------------- */ template class Communications::tag_iterator { using internal_iterator = std::map::const_iterator; public: tag_iterator(const internal_iterator & it) : it(it) {} tag_iterator & operator++() { ++it; return *this; } SynchronizationTag operator*() { return it->first; } bool operator==(const tag_iterator & other) const { return it == other.it; } bool operator!=(const tag_iterator & other) const { return it != other.it; } private: internal_iterator it; }; /* -------------------------------------------------------------------------- */ template typename Communications::CommunicationPerProcs & Communications::getCommunications(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto comm_it = this->communications[sr].find(tag); if (comm_it == this->communications[sr].end()) AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "No known communications for the tag: " << tag); return comm_it->second; } /* ---------------------------------------------------------------------- */ template UInt Communications::getPending( const SynchronizationTag & tag, const CommunicationSendRecv & sr) const { const std::map & pending = pending_communications[sr]; auto it = pending.find(tag); if (it == pending.end()) return 0; return it->second; } /* -------------------------------------------------------------------------- */ template bool Communications::hasPending( const SynchronizationTag & tag, const CommunicationSendRecv & sr) const { return this->hasCommunication(tag) && (this->getPending(tag, sr) != 0); } /* ---------------------------------------------------------------------- */ template typename Communications::iterator Communications::begin(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto & comms = this->getCommunications(tag, sr); return iterator(this->schemes[sr].begin(), comms.begin(), *this, tag); } template typename Communications::iterator Communications::end(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto & comms = this->getCommunications(tag, sr); return iterator(this->schemes[sr].end(), comms.end(), *this, tag); } /* ---------------------------------------------------------------------- */ template typename Communications::iterator Communications::waitAny(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto & comms = this->getCommunications(tag, sr); auto it = comms.begin(); auto end = comms.end(); std::vector requests; for (; it != end; ++it) { auto & request = it->second.request(); if (!request.isFreed()) requests.push_back(request); } UInt req_id = communicator.waitAny(requests); if (req_id != UInt(-1)) { auto & request = requests[req_id]; UInt proc = sr == _recv ? request.getSource() : request.getDestination(); return iterator(this->schemes[sr].find(proc), comms.find(proc), *this, tag); } else { return this->end(tag, sr); } } /* ---------------------------------------------------------------------- */ template void Communications::waitAll(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto & comms = this->getCommunications(tag, sr); auto it = comms.begin(); auto end = comms.end(); std::vector requests; for (; it != end; ++it) { requests.push_back(it->second.request()); } communicator.waitAll(requests); } template void Communications::incrementPending( const SynchronizationTag & tag, const CommunicationSendRecv & sr) { ++(pending_communications[sr][tag]); } template void Communications::decrementPending( const SynchronizationTag & tag, const CommunicationSendRecv & sr) { --(pending_communications[sr][tag]); } template void Communications::freeRequests(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { iterator it = this->begin(tag, sr); iterator end = this->end(tag, sr); for (; it != end; ++it) { (*it).freeRequest(); } } /* -------------------------------------------------------------------------- */ template typename Communications::Scheme & Communications::createScheme(UInt proc, const CommunicationSendRecv & sr) { // scheme_iterator it = schemes[sr].find(proc); // if (it != schemes[sr].end()) { // AKANTU_CUSTOM_EXCEPTION_INFO(debug::CommunicationException(), // "Communication scheme(" // << sr // << ") already created for proc: " << // proc); // } return schemes[sr][proc]; } template void Communications::resetSchemes(const CommunicationSendRecv & sr) { auto it = this->schemes[sr].begin(); auto end = this->schemes[sr].end(); for (; it != end; ++it) { it->second.resize(0); } } /* -------------------------------------------------------------------------- */ template void Communications::setCommunicationSize( const SynchronizationTag & tag, UInt proc, UInt size, const CommunicationSendRecv & sr) { // accessor that fails if it does not exists comm_size_computed[tag] = true; // TODO: need perhaps to be split based on sr auto & comms = this->communications[sr]; auto & comms_per_tag = comms.at(tag); comms_per_tag.at(proc).resize(size); } /* -------------------------------------------------------------------------- */ template void Communications::initializeCommunications( const SynchronizationTag & tag) { for (auto t = send_recv_t::begin(); t != send_recv_t::end(); ++t) { pending_communications[*t].insert(std::make_pair(tag, 0)); auto & comms = this->communications[*t]; auto & comms_per_tag = comms.insert(std::make_pair(tag, CommunicationPerProcs())) .first->second; for (auto pair : this->schemes[*t]) { comms_per_tag.emplace(std::piecewise_construct, std::forward_as_tuple(pair.first), std::forward_as_tuple(*t)); } } comm_counter.insert(std::make_pair(tag, 0)); } -/* -------------------------------------------------------------------------- */ -template -Communications::Communications(const Communicator & communicator) - : communicator(communicator) {} - /* -------------------------------------------------------------------------- */ template typename Communications::tag_iterator Communications::begin_tag() { return tag_iterator(comm_counter.begin()); } template typename Communications::tag_iterator Communications::end_tag() { return tag_iterator(comm_counter.end()); } /* -------------------------------------------------------------------------- */ template typename Communications::scheme_iterator Communications::begin_scheme(const CommunicationSendRecv & sr) { return this->schemes[sr].begin(); } template typename Communications::scheme_iterator Communications::end_scheme(const CommunicationSendRecv & sr) { return this->schemes[sr].end(); } /* -------------------------------------------------------------------------- */ template typename Communications::const_scheme_iterator Communications::begin_scheme(const CommunicationSendRecv & sr) const { return this->schemes[sr].begin(); } template typename Communications::const_scheme_iterator Communications::end_scheme(const CommunicationSendRecv & sr) const { return this->schemes[sr].end(); } /* -------------------------------------------------------------------------- */ template typename Communications::scheme_iterator Communications::begin_send_scheme() { return this->begin_scheme(_send); } template typename Communications::scheme_iterator Communications::end_send_scheme() { return this->end_scheme(_send); } /* -------------------------------------------------------------------------- */ template typename Communications::const_scheme_iterator Communications::begin_send_scheme() const { return this->begin_scheme(_send); } template typename Communications::const_scheme_iterator Communications::end_send_scheme() const { return this->end_scheme(_send); } /* -------------------------------------------------------------------------- */ template typename Communications::scheme_iterator Communications::begin_recv_scheme() { return this->begin_scheme(_recv); } template typename Communications::scheme_iterator Communications::end_recv_scheme() { return this->end_scheme(_recv); } /* -------------------------------------------------------------------------- */ template typename Communications::const_scheme_iterator Communications::begin_recv_scheme() const { return this->begin_scheme(_recv); } template typename Communications::const_scheme_iterator Communications::end_recv_scheme() const { return this->end_scheme(_recv); } /* ------------------------------------------------------------------------ */ template bool Communications::hasCommunication( const SynchronizationTag & tag) const { return (communications[_send].find(tag) != communications[_send].end()); } template void Communications::incrementCounter(const SynchronizationTag & tag) { auto it = comm_counter.find(tag); if (it == comm_counter.end()) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "No counter initialized in communications for the tags: " << tag); } ++(it->second); } template UInt Communications::getCounter(const SynchronizationTag & tag) const { auto it = comm_counter.find(tag); if (it == comm_counter.end()) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "No counter initialized in communications for the tags: " << tag); } return it->second; } template bool Communications::hasCommunicationSize( const SynchronizationTag & tag) const { auto it = comm_size_computed.find(tag); if (it == comm_size_computed.end()) { return false; } return it->second; } template void Communications::invalidateSizes() { for (auto && pair : comm_size_computed) { pair.second = true; } } template bool Communications::hasPendingRecv( const SynchronizationTag & tag) const { return this->hasPending(tag, _recv); } template bool Communications::hasPendingSend( const SynchronizationTag & tag) const { return this->hasPending(tag, _send); } template const auto & Communications::getCommunicator() const { return communicator; } /* -------------------------------------------------------------------------- */ template typename Communications::iterator Communications::waitAnyRecv(const SynchronizationTag & tag) { return this->waitAny(tag, _recv); } template typename Communications::iterator Communications::waitAnySend(const SynchronizationTag & tag) { return this->waitAny(tag, _send); } template void Communications::waitAllRecv(const SynchronizationTag & tag) { this->waitAll(tag, _recv); } template void Communications::waitAllSend(const SynchronizationTag & tag) { this->waitAll(tag, _send); } template void Communications::freeSendRequests(const SynchronizationTag & tag) { this->freeRequests(tag, _send); } template void Communications::freeRecvRequests(const SynchronizationTag & tag) { this->freeRequests(tag, _recv); } /* -------------------------------------------------------------------------- */ template typename Communications::Scheme & Communications::createSendScheme(UInt proc) { return createScheme(proc, _send); } template typename Communications::Scheme & Communications::createRecvScheme(UInt proc) { return createScheme(proc, _recv); } /* -------------------------------------------------------------------------- */ template void Communications::resetSchemes() { resetSchemes(_send); resetSchemes(_recv); } /* -------------------------------------------------------------------------- */ template typename Communications::Scheme & Communications::getScheme(UInt proc, const CommunicationSendRecv & sr) { return this->schemes[sr].find(proc)->second; } /* -------------------------------------------------------------------------- */ template const typename Communications::Scheme & Communications::getScheme(UInt proc, const CommunicationSendRecv & sr) const { return this->schemes[sr].find(proc)->second; } /* -------------------------------------------------------------------------- */ template void Communications::setSendCommunicationSize( const SynchronizationTag & tag, UInt proc, UInt size) { this->setCommunicationSize(tag, proc, size, _send); } template void Communications::setRecvCommunicationSize( const SynchronizationTag & tag, UInt proc, UInt size) { this->setCommunicationSize(tag, proc, size, _recv); } } // namespace akantu #endif /* __AKANTU_COMMUNICATIONS_TMPL_HH__ */ diff --git a/src/synchronizer/data_accessor.cc b/src/synchronizer/data_accessor.cc index c4855e9ab..3a51a7016 100644 --- a/src/synchronizer/data_accessor.cc +++ b/src/synchronizer/data_accessor.cc @@ -1,163 +1,155 @@ /** * @file data_accessor.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Sun Oct 19 2014 * * @brief data accessors constructor functions * * @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 "data_accessor.hh" #include "fe_engine.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template void DataAccessor::packUnpackNodalDataHelper( Array & data, CommunicationBuffer & buffer, const Array & elements, const Mesh & mesh) { UInt nb_component = data.getNbComponent(); UInt nb_nodes_per_element = 0; ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt * conn = nullptr; - Array::const_iterator it = elements.begin(); - Array::const_iterator end = elements.end(); - for (; it != end; ++it) { - const Element & el = *it; + for (auto & el : elements) { if (el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; conn = mesh.getConnectivity(el.type, el.ghost_type).storage(); nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); } UInt el_offset = el.element * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; Vector data_vect(data.storage() + offset_conn * nb_component, nb_component); if (pack_helper) buffer << data_vect; else buffer >> data_vect; } } } /* ------------------------------------------------------------------------ */ template void DataAccessor::packUnpackElementalDataHelper( ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & element, bool per_quadrature_point_data, const FEEngine & fem) { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt nb_quad_per_elem = 0; UInt nb_component = 0; Array * vect = nullptr; - Array::const_iterator it = element.begin(); - Array::const_iterator end = element.end(); - for (; it != end; ++it) { - const Element & el = *it; + for (auto & el : element) { if (el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; vect = &data_to_pack(el.type, el.ghost_type); - if (per_quadrature_point_data) - nb_quad_per_elem = fem.getNbIntegrationPoints(el.type, el.ghost_type); - else - nb_quad_per_elem = 1; + + nb_quad_per_elem = + per_quadrature_point_data + ? fem.getNbIntegrationPoints(el.type, el.ghost_type) + : 1; nb_component = vect->getNbComponent(); } Vector data(vect->storage() + el.element * nb_component * nb_quad_per_elem, nb_component * nb_quad_per_elem); if (pack_helper) buffer << data; else buffer >> data; } } /* -------------------------------------------------------------------------- */ template void DataAccessor::packUnpackDOFDataHelper(Array & data, CommunicationBuffer & buffer, const Array & dofs) { - Array::const_scalar_iterator it_dof = dofs.begin(); - Array::const_scalar_iterator end_dof = dofs.end(); T * data_ptr = data.storage(); - - for (; it_dof != end_dof; ++it_dof) { + for (const auto & dof : dofs) { if (pack_helper) - buffer << data_ptr[*it_dof]; + buffer << data_ptr[dof]; else - buffer >> data_ptr[*it_dof]; + buffer >> data_ptr[dof]; } } /* -------------------------------------------------------------------------- */ #define DECLARE_HELPERS(T) \ template void DataAccessor::packUnpackNodalDataHelper( \ Array & data, CommunicationBuffer & buffer, \ const Array & elements, const Mesh & mesh); \ template void DataAccessor::packUnpackNodalDataHelper( \ Array & data, CommunicationBuffer & buffer, \ const Array & elements, const Mesh & mesh); \ template void \ DataAccessor::packUnpackElementalDataHelper( \ ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, \ const Array & element, bool per_quadrature_point_data, \ const FEEngine & fem); \ template void DataAccessor::packUnpackElementalDataHelper( \ ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, \ const Array & element, bool per_quadrature_point_data, \ const FEEngine & fem); \ template void DataAccessor::packUnpackDOFDataHelper( \ Array & data, CommunicationBuffer & buffer, \ const Array & dofs); \ template void DataAccessor::packUnpackDOFDataHelper( \ Array & data, CommunicationBuffer & buffer, const Array & dofs) /* -------------------------------------------------------------------------- */ DECLARE_HELPERS(Real); DECLARE_HELPERS(UInt); DECLARE_HELPERS(bool); /* -------------------------------------------------------------------------- */ } // akantu diff --git a/src/synchronizer/element_synchronizer.cc b/src/synchronizer/element_synchronizer.cc index 9ffc5f5f8..32698e302 100644 --- a/src/synchronizer/element_synchronizer.cc +++ b/src/synchronizer/element_synchronizer.cc @@ -1,353 +1,369 @@ /** * @file element_synchronizer.cc * * @author Guillaume Anciaux * @author Dana Christen * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Sep 01 2010 * @date last modification: Fri Jan 22 2016 * * @brief implementation of a communicator using a static_communicator for * real * send/receive * * @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 "element_synchronizer.hh" #include "aka_common.hh" #include "mesh.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ElementSynchronizer::ElementSynchronizer(Mesh & mesh, const ID & id, MemoryID memory_id, bool register_to_event_manager, EventHandlerPriority event_priority) : SynchronizerImpl(mesh.getCommunicator(), id, memory_id), mesh(mesh), element_to_prank("element_to_prank", id, memory_id) { + AKANTU_DEBUG_IN(); + + if (register_to_event_manager) + this->mesh.registerEventHandler(*this, event_priority); + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +ElementSynchronizer::ElementSynchronizer(const ElementSynchronizer & other, + const ID & id, + bool register_to_event_manager, + EventHandlerPriority event_priority) + : SynchronizerImpl(other, id), mesh(other.mesh), + element_to_prank("element_to_prank", id, other.memory_id) { AKANTU_DEBUG_IN(); + element_to_prank.copy(other.element_to_prank); + if (register_to_event_manager) this->mesh.registerEventHandler(*this, event_priority); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementSynchronizer::~ElementSynchronizer() = default; /* -------------------------------------------------------------------------- */ void ElementSynchronizer::substituteElements( const std::map & old_to_new_elements) { auto found_element_end = old_to_new_elements.end(); // substitute old elements with new ones for (auto && sr : iterate_send_recv) { for (auto && scheme_pair : communications.iterateSchemes(sr)) { auto & list = scheme_pair.second; for (auto & el : list) { auto found_element_it = old_to_new_elements.find(el); if (found_element_it != found_element_end) el = found_element_it->second; } } } } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::onElementsChanged( const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray &, const ChangedElementsEvent &) { // create a map to link old elements to new ones std::map old_to_new_elements; for (UInt el = 0; el < old_elements_list.size(); ++el) { AKANTU_DEBUG_ASSERT(old_to_new_elements.find(old_elements_list(el)) == old_to_new_elements.end(), "The same element cannot appear twice in the list"); old_to_new_elements[old_elements_list(el)] = new_elements_list(el); } substituteElements(old_to_new_elements); communications.invalidateSizes(); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::onElementsRemoved( const Array & element_to_remove, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent &) { AKANTU_DEBUG_IN(); this->removeElements(element_to_remove); this->renumberElements(new_numbering); communications.invalidateSizes(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::buildElementToPrank() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); element_to_prank.initialize(mesh, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined, _with_nb_element = true, _default_value = rank); /// assign prank to all ghost elements for (auto && scheme : communications.iterateSchemes(_recv)) { auto & recv = scheme.second; auto proc = scheme.first; for (auto & element : recv) { element_to_prank(element) = proc; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Int ElementSynchronizer::getRank(const Element & element) const { if (not element_to_prank.exists(element.type, element.ghost_type)) { // Nicolas: Ok This is nasty I know.... const_cast(this)->buildElementToPrank(); } return element_to_prank(element); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::reset() { AKANTU_DEBUG_IN(); communications.resetSchemes(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::removeElements( const Array & element_to_remove) { AKANTU_DEBUG_IN(); std::vector send_requests; std::map> list_of_elements_per_proc; auto filter_list = [](const Array & keep, Array & list) { Array new_list; for (UInt e = 0; e < keep.size() - 1; ++e) { Element & el = list(keep(e)); new_list.push_back(el); } list.copy(new_list); }; // Handling ghost elements for (auto && scheme_pair : communications.iterateSchemes(_recv)) { auto proc = scheme_pair.first; auto & recv = scheme_pair.second; Array & keep_list = list_of_elements_per_proc[proc]; auto rem_it = element_to_remove.begin(); auto rem_end = element_to_remove.end(); for (auto && pair : enumerate(recv)) { const auto & el = std::get<1>(pair); auto pos = std::find(rem_it, rem_end, el); if (pos == rem_end) { keep_list.push_back(std::get<0>(pair)); } } keep_list.push_back(UInt(-1)); // To no send empty arrays auto && tag = Tag::genTag(proc, 0, Tag::_MODIFY_SCHEME, this->hash_id); AKANTU_DEBUG_INFO("Sending a message of size " << keep_list.size() << " to proc " << proc << " TAG(" << tag << ")"); send_requests.push_back(this->communicator.asyncSend(keep_list, proc, tag)); #ifndef AKANTU_NDEBUG auto old_size = recv.size(); #endif filter_list(keep_list, recv); AKANTU_DEBUG_INFO("I had " << old_size << " elements to recv from proc " << proc << " and " << keep_list.size() << " elements to keep. I have " << recv.size() << " elements left."); } for (auto && scheme_pair : communications.iterateSchemes(_send)) { auto proc = scheme_pair.first; auto & send = scheme_pair.second; CommunicationStatus status; auto && tag = Tag::genTag(rank, 0, Tag::_MODIFY_SCHEME, this->hash_id); AKANTU_DEBUG_INFO("Getting number of elements of proc " << proc << " to keep - TAG(" << tag << ")"); this->communicator.probe(proc, tag, status); Array keep_list(status.size()); AKANTU_DEBUG_INFO("Receiving list of elements (" << keep_list.size() << " elements) to keep for proc " << proc << " TAG(" << tag << ")"); this->communicator.receive(keep_list, proc, tag); #ifndef AKANTU_NDEBUG auto old_size = send.size(); #endif filter_list(keep_list, send); AKANTU_DEBUG_INFO("I had " << old_size << " elements to send to proc " << proc << " and " << keep_list.size() << " elements to keep. I have " << send.size() << " elements left."); } this->communicator.waitAll(send_requests); this->communicator.freeCommunicationRequest(send_requests); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::renumberElements( const ElementTypeMapArray & new_numbering) { for (auto && sr : iterate_send_recv) { for (auto && scheme_pair : communications.iterateSchemes(sr)) { auto & list = scheme_pair.second; for (auto && el : list) { if (new_numbering.exists(el.type, el.ghost_type)) el.element = new_numbering(el); } } } } /* -------------------------------------------------------------------------- */ UInt ElementSynchronizer::sanityCheckDataSize( const Array & elements, const SynchronizationTag &) const { UInt size = 0; size += sizeof(SynchronizationTag); // tag size += sizeof(UInt); // comm_desc.getNbData(); size += sizeof(UInt); // comm_desc.getProc(); size += sizeof(UInt); // mesh.getCommunicator().whoAmI(); // barycenters size += (elements.size() * mesh.getSpatialDimension() * sizeof(Real)); return size; } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::packSanityCheckData( CommunicationDescriptor & comm_desc) const { auto & buffer = comm_desc.getBuffer(); buffer << comm_desc.getTag(); buffer << comm_desc.getNbData(); buffer << comm_desc.getProc(); buffer << mesh.getCommunicator().whoAmI(); auto & send_element = comm_desc.getScheme(); /// pack barycenters in debug mode for (auto && element : send_element) { Vector barycenter(mesh.getSpatialDimension()); mesh.getBarycenter(element, barycenter); buffer << barycenter; } } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::unpackSanityCheckData( CommunicationDescriptor & comm_desc) const { auto & buffer = comm_desc.getBuffer(); const auto & tag = comm_desc.getTag(); auto nb_data = comm_desc.getNbData(); auto proc = comm_desc.getProc(); auto rank = mesh.getCommunicator().whoAmI(); decltype(nb_data) recv_nb_data; decltype(proc) recv_proc; decltype(rank) recv_rank; SynchronizationTag t; buffer >> t; buffer >> recv_nb_data; buffer >> recv_proc; buffer >> recv_rank; AKANTU_DEBUG_ASSERT( t == tag, "The tag received does not correspond to the tag expected"); AKANTU_DEBUG_ASSERT( nb_data == recv_nb_data, "The nb_data received does not correspond to the nb_data expected"); AKANTU_DEBUG_ASSERT(UInt(recv_rank) == proc, "The rank received does not correspond to the proc"); AKANTU_DEBUG_ASSERT(recv_proc == UInt(rank), "The proc received does not correspond to the rank"); auto & recv_element = comm_desc.getScheme(); auto spatial_dimension = mesh.getSpatialDimension(); for (auto && element : recv_element) { Vector barycenter_loc(spatial_dimension); mesh.getBarycenter(element, barycenter_loc); Vector barycenter(spatial_dimension); buffer >> barycenter; auto dist = barycenter_loc.distance(barycenter); if (not Math::are_float_equal(dist, 0.)) AKANTU_EXCEPTION("Unpacking an unknown value for the element " << element << "(barycenter " << barycenter_loc << " != buffer " << barycenter << ") [" << dist << "] - tag: " << tag << " comm from " << proc << " to " << rank); } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/element_synchronizer.hh b/src/synchronizer/element_synchronizer.hh index 53fbc6e76..8764beb79 100644 --- a/src/synchronizer/element_synchronizer.hh +++ b/src/synchronizer/element_synchronizer.hh @@ -1,211 +1,216 @@ /** * @file element_synchronizer.hh * * @author Guillaume Anciaux * @author Dana Christen * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Dec 08 2015 * * @brief Main element synchronizer * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_SYNCHRONIZER_HH__ #define __AKANTU_ELEMENT_SYNCHRONIZER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_common.hh" #include "mesh_partition.hh" #include "synchronizer_impl.hh" namespace akantu { class Mesh; } /* -------------------------------------------------------------------------- */ namespace akantu { class ElementSynchronizer : public SynchronizerImpl, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ElementSynchronizer(Mesh & mesh, const ID & id = "element_synchronizer", MemoryID memory_id = 0, bool register_to_event_manager = true, EventHandlerPriority event_priority = _ehp_synchronizer); + ElementSynchronizer(const ElementSynchronizer & other, + const ID & id = "element_synchronizer_copy", + bool register_to_event_manager = true, + EventHandlerPriority event_priority = _ehp_synchronizer); + public: ~ElementSynchronizer() override; friend class ElementInfoPerProc; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /// mesh event handler onElementsChanged void onElementsChanged(const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray & new_numbering, const ChangedElementsEvent & event) override; /// mesh event handler onRemovedElement void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; /// mesh event handler onNodesAdded void onNodesAdded(const Array & /* nodes_list*/, const NewNodesEvent & /*event*/) override{}; /// mesh event handler onRemovedNodes void onNodesRemoved(const Array & /*nodes_list*/, const Array & /*new_numbering*/, const RemovedNodesEvent & /*event*/) override{}; /// mesh event handler onElementsAdded void onElementsAdded(const Array & /*elements_list*/, const NewElementsEvent & /*event*/) override{}; protected: /// reset send and recv element lists void reset(); /// remove elements from the synchronizer without renumbering them void removeElements(const Array & element_to_remove); /// renumber the elements in the synchronizer void renumberElements(const ElementTypeMapArray & new_numbering); /// build processor to element correspondence void buildElementToPrank(); protected: /// fill the nodes type vector void fillNodesType(const MeshData & mesh_data, DynamicCommunicationBuffer * buffers, const std::string & tag_name, const ElementType & el_type, const Array & partition_num); template void fillTagBufferTemplated(const MeshData & mesh_data, DynamicCommunicationBuffer * buffers, const std::string & tag_name, const ElementType & el_type, const Array & partition_num, const CSR & ghost_partition); void fillTagBuffer(const MeshData & mesh_data, DynamicCommunicationBuffer * buffers, const std::string & tag_name, const ElementType & el_type, const Array & partition_num, const CSR & ghost_partition); /// function that handels the MeshData to be split (root side) static void synchronizeTagsSend(ElementSynchronizer & communicator, UInt root, Mesh & mesh, UInt nb_tags, const ElementType & type, const Array & partition_num, const CSR & ghost_partition, UInt nb_local_element, UInt nb_ghost_element); /// function that handles the MeshData to be split (other nodes) static void synchronizeTagsRecv(ElementSynchronizer & communicator, UInt root, Mesh & mesh, UInt nb_tags, const ElementType & type, UInt nb_local_element, UInt nb_ghost_element); /// function that handles the preexisting groups in the mesh static void synchronizeElementGroups(ElementSynchronizer & communicator, UInt root, Mesh & mesh, const ElementType & type, const Array & partition_num, const CSR & ghost_partition, UInt nb_element); /// function that handles the preexisting groups in the mesh static void synchronizeElementGroups(ElementSynchronizer & communicator, UInt root, Mesh & mesh, const ElementType & type); /// function that handles the preexisting groups in the mesh static void synchronizeNodeGroupsMaster(ElementSynchronizer & communicator, UInt root, Mesh & mesh); /// function that handles the preexisting groups in the mesh static void synchronizeNodeGroupsSlaves(ElementSynchronizer & communicator, UInt root, Mesh & mesh); template static void fillNodeGroupsFromBuffer(ElementSynchronizer & communicator, Mesh & mesh, CommunicationBuffer & buffer); /// substitute elements in the send and recv arrays void substituteElements(const std::map & old_to_new_elements); /* ------------------------------------------------------------------------ */ /* Sanity checks */ /* ------------------------------------------------------------------------ */ UInt sanityCheckDataSize(const Array & elements, const SynchronizationTag & tag) const override; /* ------------------------------------------------------------------------ */ void packSanityCheckData( CommunicationDescriptor & comm_desc) const override; /* ------------------------------------------------------------------------ */ void unpackSanityCheckData( CommunicationDescriptor & comm_desc) const override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Mesh, mesh, Mesh &); AKANTU_GET_MACRO(ElementToRank, element_to_prank, const ElementTypeMapArray &); Int getRank(const Element & element) const final; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// reference to the underlying mesh Mesh & mesh; friend class FacetSynchronizer; ElementTypeMapArray element_to_prank; }; /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_ELEMENT_SYNCHRONIZER_HH__ */ diff --git a/src/synchronizer/synchronizer.hh b/src/synchronizer/synchronizer.hh index 0f809526e..8f3501158 100644 --- a/src/synchronizer/synchronizer.hh +++ b/src/synchronizer/synchronizer.hh @@ -1,123 +1,125 @@ /** * @file synchronizer.hh * * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Dec 10 2015 * * @brief Common interface for synchronizers * * @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_memory.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SYNCHRONIZER_HH__ #define __AKANTU_SYNCHRONIZER_HH__ namespace akantu { class Communicator; } namespace akantu { /* -------------------------------------------------------------------------- */ /* Base class for synchronizers */ /* -------------------------------------------------------------------------- */ class Synchronizer : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Synchronizer(const Communicator & comm, const ID & id = "synchronizer", MemoryID memory_id = 0); + Synchronizer(const Synchronizer & other) = default; + ~Synchronizer() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// synchronize ghosts without state template void synchronizeOnce(DataAccessor & data_accessor, const SynchronizationTag & tag) const; /// synchronize ghosts template void synchronize(DataAccessor & data_accessor, const SynchronizationTag & tag); /// asynchronous synchronization of ghosts template void asynchronousSynchronize(const DataAccessor & data_accessor, const SynchronizationTag & tag); /// wait end of asynchronous synchronization of ghosts template void waitEndSynchronize(DataAccessor & data_accessor, const SynchronizationTag & tag); /// compute buffer size for a given tag and data accessor template void computeBufferSize(const DataAccessor & data_accessor, const SynchronizationTag & tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Communicator, communicator, const Communicator &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the synchronizer ID id; /// hashed version of the id int hash_id; /// message counter per tag std::map tag_counter; /// the static memory instance const Communicator & communicator; /// nb processors in the communicator UInt nb_proc; /// rank in the communicator UInt rank; }; } // namespace akantu #include "synchronizer_tmpl.hh" #endif /* __AKANTU_SYNCHRONIZER_HH__ */ diff --git a/src/synchronizer/synchronizer_impl.hh b/src/synchronizer/synchronizer_impl.hh index b4b87c96d..fcfedb3e8 100644 --- a/src/synchronizer/synchronizer_impl.hh +++ b/src/synchronizer/synchronizer_impl.hh @@ -1,121 +1,123 @@ /** * @file synchronizer_impl.hh * * @author Nicolas Richart * * @date Wed Aug 24 13:26:47 2016 * * @brief Implementation of the generic part of synchronizers * * @section LICENSE * * Copyright (©) 2010-2011 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 "communications.hh" #include "synchronizer.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SYNCHRONIZER_IMPL_HH__ #define __AKANTU_SYNCHRONIZER_IMPL_HH__ namespace akantu { template class SynchronizerImpl : public Synchronizer { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SynchronizerImpl(const Communicator & communicator, const ID & id = "synchronizer", MemoryID memory_id = 0); + SynchronizerImpl(const SynchronizerImpl & other, const ID & id); + ~SynchronizerImpl() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// synchronous synchronization without state virtual void synchronizeOnceImpl(DataAccessor & data_accessor, const SynchronizationTag & tag) const; /// asynchronous synchronization of ghosts virtual void asynchronousSynchronizeImpl(const DataAccessor & data_accessor, const SynchronizationTag & tag); /// wait end of asynchronous synchronization of ghosts virtual void waitEndSynchronizeImpl(DataAccessor & data_accessor, const SynchronizationTag & tag); /// compute all buffer sizes virtual void computeAllBufferSizes(const DataAccessor & data_accessor); /// compute buffer size for a given tag and data accessor virtual void computeBufferSizeImpl(const DataAccessor & data_accessor, const SynchronizationTag & tag); /* ------------------------------------------------------------------------ */ virtual void synchronizeImpl(DataAccessor & data_accessor, const SynchronizationTag & tag) { this->asynchronousSynchronizeImpl(data_accessor, tag); this->waitEndSynchronizeImpl(data_accessor, tag); } /* ------------------------------------------------------------------------ */ /// extract the elements that have a true predicate from in_synchronizer and /// store them in the current synchronizer template void split(SynchronizerImpl & in_synchronizer, Pred && pred); /// update schemes in a synchronizer template void updateSchemes(Updater && scheme_updater); /* ------------------------------------------------------------------------ */ virtual UInt sanityCheckDataSize(const Array & elements, const SynchronizationTag & tag) const; virtual void packSanityCheckData(CommunicationDescriptor & comm_desc) const; virtual void unpackSanityCheckData(CommunicationDescriptor & comm_desc) const; public: AKANTU_GET_MACRO(Communications, communications, const Communications &); protected: AKANTU_GET_MACRO_NOT_CONST(Communications, communications, Communications &); virtual Int getRank(const Entity & entity) const = 0; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// information on the communications Communications communications; }; } // namespace akantu #include "synchronizer_impl_tmpl.hh" #endif /* __AKANTU_SYNCHRONIZER_IMPL_HH__ */ diff --git a/src/synchronizer/synchronizer_impl_tmpl.hh b/src/synchronizer/synchronizer_impl_tmpl.hh index 8cd7a2cdf..b8129ce99 100644 --- a/src/synchronizer/synchronizer_impl_tmpl.hh +++ b/src/synchronizer/synchronizer_impl_tmpl.hh @@ -1,305 +1,314 @@ /** * @file synchronizer_impl_tmpl.hh * * @author Nicolas Richart * * @date Wed Aug 24 13:29:47 2016 * * @brief Implementation of the SynchronizerImpl * * @section LICENSE * * Copyright (©) 2010-2011 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 "synchronizer_impl.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__ #define __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ template -SynchronizerImpl::SynchronizerImpl(const Communicator & comm, const ID & id, MemoryID memory_id) +SynchronizerImpl::SynchronizerImpl(const Communicator & comm, + const ID & id, MemoryID memory_id) : Synchronizer(comm, id, memory_id), communications(comm) {} +/* -------------------------------------------------------------------------- */ +template +SynchronizerImpl::SynchronizerImpl(const SynchronizerImpl & other, + const ID & id) + : Synchronizer(other), communications(other.communications) { + this->id = id; +} + /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::synchronizeOnceImpl( DataAccessor & data_accessor, const SynchronizationTag & tag) const { // no need to synchronize if (this->nb_proc == 1) return; using CommunicationRequests = std::vector; using CommunicationBuffers = std::map; CommunicationRequests send_requests, recv_requests; CommunicationBuffers send_buffers, recv_buffers; auto postComm = [&](const CommunicationSendRecv & sr, CommunicationBuffers & buffers, CommunicationRequests & requests) -> void { for (auto && pair : communications.iterateSchemes(sr)) { auto & proc = pair.first; const auto & scheme = pair.second; auto & buffer = buffers[proc]; UInt buffer_size = data_accessor.getNbData(scheme, tag); buffer.resize(buffer_size); if (sr == _recv) { requests.push_back(communicator.asyncReceive( buffer, proc, Tag::genTag(this->rank, 0, Tag::_SYNCHRONIZE, this->hash_id))); } else { data_accessor.packData(buffer, scheme, tag); send_requests.push_back(communicator.asyncSend( buffer, proc, Tag::genTag(proc, 0, Tag::_SYNCHRONIZE, this->hash_id))); } } }; // post the receive requests postComm(_recv, recv_buffers, recv_requests); // post the send data requests postComm(_send, send_buffers, send_requests); // treat the receive requests UInt request_ready; while ((request_ready = communicator.waitAny(recv_requests)) != UInt(-1)) { CommunicationRequest & req = recv_requests[request_ready]; UInt proc = req.getSource(); CommunicationBuffer & buffer = recv_buffers[proc]; const auto & scheme = this->communications.getScheme(proc, _recv); data_accessor.unpackData(buffer, scheme, tag); req.free(); recv_requests.erase(recv_requests.begin() + request_ready); } communicator.waitAll(send_requests); communicator.freeCommunicationRequest(send_requests); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::asynchronousSynchronizeImpl( const DataAccessor & data_accessor, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); if (not this->communications.hasCommunicationSize(tag)) this->computeBufferSize(data_accessor, tag); this->communications.incrementCounter(tag); // Posting the receive ------------------------------------------------------- if (this->communications.hasPendingRecv(tag)) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "There must still be some pending receive communications." << " Tag is " << tag << " Cannot start new ones"); } for (auto && comm_desc : this->communications.iterateRecv(tag)) { comm_desc.postRecv(this->hash_id); } // Posting the sends ------------------------------------------------------- if (communications.hasPendingSend(tag)) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "There must be some pending sending communications." << " Tag is " << tag); } for (auto && comm_desc : this->communications.iterateSend(tag)) { comm_desc.resetBuffer(); #ifndef AKANTU_NDEBUG this->packSanityCheckData(comm_desc); #endif comm_desc.packData(data_accessor); comm_desc.postSend(this->hash_id); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::waitEndSynchronizeImpl( DataAccessor & data_accessor, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); #ifndef AKANTU_NDEBUG if (this->communications.begin(tag, _recv) != - this->communications.end(tag, _recv) && + this->communications.end(tag, _recv) && !this->communications.hasPendingRecv(tag)) AKANTU_CUSTOM_EXCEPTION_INFO(debug::CommunicationException(), "No pending communication with the tag \"" << tag); #endif auto recv_end = this->communications.end(tag, _recv); decltype(recv_end) recv_it; while ((recv_it = this->communications.waitAnyRecv(tag)) != recv_end) { auto && comm_desc = *recv_it; #ifndef AKANTU_NDEBUG this->unpackSanityCheckData(comm_desc); #endif comm_desc.unpackData(data_accessor); comm_desc.resetBuffer(); comm_desc.freeRequest(); } this->communications.waitAllSend(tag); this->communications.freeSendRequests(tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::computeAllBufferSizes( const DataAccessor & data_accessor) { for (auto && tag : this->communications.iterateTags()) { this->computeBufferSize(data_accessor, tag); } } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::computeBufferSizeImpl( const DataAccessor & data_accessor, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); if (not this->communications.hasCommunication(tag)) { this->communications.initializeCommunications(tag); AKANTU_DEBUG_ASSERT(communications.hasCommunication(tag) == true, "Communications where not properly initialized"); } for (auto sr : iterate_send_recv) { for (auto && pair : this->communications.iterateSchemes(sr)) { auto proc = pair.first; const auto & scheme = pair.second; UInt size = 0; #ifndef AKANTU_NDEBUG size += this->sanityCheckDataSize(scheme, tag); #endif size += data_accessor.getNbData(scheme, tag); AKANTU_DEBUG_INFO("I have " << size << "(" << printMemorySize(size) << " - " << scheme.size() << " element(s)) data to " << std::string(sr == _recv ? "receive from" : "send to") << proc << " for tag " << tag); this->communications.setCommunicationSize(tag, proc, size, sr); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void SynchronizerImpl::split(SynchronizerImpl & in_synchronizer, Pred && pred) { AKANTU_DEBUG_IN(); auto filter_list = [&](auto & list, auto & new_list) { auto copy = list; list.resize(0); new_list.resize(0); for (auto && entity : copy) { if (std::forward(pred)(entity)) { new_list.push_back(entity); } else { list.push_back(entity); } } }; for (auto sr : iterate_send_recv) { for (auto & scheme_pair : in_synchronizer.communications.iterateSchemes(sr)) { auto proc = scheme_pair.first; auto & scheme = scheme_pair.second; auto & new_scheme = communications.createScheme(proc, sr); - filter_list(new_scheme, scheme); + filter_list(scheme, new_scheme); } } in_synchronizer.communications.invalidateSizes(); communications.invalidateSizes(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void SynchronizerImpl::updateSchemes(Updater && scheme_updater) { for (auto sr : iterate_send_recv) { for (auto & scheme_pair : communications.iterateSchemes(sr)) { auto proc = scheme_pair.first; auto & scheme = scheme_pair.second; std::forward(scheme_updater)(scheme, proc, sr); } } communications.invalidateSizes(); } /* -------------------------------------------------------------------------- */ template UInt SynchronizerImpl::sanityCheckDataSize( const Array &, const SynchronizationTag &) const { return 0; } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::packSanityCheckData( CommunicationDescriptor &) const {} /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::unpackSanityCheckData( CommunicationDescriptor &) const {} /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh index 427f9444a..9d97a7293 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh @@ -1,258 +1,274 @@ /* -------------------------------------------------------------------------- */ +#include "communicator.hh" #include "solid_mechanics_model_cohesive.hh" #include "test_gtest_utils.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TEST_COHESIVE_FIXTURE_HH__ #define __AKANTU_TEST_COHESIVE_FIXTURE_HH__ using namespace akantu; template <::akantu::AnalysisMethod t> using analysis_method_t = std::integral_constant<::akantu::AnalysisMethod, t>; template class TestSMMCFixture : public ::testing::Test { public: static constexpr ElementType cohesive_type = std::tuple_element_t<0, param_>::value; static constexpr ElementType type_1 = std::tuple_element_t<1, param_>::value; static constexpr ElementType type_2 = std::tuple_element_t<2, param_>::value; static constexpr size_t dim = ElementClass::getSpatialDimension(); void SetUp() override { normal = Vector(this->dim, 0.); - normal(this->dim - 1) = 1.; + if (dim == 1) + normal(_x) = 1.; + else + normal(_y) = 1.; mesh = std::make_unique(this->dim); - EXPECT_NO_THROW({ mesh->read(this->mesh_name); }); + if (Communicator::getStaticCommunicator().whoAmI() == 0) { + EXPECT_NO_THROW({ mesh->read(this->mesh_name); }); + } + mesh->distribute(); } void TearDown() override { model.reset(nullptr); mesh.reset(nullptr); } void createModel() { model = std::make_unique(*mesh); model->initFull(_analysis_method = this->analysis_method, _is_extrinsic = this->is_extrinsic); // auto stable_time_step = this->model->getStableTimeStep(); this->model->setTimeStep(4e-6); // std::cout << stable_time_step *0.0 << std::endl; if (dim == 1) { surface = 1; return; } auto facet_type = mesh->getFacetType(this->cohesive_type); auto & fe_engine = model->getFEEngineBoundary(); auto & group = mesh->getElementGroup("insertion").getElements(facet_type); Array ones(fe_engine.getNbIntegrationPoints(facet_type) * group.size()); ones.set(1.); surface = fe_engine.integrate(ones, facet_type, _not_ghost, group); + mesh->getCommunicator().allReduce(surface, SynchronizerOperation::_sum); } void setInitialCondition(const Vector & direction, Real strain) { auto lower = this->mesh->getLowerBounds().dot(normal); auto upper = this->mesh->getUpperBounds().dot(normal); Real L = upper - lower; for (auto && data : zip(make_view(this->mesh->getNodes(), this->dim), make_view(this->model->getDisplacement(), this->dim))) { const auto & pos = std::get<0>(data); auto & disp = std::get<1>(data); auto x = pos.dot(normal) - (upper + lower) / 2.; disp = direction * (x * 2. * strain / L); } } #define debug 1 void steps(const Vector & displacement_max) { #if debug this->model->addDumpFieldVector("displacement"); this->model->addDumpFieldVector("velocity"); this->model->addDumpField("stress"); this->model->addDumpField("strain"); this->model->assembleInternalForces(); this->model->setBaseNameToDumper("cohesive elements", "cohesive_elements"); this->model->addDumpFieldVectorToDumper("cohesive elements", "displacement"); this->model->addDumpFieldToDumper("cohesive elements", "damage"); this->model->addDumpFieldToDumper("cohesive elements", "tractions"); this->model->addDumpFieldToDumper("cohesive elements", "opening"); this->model->dump(); this->model->dump("cohesive elements"); #endif auto inc_load = BC::Dirichlet::Increment(displacement_max / Real(nb_steps)); auto inc_fix = BC::Dirichlet::Increment(-1. / Real(nb_steps) * displacement_max); for (auto _[[gnu::unused]] : arange(nb_steps)) { this->model->applyBC(inc_load, "loading"); this->model->applyBC(inc_fix, "fixed"); if (this->is_extrinsic) this->model->checkCohesiveStress(); this->model->solveStep(); #if debug this->model->dump(); this->model->dump("cohesive elements"); #endif } } void checkInsertion() { auto nb_cohesive_element = this->mesh->getNbElement(cohesive_type); + mesh->getCommunicator().allReduce(nb_cohesive_element, + SynchronizerOperation::_sum); + auto facet_type = this->mesh->getFacetType(this->cohesive_type); const auto & group = this->mesh->getElementGroup("insertion").getElements(facet_type); + auto group_size = group.size(); + mesh->getCommunicator().allReduce(group_size, SynchronizerOperation::_sum); - EXPECT_EQ(nb_cohesive_element, group.size()); + EXPECT_EQ(nb_cohesive_element, group_size); } void checkDissipated(Real expected_density) { Real edis = this->model->getEnergy("dissipated"); EXPECT_NEAR(this->surface * expected_density, edis, 4e-1); } void testModeI() { Vector direction(this->dim, 0.); - direction(this->dim - 1) = 1.; + if (dim == 1) + direction(_x) = 1.; + else + direction(_y) = 1.; // EXPECT_NO_THROW(this->createModel()); this->createModel(); if (this->dim > 1) this->model->applyBC(BC::Dirichlet::FlagOnly(_x), "sides"); if (this->dim > 2) this->model->applyBC(BC::Dirichlet::FlagOnly(_z), "sides"); auto & mat_co = this->model->getMaterial("insertion"); Real sigma_c = mat_co.get("sigma_c"); Real G_c = mat_co.get("G_c"); auto & mat_el = this->model->getMaterial("body"); Real E = mat_el.get("E"); Real nu = mat_el.get("nu"); auto delta_c = 2. * G_c / sigma_c; Real strain = sigma_c / E; if (dim == 1) strain *= .9999; else strain *= .935; // there must be an error in my computations if (this->dim == 2) strain *= (1. - nu) * (1. + nu); auto max_travel = .5 * delta_c; this->setInitialCondition(direction, strain); this->steps(direction * max_travel); } void testModeII() { Vector direction(this->dim, 0.); - direction(0) = 1.; + direction(_x) = 1.; EXPECT_NO_THROW(this->createModel()); if (this->dim > 1) this->model->applyBC(BC::Dirichlet::FlagOnly(_y), "sides"); if (this->dim > 2) this->model->applyBC(BC::Dirichlet::FlagOnly(_z), "sides"); auto & mat_co = this->model->getMaterial("insertion"); Real sigma_c = mat_co.get("sigma_c"); Real G_c = mat_co.get("G_c"); Real beta = mat_co.get("beta"); auto & mat_el = this->model->getMaterial("body"); Real E = mat_el.get("E"); Real nu = mat_el.get("nu"); auto L = this->mesh->getUpperBounds().dot(direction) - this->mesh->getLowerBounds().dot(direction); auto delta_c = 2. * G_c / sigma_c; Real strain = .99999 * L * beta * beta * sigma_c / E; - if(this->dim > 1) { + if (this->dim > 1) { strain *= (1. + nu); } auto max_travel = 1.2 * delta_c; this->setInitialCondition(direction, strain); this->steps(direction * max_travel); } protected: std::unique_ptr mesh; std::unique_ptr model; std::string mesh_name{aka::to_string(cohesive_type) + aka::to_string(type_1) + (type_1 == type_2 ? "" : aka::to_string(type_2)) + ".msh"}; bool is_extrinsic; AnalysisMethod analysis_method; Vector normal; Real surface{0}; UInt nb_steps{300}; }; /* -------------------------------------------------------------------------- */ template constexpr ElementType TestSMMCFixture::cohesive_type; template constexpr ElementType TestSMMCFixture::type_1; template constexpr ElementType TestSMMCFixture::type_2; template constexpr size_t TestSMMCFixture::dim; /* -------------------------------------------------------------------------- */ using IsExtrinsicTypes = std::tuple; using AnalysisMethodTypes = std::tuple>; using types = gtest_list_t, element_type_t<_segment_2>, element_type_t<_segment_2>>, std::tuple, element_type_t<_triangle_3>, element_type_t<_triangle_3>>, std::tuple, element_type_t<_quadrangle_4>, element_type_t<_quadrangle_4>>, std::tuple, element_type_t<_triangle_3>, element_type_t<_quadrangle_4>>, std::tuple, element_type_t<_triangle_6>, element_type_t<_triangle_6>>, std::tuple, element_type_t<_quadrangle_8>, element_type_t<_quadrangle_8>>, std::tuple, element_type_t<_triangle_6>, element_type_t<_quadrangle_8>>, std::tuple, element_type_t<_tetrahedron_4>, element_type_t<_tetrahedron_4>>, std::tuple, element_type_t<_tetrahedron_10>, element_type_t<_tetrahedron_10>>, std::tuple, element_type_t<_hexahedron_8>, element_type_t<_hexahedron_8>>, std::tuple, element_type_t<_hexahedron_20>, element_type_t<_hexahedron_20>>>>; TYPED_TEST_CASE(TestSMMCFixture, types); #endif /* __AKANTU_TEST_COHESIVE_FIXTURE_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc index df4124325..2be5f8e32 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc @@ -1,314 +1,314 @@ /** * @file test_solid_mechanics_model_cube3d.cc * * @author Guillaume Anciaux * * @date creation: Wed Aug 04 2010 * @date last modification: Thu Aug 06 2015 * * @brief test of the class SolidMechanicsModel on the 3d cube * * @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 "boundary_condition_functor.hh" #include "test_solid_mechanics_model_fixture.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { const Real A = 1e-1; const Real E = 1.; const Real poisson = 3. / 10; const Real lambda = E * poisson / (1 + poisson) / (1 - 2 * poisson); const Real mu = E / 2 / (1. + poisson); const Real rho = 1.; const Real cp = std::sqrt((lambda + 2 * mu) / rho); const Real cs = std::sqrt(mu / rho); const Real c = std::sqrt(E / rho); const Vector k = {.5, 0., 0.}; const Vector psi1 = {0., 0., 1.}; const Vector psi2 = {0., 1., 0.}; const Real knorm = k.norm(); /* -------------------------------------------------------------------------- */ template struct Verification {}; /* -------------------------------------------------------------------------- */ template <> struct Verification<1> { void displacement(Vector & disp, const Vector & coord, Real current_time) { const auto & x = coord(_x); const Real omega = c * k[0]; disp(_x) = A * cos(k[0] * x - omega * current_time); } void velocity(Vector & vel, const Vector & coord, Real current_time) { const auto & x = coord(_x); const Real omega = c * k[0]; vel(_x) = omega * A * sin(k[0] * x - omega * current_time); } }; /* -------------------------------------------------------------------------- */ template <> struct Verification<2> { void displacement(Vector & disp, const Vector & X, Real current_time) { Vector kshear = {k[1], k[0]}; Vector kpush = {k[0], k[1]}; const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; disp = A * (kpush * cos(phase_p) + kshear * cos(phase_s)); } void velocity(Vector & vel, const Vector & X, Real current_time) { Vector kshear = {k[1], k[0]}; Vector kpush = {k[0], k[1]}; const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; vel = A * (kpush * omega_p * cos(phase_p) + kshear * omega_s * cos(phase_s)); } }; /* -------------------------------------------------------------------------- */ template <> struct Verification<3> { void displacement(Vector & disp, const Vector & coord, Real current_time) { const auto & X = coord; Vector kpush = k; Vector kshear1(3); Vector kshear2(3); kshear1.crossProduct(k, psi1); kshear2.crossProduct(k, psi2); const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; disp = A * (kpush * cos(phase_p) + kshear1 * cos(phase_s) + kshear2 * cos(phase_s)); } void velocity(Vector & vel, const Vector & coord, Real current_time) { const auto & X = coord; Vector kpush = k; Vector kshear1(3); Vector kshear2(3); kshear1.crossProduct(k, psi1); kshear2.crossProduct(k, psi2); const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; vel = A * (kpush * omega_p * cos(phase_p) + kshear1 * omega_s * cos(phase_s) + kshear2 * omega_s * cos(phase_s)); } }; /* -------------------------------------------------------------------------- */ template class SolutionFunctor : public BC::Dirichlet::DirichletFunctor { public: SolutionFunctor(Real current_time, SolidMechanicsModel & model) : current_time(current_time), model(model) {} public: static constexpr UInt dim = ElementClass<_type>::getSpatialDimension(); inline void operator()(UInt node, Vector & flags, Vector & primal, const Vector & coord) const { flags(0) = true; auto & vel = model.getVelocity(); auto it = vel.begin(model.getSpatialDimension()); Vector v = it[node]; Verification verif; verif.displacement(primal, coord, current_time); verif.velocity(v, coord, current_time); } private: Real current_time; SolidMechanicsModel & model; }; /* -------------------------------------------------------------------------- */ // This fixture uses somewhat finer meshes representing bars. template class TestSMMFixtureBar : public TestSMMFixture> { using parent = TestSMMFixture>; public: void SetUp() override { - this->mesh_file = "bar" + aka::to_string(this->type) + ".msh"; + this->mesh_file = "../patch_tests/data/bar" + aka::to_string(this->type) + ".msh"; parent::SetUp(); getStaticParser().parse("test_solid_mechanics_model_" "dynamics_material.dat"); auto analysis_method = std::tuple_element_t<1, type_>::value; this->model->initFull(_analysis_method = analysis_method); if (this->dump_paraview) { std::stringstream base_name; base_name << "bar" << analysis_method << this->type; this->model->setBaseName(base_name.str()); this->model->addDumpFieldVector("displacement"); this->model->addDumpField("mass"); this->model->addDumpField("velocity"); this->model->addDumpField("acceleration"); this->model->addDumpFieldVector("external_force"); this->model->addDumpFieldVector("internal_force"); this->model->addDumpField("stress"); this->model->addDumpField("strain"); } time_step = this->model->getStableTimeStep() / 10.; this->model->setTimeStep(time_step); // std::cout << "timestep: " << time_step << std::endl; const auto & position = this->mesh->getNodes(); auto & displacement = this->model->getDisplacement(); auto & velocity = this->model->getVelocity(); constexpr auto dim = parent::spatial_dimension; Verification verif; for (auto && tuple : zip(make_view(position, dim), make_view(displacement, dim), make_view(velocity, dim))) { verif.displacement(std::get<1>(tuple), std::get<0>(tuple), 0.); verif.velocity(std::get<2>(tuple), std::get<0>(tuple), 0.); } if (dump_paraview) this->model->dump(); /// boundary conditions this->model->applyBC(SolutionFunctor(0., *this->model), "BC"); wave_velocity = 1.; // sqrt(E/rho) = sqrt(1/1) = 1 simulation_time = 5 / wave_velocity; max_steps = simulation_time / time_step; // 100 auto ndump = 50; dump_freq = max_steps / ndump; } protected: Real time_step; Real wave_velocity; Real simulation_time; UInt max_steps; UInt dump_freq; bool dump_paraview{false}; }; template using analysis_method_t = std::integral_constant; #ifdef AKANTU_IMPLICIT using TestAnalysisTypes = std::tuple, analysis_method_t<_explicit_lumped_mass>>; #else using TestAnalysisTypes = std::tuple>; #endif using TestTypes = gtest_list_t>; TYPED_TEST_CASE(TestSMMFixtureBar, TestTypes); /* -------------------------------------------------------------------------- */ TYPED_TEST(TestSMMFixtureBar, DynamicsExplicit) { constexpr auto dim = TestFixture::spatial_dimension; Real current_time = 0.; const auto & position = this->mesh->getNodes(); const auto & displacement = this->model->getDisplacement(); UInt nb_nodes = this->mesh->getNbNodes(); UInt nb_global_nodes = this->mesh->getNbGlobalNodes(); Real max_error{0.}; Array displacement_solution(nb_nodes, dim); Verification verif; for (UInt s = 0; s < this->max_steps; ++s, current_time += this->time_step) { if (s % this->dump_freq == 0 && this->dump_paraview) this->model->dump(); /// boundary conditions this->model->applyBC( SolutionFunctor(current_time, *this->model), "BC"); // compute the disp solution for (auto && tuple : zip(make_view(position, dim), make_view(displacement_solution, dim))) { verif.displacement(std::get<1>(tuple), std::get<0>(tuple), current_time); } // compute the error solution Real disp_error = 0.; for (auto && tuple : zip(make_view(displacement, dim), make_view(displacement_solution, dim))) { auto diff = std::get<1>(tuple) - std::get<0>(tuple); disp_error += diff.dot(diff); } this->mesh->getCommunicator().allReduce(disp_error, SynchronizerOperation::_sum); disp_error = sqrt(disp_error) / nb_global_nodes; max_error = std::max(disp_error, max_error); ASSERT_NEAR(disp_error, 0., 1e-3); this->model->solveStep(); } // std::cout << "max error: " << max_error << std::endl; } } diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh index 57400aa38..ff2979628 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh @@ -1,61 +1,55 @@ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "test_gtest_utils.hh" #include "communicator.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TEST_SOLID_MECHANICS_MODEL_FIXTURE_HH__ #define __AKANTU_TEST_SOLID_MECHANICS_MODEL_FIXTURE_HH__ using namespace akantu; // This fixture uses very small meshes with a volume of 1. template class TestSMMFixture : public ::testing::Test { public: static constexpr ElementType type = type_::value; static constexpr size_t spatial_dimension = ElementClass::getSpatialDimension(); void SetUp() override { this->mesh = std::make_unique(this->spatial_dimension); -#if defined(AKANTU_PARALLEL) if(Communicator::getStaticCommunicator().whoAmI() == 0) { -#endif - this->mesh->read(this->mesh_file); - -#if defined(AKANTU_PARALLEL) } mesh->distribute(); -#endif SCOPED_TRACE(aka::to_string(this->type).c_str()); model = std::make_unique(*mesh, _all_dimensions, aka::to_string(this->type)); } void TearDown() override { model.reset(nullptr); mesh.reset(nullptr); } protected: std::string mesh_file{aka::to_string(this->type) + ".msh"}; std::unique_ptr mesh; std::unique_ptr model; }; template constexpr ElementType TestSMMFixture::type; template constexpr size_t TestSMMFixture::spatial_dimension; using gtest_element_types = gtest_list_t; TYPED_TEST_CASE(TestSMMFixture, gtest_element_types); #endif /* __AKANTU_TEST_SOLID_MECHANICS_MODEL_FIXTURE_HH__ */