diff --git a/src/mesh/element_type_map.cc b/src/mesh/element_type_map.cc
index 10ed4c342..ee0fc8cca 100644
--- a/src/mesh/element_type_map.cc
+++ b/src/mesh/element_type_map.cc
@@ -1,59 +1,71 @@
/**
* @file element_type_map.cc
*
* @author Nicolas Richart
*
* @date creation Tue Jun 27 2017
*
* @brief A Documented file.
*
* @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 "fe_engine.hh"
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
-FEEngineElementTypeMapArrayInializer::FEEngineElementTypeMapArrayInializer(
+FEEngineElementTypeMapArrayInitializer::FEEngineElementTypeMapArrayInitializer(
const FEEngine & fe_engine, UInt nb_component, UInt spatial_dimension,
const GhostType & ghost_type, const ElementKind & element_kind)
- : MeshElementTypeMapArrayInializer(
- fe_engine.getMesh(),
- nb_component,
+ : MeshElementTypeMapArrayInitializer(
+ fe_engine.getMesh(), nb_component,
spatial_dimension == UInt(-2)
? fe_engine.getMesh().getSpatialDimension()
: spatial_dimension,
ghost_type, element_kind, true, false),
fe_engine(fe_engine) {}
-UInt FEEngineElementTypeMapArrayInializer::size(
+FEEngineElementTypeMapArrayInitializer::FEEngineElementTypeMapArrayInitializer(
+ const FEEngine & fe_engine,
+ const ElementTypeMapArrayInitializer::CompFunc & nb_component,
+ UInt spatial_dimension, const GhostType & ghost_type,
+ const ElementKind & element_kind)
+ : MeshElementTypeMapArrayInitializer(
+ fe_engine.getMesh(), nb_component,
+ spatial_dimension == UInt(-2)
+ ? fe_engine.getMesh().getSpatialDimension()
+ : spatial_dimension,
+ ghost_type, element_kind, true, false),
+ fe_engine(fe_engine) {}
+
+UInt FEEngineElementTypeMapArrayInitializer::size(
const ElementType & type) const {
- return MeshElementTypeMapArrayInializer::size(type) *
+ return MeshElementTypeMapArrayInitializer::size(type) *
fe_engine.getNbIntegrationPoints(type, this->ghost_type);
}
-FEEngineElementTypeMapArrayInializer::ElementTypesIteratorHelper
-FEEngineElementTypeMapArrayInializer::elementTypes() const {
+FEEngineElementTypeMapArrayInitializer::ElementTypesIteratorHelper
+FEEngineElementTypeMapArrayInitializer::elementTypes() const {
return this->fe_engine.elementTypes(spatial_dimension, ghost_type,
element_kind);
}
} // namespace akantu
diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh
index 2b132a0ee..f3aefd368 100644
--- a/src/mesh/element_type_map.hh
+++ b/src/mesh/element_type_map.hh
@@ -1,444 +1,449 @@
/**
* @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(nb_component_functor);
DECLARE_NAMED_ARGUMENT(with_nb_element);
DECLARE_NAMED_ARGUMENT(with_nb_nodes_per_element);
DECLARE_NAMED_ARGUMENT(spatial_dimension);
DECLARE_NAMED_ARGUMENT(do_not_default);
} // 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:
+ /*!
+ * \param _spatial_dimension optional: filter for elements of given spatial
+ * dimension
+ * \param _ghost_type optional: filter for a certain ghost_type
+ * \param _element_kind optional: filter for elements of given kind
+ */
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,
- CompFunc && func);
+ 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 cf1cdc109..23a34e1fe 100644
--- a/src/mesh/element_type_map_tmpl.hh
+++ b/src/mesh/element_type_map_tmpl.hh
@@ -1,696 +1,747 @@
/**
* @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"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__
#define __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__
namespace akantu {
/* -------------------------------------------------------------------------- */
/* ElementTypeMap */
/* -------------------------------------------------------------------------- */
template
inline std::string
ElementTypeMap::printType(const SupportType & type,
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) {
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
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 {
+class ElementTypeMapArrayInitializer {
+protected:
+ using CompFunc = std::function;
+
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),
+ ElementTypeMapArrayInitializer(const CompFunc & comp_func,
+ UInt spatial_dimension = _all_dimensions,
+ const GhostType & ghost_type = _not_ghost,
+ const ElementKind & element_kind = _ek_regular)
+ : comp_func(comp_func), spatial_dimension(spatial_dimension),
ghost_type(ghost_type), element_kind(element_kind) {}
const GhostType & ghostType() const { return ghost_type; }
+ virtual UInt nbComponent(const ElementType & type) const {
+ return comp_func(type, ghostType());
+ }
+
protected:
+ CompFunc comp_func;
UInt spatial_dimension;
- UInt nb_component;
GhostType ghost_type;
ElementKind element_kind;
};
/* -------------------------------------------------------------------------- */
-class MeshElementTypeMapArrayInializer : public ElementTypeMapArrayInializer {
+class MeshElementTypeMapArrayInitializer
+ : public ElementTypeMapArrayInitializer {
+ using CompFunc = ElementTypeMapArrayInitializer::CompFunc;
+
public:
- MeshElementTypeMapArrayInializer(
+ MeshElementTypeMapArrayInitializer(
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),
+ : MeshElementTypeMapArrayInitializer(
+ mesh,
+ [nb_component](const ElementType &, const GhostType &) -> UInt {
+ return nb_component;
+ },
+ spatial_dimension, ghost_type, element_kind, with_nb_element,
+ with_nb_nodes_per_element) {}
+
+ MeshElementTypeMapArrayInitializer(
+ const Mesh & mesh, const CompFunc & comp_func,
+ 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)
+ : ElementTypeMapArrayInitializer(comp_func, spatial_dimension, ghost_type,
+ element_kind),
mesh(mesh), with_nb_element(with_nb_element),
with_nb_nodes_per_element(with_nb_nodes_per_element) {}
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 {
+ UInt nbComponent(const ElementType & type) const override {
+ UInt res = ElementTypeMapArrayInitializer::nbComponent(type);
if (with_nb_nodes_per_element)
- return (this->nb_component * mesh.getNbNodesPerElement(type));
+ return (res * mesh.getNbNodesPerElement(type));
- return this->nb_component;
+ return res;
}
protected:
const Mesh & mesh;
bool with_nb_element;
bool with_nb_nodes_per_element;
};
/* -------------------------------------------------------------------------- */
-class FEEngineElementTypeMapArrayInializer
- : public MeshElementTypeMapArrayInializer {
+class FEEngineElementTypeMapArrayInitializer
+ : public MeshElementTypeMapArrayInitializer {
public:
- FEEngineElementTypeMapArrayInializer(
+ FEEngineElementTypeMapArrayInitializer(
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);
+ FEEngineElementTypeMapArrayInitializer(
+ const FEEngine & fe_engine,
+ const ElementTypeMapArrayInitializer::CompFunc & nb_component,
+ 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
+template
void ElementTypeMapArray::initialize(const Func & f,
const T & default_value,
- bool do_not_default,
- CompFunc && comp_func) {
+ bool do_not_default) {
auto ghost_type = f.ghostType();
for (auto & type : f.elementTypes()) {
if (not this->exists(type, ghost_type))
if (do_not_default) {
- auto & array =
- this->alloc(0, comp_func(type, ghost_type), type, ghost_type);
+ auto & array = this->alloc(0, f.nbComponent(type), type, ghost_type);
array.resize(f.size(type));
} else {
- this->alloc(f.size(type), comp_func(type, ghost_type), type, ghost_type,
+ this->alloc(f.size(type), f.nbComponent(type), type, ghost_type,
default_value);
}
else {
auto & array = this->operator()(type, ghost_type);
if (not do_not_default)
array.resize(f.size(type), default_value);
else
array.resize(f.size(type));
}
}
}
/* -------------------------------------------------------------------------- */
+/**
+ * All parameters are named optionals
+ * \param _nb_component a functor giving the number of components per
+ * (ElementType, GhostType) pair or a scalar giving a unique number of
+ * components
+ * regardless of type
+ * \param _spatial_dimension a filter for the elements of a specific dimension
+ * \param _element_kind filter with element kind (_ek_regular, _ek_structural,
+ * ...)
+ * \param _with_nb_element allocate the arrays with the number of elements for
+ * each
+ * type in the mesh
+ * \param _with_nb_nodes_per_element multiply the number of components by the
+ * number of nodes per element
+ * \param _default_value default inital value
+ * \param _do_not_default do not initialize the allocated arrays
+ * \param _ghost_type filter a type of ghost
+ */
template
template
void ElementTypeMapArray::initialize(const Mesh & mesh,
pack &&... _pack) {
GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _casper);
bool all_ghost_types = requested_ghost_type == _casper;
for (auto ghost_type : ghost_types) {
if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types))
continue;
- auto functor = MeshElementTypeMapArrayInializer(
- mesh, OPTIONAL_NAMED_ARG(nb_component, 1),
+ // This thing should have a UInt or functor type
+ auto && nb_components = OPTIONAL_NAMED_ARG(nb_component, 1);
+ auto functor = MeshElementTypeMapArrayInitializer(
+ mesh, std::forward(nb_components),
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));
- std::function
- nb_component_functor =
- [&](const ElementType & type, const GhostType &
- /*ghost_type*/) -> UInt { return functor.nbComponent(type); };
-
- this->initialize(
- functor, OPTIONAL_NAMED_ARG(default_value, T()),
- OPTIONAL_NAMED_ARG(do_not_default, false),
- OPTIONAL_NAMED_ARG(nb_component_functor, nb_component_functor));
+ this->initialize(functor, OPTIONAL_NAMED_ARG(default_value, T()),
+ OPTIONAL_NAMED_ARG(do_not_default, false));
}
}
/* -------------------------------------------------------------------------- */
+/**
+ * All parameters are named optionals
+ * \param _nb_component a functor giving the number of components per
+ * (ElementType, GhostType) pair or a scalar giving a unique number of
+ * components
+ * regardless of type
+ * \param _spatial_dimension a filter for the elements of a specific dimension
+ * \param _element_kind filter with element kind (_ek_regular, _ek_structural,
+ * ...)
+ * \param _default_value default inital value
+ * \param _do_not_default do not initialize the allocated arrays
+ * \param _ghost_type filter a specific ghost type
+ * \param _all_ghost_types get all ghost types
+ */
template
template
void ElementTypeMapArray::initialize(const FEEngine & fe_engine,
pack &&... _pack) {
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;
- auto functor = FEEngineElementTypeMapArrayInializer(
- fe_engine, OPTIONAL_NAMED_ARG(nb_component, 1),
+ auto && nb_components = OPTIONAL_NAMED_ARG(nb_component, 1);
+ auto functor = FEEngineElementTypeMapArrayInitializer(
+ fe_engine, std::forward(nb_components),
OPTIONAL_NAMED_ARG(spatial_dimension, UInt(-2)), ghost_type,
OPTIONAL_NAMED_ARG(element_kind, _ek_regular));
- std::function
- nb_component_functor =
- [&](const ElementType & type, const GhostType &
- /*ghost_type*/) -> UInt { return functor.nbComponent(type); };
-
- this->initialize(
- functor, OPTIONAL_NAMED_ARG(default_value, T()),
- OPTIONAL_NAMED_ARG(do_not_default, false),
- OPTIONAL_NAMED_ARG(nb_component_functor, nb_component_functor));
+ this->initialize(functor, OPTIONAL_NAMED_ARG(default_value, T()),
+ OPTIONAL_NAMED_ARG(do_not_default, false));
}
}
/* -------------------------------------------------------------------------- */
template
inline T & ElementTypeMapArray::
operator()(const Element & element, UInt component) {
return this->operator()(element.type, element.ghost_type)(element.element,
component);
}
/* -------------------------------------------------------------------------- */
template
inline const T & ElementTypeMapArray::
operator()(const Element & element, UInt component) const {
return this->operator()(element.type, element.ghost_type)(element.element,
component);
}
} // namespace akantu
#endif /* __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh
index 09caa4037..4dc96d2e7 100644
--- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh
+++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh
@@ -1,818 +1,803 @@
/**
* @file material_reinforcement_tmpl.hh
*
* @author Lucas Frerot
*
* @date creation: Thu Feb 1 2018
*
* @brief Reinforcement material
*
* @section LICENSE
*
* Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory
* (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_voigthelper.hh"
#include "material_reinforcement.hh"
namespace akantu {
/* -------------------------------------------------------------------------- */
template
MaterialReinforcement::MaterialReinforcement(
EmbeddedInterfaceModel & model, const ID & id)
: Mat(model, 1,
model.getInterfaceMesh(),
model.getFEEngine("EmbeddedInterfaceFEEngine"), id),
emodel(model),
gradu_embedded("gradu_embedded", *this, 1,
model.getFEEngine("EmbeddedInterfaceFEEngine"),
this->element_filter),
directing_cosines("directing_cosines", *this, 1,
model.getFEEngine("EmbeddedInterfaceFEEngine"),
this->element_filter),
pre_stress("pre_stress", *this, 1,
model.getFEEngine("EmbeddedInterfaceFEEngine"),
this->element_filter),
area(1.0), shape_derivatives() {
AKANTU_DEBUG_IN();
this->initialize();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::initialize() {
AKANTU_DEBUG_IN();
this->registerParam("area", area, _pat_parsable | _pat_modifiable,
"Reinforcement cross-sectional area");
this->registerParam("pre_stress", pre_stress, _pat_parsable | _pat_modifiable,
"Uniform pre-stress");
// this->unregisterInternal(this->stress);
// Fool the AvgHomogenizingFunctor
// stress.initialize(dim * dim);
// Reallocate the element filter
this->element_filter.initialize(this->emodel.getInterfaceMesh(),
_spatial_dimension = 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
MaterialReinforcement::~MaterialReinforcement() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::initMaterial() {
Mat::initMaterial();
gradu_embedded.initialize(dim * dim);
pre_stress.initialize(1);
/// We initialise the stuff that is not going to change during the simulation
this->initFilters();
this->allocBackgroundShapeDerivatives();
this->initBackgroundShapeDerivatives();
this->initDirectingCosines();
}
-/* -------------------------------------------------------------------------- */
-
-namespace detail {
- class FilterInitializer : public MeshElementTypeMapArrayInializer {
- public:
- FilterInitializer(EmbeddedInterfaceModel & emodel,
- const GhostType & ghost_type)
- : MeshElementTypeMapArrayInializer(emodel.getMesh(),
- 1, emodel.getSpatialDimension(),
- ghost_type, _ek_regular) {}
-
- UInt size(const ElementType & /*bgtype*/) const override { return 0; }
- };
-}
-
/* -------------------------------------------------------------------------- */
/// Initialize the filter for background elements
template
void MaterialReinforcement::initFilters() {
for (auto gt : ghost_types) {
for (auto && type : emodel.getInterfaceMesh().elementTypes(1, gt)) {
std::string shaped_id = "filter";
if (gt == _ghost)
shaped_id += ":ghost";
auto & background =
background_filter(std::make_unique>(
"bg_" + shaped_id, this->name),
type, gt);
auto & foreground = foreground_filter(
std::make_unique>(shaped_id, this->name),
type, gt);
- foreground->initialize(detail::FilterInitializer(emodel, gt), 0, true);
- background->initialize(detail::FilterInitializer(emodel, gt), 0, true);
+ foreground->initialize(emodel.getMesh(), _nb_component = 1,
+ _ghost_type = gt);
+ background->initialize(emodel.getMesh(), _nb_component = 1,
+ _ghost_type = gt);
// Computing filters
for (auto && bg_type : background->elementTypes(dim, gt)) {
filterInterfaceBackgroundElements(
(*foreground)(bg_type), (*background)(bg_type), bg_type, type, gt);
}
}
}
}
/* -------------------------------------------------------------------------- */
/// Construct a filter for a (interface_type, background_type) pair
template
void MaterialReinforcement::filterInterfaceBackgroundElements(
Array & foreground, Array & background,
const ElementType & type, const ElementType & interface_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
foreground.resize(0);
background.resize(0);
Array & elements =
emodel.getInterfaceAssociatedElements(interface_type, ghost_type);
Array & elem_filter = this->element_filter(interface_type, ghost_type);
for (auto & elem_id : elem_filter) {
Element & elem = elements(elem_id);
if (elem.type == type) {
background.push_back(elem.element);
foreground.push_back(elem_id);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
namespace detail {
- class BackgroundShapeDInitializer : public ElementTypeMapArrayInializer {
+ class BackgroundShapeDInitializer : public ElementTypeMapArrayInitializer {
public:
- BackgroundShapeDInitializer(UInt spatial_dimension,
- FEEngine & engine,
- const ElementType & foreground_type,
- ElementTypeMapArray & filter,
+ BackgroundShapeDInitializer(UInt spatial_dimension, FEEngine & engine,
+ const ElementType & foreground_type,
+ const ElementTypeMapArray & filter,
const GhostType & ghost_type)
- : ElementTypeMapArrayInializer(spatial_dimension, 0, ghost_type,
- _ek_regular) {
+ : ElementTypeMapArrayInitializer(
+ [](const ElementType & bgtype, const GhostType &) {
+ return ShapeFunctions::getShapeDerivativesSize(bgtype);
+ },
+ spatial_dimension, ghost_type, _ek_regular) {
auto nb_quad = engine.getNbIntegrationPoints(foreground_type);
// Counting how many background elements are affected by elements of
// interface_type
for (auto type : filter.elementTypes(this->spatial_dimension)) {
// Inserting size
array_size_per_bg_type(filter(type).size() * nb_quad, type,
this->ghost_type);
}
}
auto elementTypes() const -> decltype(auto) {
return array_size_per_bg_type.elementTypes();
}
UInt size(const ElementType & bgtype) const {
return array_size_per_bg_type(bgtype, this->ghost_type);
}
- UInt nbComponent(const ElementType & bgtype) const {
- return ShapeFunctions::getShapeDerivativesSize(bgtype);
- }
-
protected:
ElementTypeMap array_size_per_bg_type;
};
}
/**
* Background shape derivatives need to be stored per background element
* types but also per embedded element type, which is why they are stored
* in an ElementTypeMap *>. The outer ElementTypeMap
* refers to the embedded types, and the inner refers to the background types.
*/
template
void MaterialReinforcement::allocBackgroundShapeDerivatives() {
AKANTU_DEBUG_IN();
for (auto gt : ghost_types) {
for (auto && type : emodel.getInterfaceMesh().elementTypes(1, gt)) {
std::string shaped_id = "embedded_shape_derivatives";
if (gt == _ghost)
shaped_id += ":ghost";
auto & shaped_etma = shape_derivatives(
std::make_unique>(shaped_id, this->name),
type, gt);
shaped_etma->initialize(
detail::BackgroundShapeDInitializer(
emodel.getSpatialDimension(),
emodel.getFEEngine("EmbeddedInterfaceFEEngine"), type,
*background_filter(type, gt), gt),
0, true);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::initBackgroundShapeDerivatives() {
AKANTU_DEBUG_IN();
for (auto interface_type :
this->element_filter.elementTypes(this->spatial_dimension)) {
for (auto type : background_filter(interface_type)->elementTypes(dim)) {
computeBackgroundShapeDerivatives(interface_type, type, _not_ghost,
this->element_filter(interface_type));
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::computeBackgroundShapeDerivatives(
const ElementType & interface_type, const ElementType & bg_type,
GhostType ghost_type, const Array & filter) {
auto & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine");
auto & engine = emodel.getFEEngine();
auto & interface_mesh = emodel.getInterfaceMesh();
const auto nb_nodes_elem_bg = Mesh::getNbNodesPerElement(bg_type);
// const auto nb_strss = VoigtHelper::size;
const auto nb_quads_per_elem =
interface_engine.getNbIntegrationPoints(interface_type);
Array quad_pos(0, dim, "interface_quad_pos");
interface_engine.interpolateOnIntegrationPoints(interface_mesh.getNodes(),
quad_pos, dim, interface_type,
ghost_type, filter);
auto & background_shapesd =
(*shape_derivatives(interface_type, ghost_type))(bg_type, ghost_type);
auto & background_elements =
(*background_filter(interface_type, ghost_type))(bg_type, ghost_type);
auto & foreground_elements =
(*foreground_filter(interface_type, ghost_type))(bg_type, ghost_type);
auto shapesd_begin =
background_shapesd.begin(dim, nb_nodes_elem_bg, nb_quads_per_elem);
auto quad_begin = quad_pos.begin(dim, nb_quads_per_elem);
for (auto && tuple : zip(background_elements, foreground_elements)) {
UInt bg = std::get<0>(tuple), fg = std::get<1>(tuple);
for (UInt i = 0; i < nb_quads_per_elem; ++i) {
Matrix shapesd = Tensor3(shapesd_begin[fg])(i);
Vector quads = Matrix(quad_begin[fg])(i);
engine.computeShapeDerivatives(quads, bg, bg_type, shapesd,
ghost_type);
}
}
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::initDirectingCosines() {
AKANTU_DEBUG_IN();
Mesh & mesh = emodel.getInterfaceMesh();
Mesh::type_iterator type_it = mesh.firstType(1, _not_ghost);
Mesh::type_iterator type_end = mesh.lastType(1, _not_ghost);
const UInt voigt_size = VoigtHelper::size;
directing_cosines.initialize(voigt_size);
for (; type_it != type_end; ++type_it) {
computeDirectingCosines(*type_it, _not_ghost);
// computeDirectingCosines(*type_it, _ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::assembleStiffnessMatrix(
GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & interface_mesh = emodel.getInterfaceMesh();
Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost);
Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost);
for (; type_it != type_end; ++type_it) {
assembleStiffnessMatrix(*type_it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::assembleInternalForces(
GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & interface_mesh = emodel.getInterfaceMesh();
Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost);
Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost);
for (; type_it != type_end; ++type_it) {
this->assembleInternalForces(*type_it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::computeAllStresses(GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh::type_iterator it = emodel.getInterfaceMesh().firstType();
Mesh::type_iterator last_type = emodel.getInterfaceMesh().lastType();
for (; it != last_type; ++it) {
computeGradU(*it, ghost_type);
this->computeStress(*it, ghost_type);
addPrestress(*it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::addPrestress(const ElementType & type,
GhostType ghost_type) {
auto & stress = this->stress(type, ghost_type);
auto & sigma_p = this->pre_stress(type, ghost_type);
for (auto && tuple : zip(stress, sigma_p)) {
std::get<0>(tuple) += std::get<1>(tuple);
}
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::assembleInternalForces(
const ElementType & type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & mesh = emodel.getMesh();
Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type);
Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type);
for (; type_it != type_end; ++type_it) {
assembleInternalForcesInterface(type, *type_it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Computes and assemble the residual. Residual in reinforcement is computed as:
*
* \f[
* \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s}
* \f]
*/
template
void MaterialReinforcement::assembleInternalForcesInterface(
const ElementType & interface_type, const ElementType & background_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt voigt_size = VoigtHelper::size;
FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine");
Array & elem_filter = this->element_filter(interface_type, ghost_type);
UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type);
UInt nb_quadrature_points =
interface_engine.getNbIntegrationPoints(interface_type, ghost_type);
UInt nb_element = elem_filter.size();
UInt back_dof = dim * nodes_per_background_e;
Array & shapesd = (*shape_derivatives(interface_type, ghost_type))(
background_type, ghost_type);
Array integrant(nb_quadrature_points * nb_element, back_dof,
"integrant");
Array::vector_iterator integrant_it = integrant.begin(back_dof);
Array::vector_iterator integrant_end = integrant.end(back_dof);
Array::matrix_iterator B_it =
shapesd.begin(dim, nodes_per_background_e);
Array::vector_iterator C_it =
directing_cosines(interface_type, ghost_type).begin(voigt_size);
auto sigma_it = this->stress(interface_type, ghost_type).begin();
Matrix Bvoigt(voigt_size, back_dof);
for (; integrant_it != integrant_end;
++integrant_it, ++B_it, ++C_it, ++sigma_it) {
VoigtHelper::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt,
nodes_per_background_e);
Vector & C = *C_it;
Vector & BtCt_sigma = *integrant_it;
BtCt_sigma.mul(Bvoigt, C);
BtCt_sigma *= *sigma_it * area;
}
Array residual_interface(nb_element, back_dof, "residual_interface");
interface_engine.integrate(integrant, residual_interface, back_dof,
interface_type, ghost_type, elem_filter);
integrant.resize(0);
Array background_filter(nb_element, 1, "background_filter");
auto & filter =
getBackgroundFilter(interface_type, background_type, ghost_type);
emodel.getDOFManager().assembleElementalArrayLocalArray(
residual_interface, emodel.getInternalForce(), background_type,
ghost_type, -1., filter);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::computeDirectingCosines(
const ElementType & type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & interface_mesh = emodel.getInterfaceMesh();
const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const UInt steel_dof = dim * nb_nodes_per_element;
const UInt voigt_size = VoigtHelper::size;
const UInt nb_quad_points = emodel.getFEEngine("EmbeddedInterfaceFEEngine")
.getNbIntegrationPoints(type, ghost_type);
Array node_coordinates(this->element_filter(type, ghost_type).size(),
steel_dof);
this->emodel.getFEEngine().template extractNodalToElementField(
interface_mesh, interface_mesh.getNodes(), node_coordinates, type,
ghost_type, this->element_filter(type, ghost_type));
Array::matrix_iterator directing_cosines_it =
directing_cosines(type, ghost_type).begin(1, voigt_size);
Array::matrix_iterator node_coordinates_it =
node_coordinates.begin(dim, nb_nodes_per_element);
Array::matrix_iterator node_coordinates_end =
node_coordinates.end(dim, nb_nodes_per_element);
for (; node_coordinates_it != node_coordinates_end; ++node_coordinates_it) {
for (UInt i = 0; i < nb_quad_points; i++, ++directing_cosines_it) {
Matrix & nodes = *node_coordinates_it;
Matrix & cosines = *directing_cosines_it;
computeDirectingCosinesOnQuad(nodes, cosines);
}
}
// Mauro: the directing_cosines internal is defined on the quadrature points
// of each element
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialReinforcement::assembleStiffnessMatrix(
const ElementType & type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Mesh & mesh = emodel.getMesh();
Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type);
Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type);
for (; type_it != type_end; ++type_it) {
assembleStiffnessMatrixInterface(type, *type_it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001)
* \f[
* \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T
* \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}}
* \f]
*/
template
void MaterialReinforcement::assembleStiffnessMatrixInterface(
const ElementType & interface_type, const ElementType & background_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt voigt_size = VoigtHelper::size;
FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine");
Array & elem_filter = this->element_filter(interface_type, ghost_type);
Array & grad_u = gradu_embedded(interface_type, ghost_type);
UInt nb_element = elem_filter.size();
UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type);
UInt nb_quadrature_points =
interface_engine.getNbIntegrationPoints(interface_type, ghost_type);
UInt back_dof = dim * nodes_per_background_e;
UInt integrant_size = back_dof;
grad_u.resize(nb_quadrature_points * nb_element);
Array tangent_moduli(nb_element * nb_quadrature_points, 1,
"interface_tangent_moduli");
this->computeTangentModuli(interface_type, tangent_moduli, ghost_type);
Array & shapesd = (*shape_derivatives(interface_type, ghost_type))(
background_type, ghost_type);
Array integrant(nb_element * nb_quadrature_points,
integrant_size * integrant_size, "B^t*C^t*D*C*B");
/// Temporary matrices for integrant product
Matrix Bvoigt(voigt_size, back_dof);
Matrix DCB(1, back_dof);
Matrix CtDCB(voigt_size, back_dof);
Array::scalar_iterator D_it = tangent_moduli.begin();
Array::scalar_iterator D_end = tangent_moduli.end();
Array::matrix_iterator C_it =
directing_cosines(interface_type, ghost_type).begin(1, voigt_size);
Array