diff --git a/src/model/common/constitutive_laws/constitutive_law.hh b/src/model/common/constitutive_laws/constitutive_law.hh
index 4edafec6f..4c4c10018 100644
--- a/src/model/common/constitutive_laws/constitutive_law.hh
+++ b/src/model/common/constitutive_laws/constitutive_law.hh
@@ -1,244 +1,255 @@
/**
* @file constitutive_law.hh
*
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Mother class for all constitutive laws
*
*
* Copyright (©) 2010-2018 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "data_accessor.hh"
#include "mesh_events.hh"
#include "parsable.hh"
#include "parser.hh"
/* -------------------------------------------------------------------------- */
#include "internal_field.hh"
#include "random_internal_field.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_CONSTITUTIVE_LAW_HH__
#define __AKANTU_CONSTITUTIVE_LAW_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
class ConstitutiveLawInternalHandler {
public:
- ConstitutiveLawInternalHandler(const ID & id) : id(id) {}
+ ConstitutiveLawInternalHandler(const ID & id, UInt dim) :
+ id(id), spatial_dimension(dim) {}
template <typename T>
inline void registerInternal(std::shared_ptr<InternalField<T>> & vect);
inline void unregisterInternal(const ID & id);
/// resize the internals arrrays
virtual void resizeInternals();
public:
/// save the internals in the previous_state if needed
virtual void savePreviousState();
/// restore the internals from previous_state if needed
virtual void restorePreviousState();
template <typename T>
const InternalField<T> & getInternal(const ID & id) const;
template <typename T> InternalField<T> & getInternal(const ID & id);
template <typename T>
inline bool isInternal(const ID & id,
const ElementKind & element_kind) const;
template <typename T>
const Array<T> & getArray(const ID & id, ElementType type,
GhostType ghost_type = _not_ghost) const;
template <typename T>
Array<T> & getArray(const ID & id, ElementType type,
GhostType ghost_type = _not_ghost);
public:
- virtual FEEngine & getFEEngine() = 0;
- virtual UInt getSpatialDimension() = 0;
+ virtual const FEEngine & getFEEngine(const ID & id = "") const = 0;
+ virtual FEEngine & getFEEngine(const ID & id = "") = 0;
+ UInt getSpatialDimension() {return spatial_dimension; }
+ virtual const ElementTypeMapArray<UInt> & getElementFilter() const = 0;
AKANTU_GET_MACRO(Name, name, const std::string &);
AKANTU_GET_MACRO(ID, id, const ID &);
+
+
private:
std::map<ID, std::shared_ptr<InternalFieldBase>> internal_vectors;
protected:
ID id;
+ // spatial dimension for constitutive law
+ UInt spatial_dimension;
+
/// constitutive law name
std::string name;
};
template <class ConstitutiveLawsHandler_>
class ConstitutiveLaw : public ConstitutiveLawInternalHandler,
public DataAccessor<Element>,
public MeshEventHandler,
public Parsable {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
using ConstitutiveLawsHandler = ConstitutiveLawsHandler_;
ConstitutiveLaw(const ConstitutiveLaw & law) = delete;
ConstitutiveLaw & operator=(const ConstitutiveLaw & law) = delete;
/// Initialize constitutive law with defaults
ConstitutiveLaw(ConstitutiveLawsHandler & handler, const ID & id = "",
- ElementKind element_kind = _ek_regular);
+ UInt spatial_dimension = _all_dimensions, ElementKind element_kind = _ek_regular);
/// Destructor
~ConstitutiveLaw() override;
protected:
void initialize();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize the constitutive law computed parameter
virtual void initConstitutiveLaw();
/// add an element to the local mesh filter
inline UInt addElement(ElementType type, UInt element, GhostType ghost_type);
inline UInt addElement(const Element & element);
/// add many elements at once
void addElements(const Array<Element> & elements_to_add);
/// remove many element at once
void removeElements(const Array<Element> & elements_to_remove);
/// function to print the contain of the class
void printself(std::ostream & stream, int indent = 0) const override;
protected:
/// function called to update the internal parameters when the
/// modifiable parameters are modified
virtual void updateInternalParameters() {}
/// converts global element to local element
inline Element convertToLocalElement(const Element & global_element) const;
/// converts local element to global element
inline Element convertToGlobalElement(const Element & local_element) const;
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
template <typename T>
inline void packElementDataHelper(const ElementTypeMapArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID()) const;
template <typename T>
inline void unpackElementDataHelper(ElementTypeMapArray<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID());
/* ------------------------------------------------------------------------ */
/* MeshEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
virtual void onNodesAdded(const Array<UInt> & /*unused*/,
const NewNodesEvent & /*unused*/) override{};
virtual void onNodesRemoved(const Array<UInt> & /*unused*/,
const Array<UInt> & /*unused*/,
const RemovedNodesEvent & /*unused*/) override{};
virtual void
onElementsChanged(const Array<Element> & /*unused*/,
const Array<Element> & /*unused*/,
const ElementTypeMapArray<UInt> & /*unused*/,
const ChangedElementsEvent & /*unused*/) override{};
virtual void onElementsAdded(const Array<Element> & /*unused*/,
const NewElementsEvent & /*unused*/);
virtual void
onElementsRemoved(const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event);
public:
template <typename T> inline void setParam(const ID & param, T value);
inline const Parameter & getParam(const ID & param) const;
template <typename T>
void flattenInternal(const std::string & field_id,
ElementTypeMapArray<T> & internal_flat,
const GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) const;
/* ------------------------------------------------------------------------
*/
/* Accessors */
/* ------------------------------------------------------------------------
*/
public:
- AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt);
-
+ const FEEngine & getFEEngine(const ID & id = "") const override {return handler.getFEEngine(); }
+ FEEngine & getFEEngine(const ID & id = "") override {return handler.getFEEngine(); }
+ const ElementTypeMapArray<UInt> & getElementFilter() const override {return element_filter;}
+
template <typename T>
ElementTypeMap<UInt> getInternalDataPerElem(const ID & id,
ElementKind element_kind) const;
protected:
bool isInit() const { return is_init; }
/* ------------------------------------------------------------------------
*/
/* Class Members */
/* ------------------------------------------------------------------------
*/
private:
/// boolean to know if the constitutive law has been initialized
bool is_init{false};
protected:
+
/// list of element handled by the constitutive law
ElementTypeMapArray<UInt> element_filter;
// Constitutive law handler for which this is a constitutive law
ConstitutiveLawsHandler & handler;
};
} // namespace akantu
#include "constitutive_law_tmpl.hh"
#endif
diff --git a/src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh b/src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh
index f56d90afe..d68184c10 100644
--- a/src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh
+++ b/src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh
@@ -1,126 +1,126 @@
/**
* @file constitutive_law_non_local_tmpl.hh
*
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Jul 06 2017
* @date last modification: Tue Nov 07 2017
*
* @brief Implementation of constitutive_law non-local
*
*
* Copyright (©) 2016-2018 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "constitutive_law.hh"
#include "constitutive_law_non_local_interface.hh"
#include "non_local_neighborhood.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
-template <UInt dim, class ConstitutiveLawParent>
-ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::ConstitutiveLawNonLocal(
+template <UInt dim, class ConstitutiveLawNonLocalInterface, class ConstitutiveLawParent>
+ConstitutiveLawNonLocal<dim,ConstitutiveLawNonLocalInterface,ConstitutiveLawParent>::ConstitutiveLawNonLocal(
typename ConstitutiveLawParent::ConstitutiveLawsHandler & handler,
const ID & id)
: ConstitutiveLawParent(handler, id) {
}
/* -------------------------------------------------------------------------- */
-template <UInt dim, class ConstitutiveLawParent>
-void ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::
+template <UInt dim, class ConstitutiveLawNonLocalInterface, class ConstitutiveLawParent>
+void ConstitutiveLawNonLocal<dim, ConstitutiveLawNonLocalInterface, ConstitutiveLawParent>::
insertIntegrationPointsInNeighborhoods(
GhostType ghost_type,
const ElementTypeMapReal & quadrature_points_coordinates) {
IntegrationPoint q;
q.ghost_type = ghost_type;
auto & neighborhood = this->model.getNonLocalManager().getNeighborhood(
this->getNeighborhoodName());
for (auto & type :
this->element_filter.elementTypes(dim, ghost_type, _ek_regular)) {
q.type = type;
const auto & elem_filter = this->element_filter(type, ghost_type);
UInt nb_element = elem_filter.size();
if (nb_element != 0U) {
UInt nb_quad =
this->getFEEngine().getNbIntegrationPoints(type, ghost_type);
const auto & quads = quadrature_points_coordinates(type, ghost_type);
auto nb_total_element =
this->model.getMesh().getNbElement(type, ghost_type);
auto quads_it = quads.begin_reinterpret(dim, nb_quad, nb_total_element);
for (auto & elem : elem_filter) {
Matrix<Real> quads = quads_it[elem];
q.element = elem;
for (UInt nq = 0; nq < nb_quad; ++nq) {
q.num_point = nq;
q.global_num = q.element * nb_quad + nq;
neighborhood.insertIntegrationPoint(q, quads(nq));
}
}
}
}
}
/* -------------------------------------------------------------------------- */
-template <UInt dim, class ConstitutiveLawParent>
-void ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::updateNonLocalInternals(
+template <UInt dim, class ConstitutiveLawNonLocalInterface, class ConstitutiveLawParent>
+void ConstitutiveLawNonLocal<dim, ConstitutiveLawNonLocalInterface, ConstitutiveLawParent>::updateNonLocalInternals(
ElementTypeMapReal & non_local_flattened, const ID & field_id,
GhostType ghost_type, ElementKind kind) {
/// loop over all types in the constitutive_law
for (auto & el_type :
this->element_filter.elementTypes(dim, ghost_type, kind)) {
Array<Real> & internal =
this->template getInternal<Real>(field_id)(el_type, ghost_type);
auto & internal_flat = non_local_flattened(el_type, ghost_type);
auto nb_component = internal_flat.getNbComponent();
auto internal_it = internal.begin(nb_component);
auto internal_flat_it = internal_flat.begin(nb_component);
/// loop all elements for the given type
const auto & filter = this->element_filter(el_type, ghost_type);
UInt nb_quads =
this->getFEEngine().getNbIntegrationPoints(el_type, ghost_type);
for (auto & elem : filter) {
for (UInt q = 0; q < nb_quads; ++q, ++internal_it) {
UInt global_quad = elem * nb_quads + q;
*internal_it = internal_flat_it[global_quad];
}
}
}
}
/* -------------------------------------------------------------------------- */
-template <UInt dim, class ConstitutiveLawParent>
-void ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::registerNeighborhood() {
+template <UInt dim, class ConstitutiveLawNonLocalInterface, class ConstitutiveLawParent>
+void ConstitutiveLawNonLocal<dim, ConstitutiveLawNonLocalInterface, ConstitutiveLawParent>::registerNeighborhood() {
ID name = this->getNeighborhoodName();
this->handler.getNonLocalManager().registerNeighborhood(name, name);
}
} // namespace akantu
diff --git a/src/model/common/constitutive_laws/constitutive_law_tmpl.hh b/src/model/common/constitutive_laws/constitutive_law_tmpl.hh
index 7d1e0c5da..d721a3c41 100644
--- a/src/model/common/constitutive_laws/constitutive_law_tmpl.hh
+++ b/src/model/common/constitutive_laws/constitutive_law_tmpl.hh
@@ -1,511 +1,512 @@
/**
* @file constitutive_law_tmpl.hh *
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Implementation of the templated part of the constitutive law class
*
*
* Copyright (©) 2010-2018 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "constitutive_law.hh"
#include "constitutive_laws_handler.hh"
#include "fe_engine.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_CONSTITUTIVE_LAW_TMPL_HH
#define AKANTU_CONSTITUTIVE_LAW_TMPL_HH
namespace akantu {
/* -------------------------------------------------------------------------- */
template <typename T>
inline void ConstitutiveLawInternalHandler::registerInternal(
std::shared_ptr<InternalField<T>> & vect) {
internal_vectors[vect.getID()] = vect;
}
/* -------------------------------------------------------------------------- */
inline void ConstitutiveLawInternalHandler::unregisterInternal(const ID & id) {
internal_vectors.erase(id);
}
/* -------------------------------------------------------------------------- */
void ConstitutiveLawInternalHandler::savePreviousState() {
for (auto pair : internal_vectors) {
if (pair.second->hasHistory()) {
pair.second->saveCurrentValues();
}
}
}
/* -------------------------------------------------------------------------- */
void ConstitutiveLawInternalHandler::restorePreviousState() {
for (auto pair : internal_vectors) {
if (pair.second->hasHistory()) {
pair.second->restorePreviousValues();
}
}
}
/* -------------------------------------------------------------------------- */
void ConstitutiveLawInternalHandler::resizeInternals() {
for (auto && pair : internal_vectors) {
pair.second->resize();
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
const InternalField<T> &
ConstitutiveLawInternalHandler::getInternal(const ID & id) const {
auto it = internal_vectors.find(this->id + ":" + id);
if (it != internal_vectors.end() and
aka::is_of_type<InternalField<T>>(*it->second)) {
return aka::as_type<InternalField<T>>(*it->second);
}
AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
<< ") does not contain an internal "
<< id << " (" << (getID() + ":" + id)
<< ")");
}
/* -------------------------------------------------------------------------- */
template <typename T>
InternalField<T> & ConstitutiveLawInternalHandler::getInternal(const ID & id) {
auto it = internal_vectors.find(getID() + ":" + id);
if (it != internal_vectors.end() and
aka::is_of_type<InternalField<T>>(*it->second)) {
return aka::as_type<InternalField<T>>(*it->second);
}
AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
<< ") does not contain an internal "
<< id << " (" << (getID() + ":" + id)
<< ")");
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline bool ConstitutiveLawInternalHandler::isInternal(
const ID & id, const ElementKind & element_kind) const {
auto it = internal_vectors.find(this->getID() + ":" + id);
return (it != internal_vectors.end() and
aka::is_of_type<InternalField<T>>(*it->second) and
aka::as_type<InternalField<T>>(*it->second).getElementKind() !=
element_kind);
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
ConstitutiveLaw<ConstitutiveLawsHandler_>::ConstitutiveLaw(
- ConstitutiveLawsHandler_ & handler, const ID & id, ElementKind element_kind)
- : ConstitutiveLawInternalHandler(id),
+ ConstitutiveLawsHandler_ & handler, const ID & id, UInt spatial_dimension,
+ ElementKind element_kind)
+ : ConstitutiveLawInternalHandler(id, spatial_dimension),
Parsable(ParserType::_constitutive_law, id), handler(handler),
element_filter("element_filter", id) {
/// for each connectivity types allocate the element filer array of the
/// constitutive law
element_filter.initialize(handler.getMesh(),
_spatial_dimension = handler.getSpatialDimension(),
_element_kind = element_kind);
this->initialize();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::initialize() {
registerParam("name", name, std::string(), _pat_parsable | _pat_readable);
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::initConstitutiveLaw() {
this->resizeInternals();
this->updateInternalParmaters();
is_init = true;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::addElements(
const Array<Element> & elements_to_add) {
AKANTU_DEBUG_IN();
UInt law_id = handler.getConstitutiveLawIndex(name);
for (const auto & element : elements_to_add) {
auto index = this->addElement(element);
handler.constitutive_law_index(element) = law_id;
handler.constitutive_law_local_numbering(element) = index;
}
this->resizeInternals();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::removeElements(
const Array<Element> & elements_to_remove) {
AKANTU_DEBUG_IN();
auto el_begin = elements_to_remove.begin();
auto el_end = elements_to_remove.end();
if (elements_to_remove.empty()) {
return;
}
auto & mesh = handler.getMesh();
ElementTypeMapArray<UInt> constitutive_law_local_new_numbering(
"remove constitutive law filter elem", id);
constitutive_law_local_new_numbering.initialize(
mesh, _element_filter = &element_filter, _element_kind = _ek_not_defined,
_with_nb_element = true);
ElementTypeMapArray<UInt> element_filter_tmp("element_filter_tmp", id);
element_filter_tmp.initialize(mesh, _element_filter = &element_filter,
_element_kind = _ek_not_defined);
ElementTypeMap<UInt> new_ids, element_ids;
for_each_element(
mesh,
[&](auto && el) {
if (not new_ids(el.type, el.ghost_type)) {
element_ids(el.type, el.ghost_type) = 0;
}
auto & element_id = element_ids(el.type, el.ghost_type);
auto l_el = Element{el.type, element_id, el.ghost_type};
if (std::find(el_begin, el_end, el) != el_end) {
constitutive_law_local_new_numbering(l_el) = UInt(-1);
return;
}
element_filter_tmp(el.type, el.ghost_type).push_back(el.element);
if (not new_ids(el.type, el.ghost_type)) {
new_ids(el.type, el.ghost_type) = 0;
}
auto & new_id = new_ids(el.type, el.ghost_type);
constitutive_law_local_new_numbering(l_el) = new_id;
handler.constitutive_law_local_numbering(el) = new_id;
++new_id;
++element_id;
},
_element_filter = &element_filter, _element_kind = _ek_not_defined);
for (auto ghost_type : ghost_types) {
for (const auto & type : element_filter.elementTypes(
_ghost_type = ghost_type, _element_kind = _ek_not_defined)) {
element_filter(type, ghost_type)
.copy(element_filter_tmp(type, ghost_type));
}
}
for (auto && pair : internal_vectors) {
pair.second->removeIntegrationPoints(constitutive_law_local_new_numbering);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::onElementsAdded(
const Array<Element> & /*unused*/, const NewElementsEvent & /*unused*/) {
this->resizeInternals();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::onElementsRemoved(
const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & /*event*/) {
UInt my_num = handler.getInternalIndexFromID(getID());
ElementTypeMapArray<UInt> constitutive_law_local_new_numbering(
"remove constitutive law filter elem", getID());
auto el_begin = element_list.begin();
auto el_end = element_list.end();
for (auto && gt : ghost_types) {
for (auto && type :
new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) {
if (not element_filter.exists(type, gt) ||
element_filter(type, gt).empty()) {
continue;
}
auto & elem_filter = element_filter(type, gt);
auto & law_indexes = this->handler.constitutive_law_index(type, gt);
auto & law_loc_num =
this->handler.constitutive_law_local_numbering(type, gt);
auto nb_element = this->handler.getMesh().getNbElement(type, gt);
// all constitutive laws will resized to the same size...
law_indexes.resize(nb_element);
law_loc_num.resize(nb_element);
if (!constitutive_law_local_new_numbering.exists(type, gt)) {
constitutive_law_local_new_numbering.alloc(elem_filter.size(), 1, type,
gt);
}
auto & law_renumbering = constitutive_law_local_new_numbering(type, gt);
const auto & renumbering = new_numbering(type, gt);
Array<UInt> 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);
law_renumbering(i) = ni;
law_indexes(new_el) = my_num;
law_loc_num(new_el) = ni;
++ni;
} else {
law_renumbering(i) = UInt(-1);
}
}
elem_filter.resize(elem_filter_tmp.size());
elem_filter.copy(elem_filter_tmp);
}
}
for (auto pair : internal_vectors) {
pair.second->removeIntegrationPoints(constitutive_law_local_new_numbering);
}
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
template <typename T>
inline void ConstitutiveLaw<ConstitutiveLawsHandler_>::packElementDataHelper(
const ElementTypeMapArray<T> & data_to_pack, CommunicationBuffer & buffer,
const Array<Element> & elements, const ID & fem_id) const {
DataAccessor::packElementalDataHelper<T>(data_to_pack, buffer, elements, true,
handler.getFEEngine(fem_id));
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
template <typename T>
inline void ConstitutiveLaw<ConstitutiveLawsHandler_>::unpackElementDataHelper(
ElementTypeMapArray<T> & data_to_unpack, CommunicationBuffer & buffer,
const Array<Element> & elements, const ID & fem_id) {
DataAccessor::unpackElementalDataHelper<T>(data_to_unpack, buffer, elements,
true, handler.getFEEngine(fem_id));
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
inline Element ConstitutiveLaw<ConstitutiveLawsHandler_>::convertToLocalElement(
const Element & global_element) const {
UInt ge = global_element.element;
#ifndef AKANTU_NDEBUG
UInt model_law_index = handler.getConstitutiveLawByElement(
global_element.type, global_element.ghost_type)(ge);
UInt law_index = handler.getConstitutiveLawIndex(this->name);
AKANTU_DEBUG_ASSERT(model_law_index == law_index,
"Conversion of a global element in a local element for "
"the wrong constitutive law "
<< this->name << std::endl);
#endif
UInt le = handler.getConstitutiveLawLocalNumbering(
global_element.type, global_element.ghost_type)(ge);
Element tmp_quad{global_element.type, le, global_element.ghost_type};
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
inline Element
ConstitutiveLaw<ConstitutiveLawsHandler_>::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;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
inline UInt ConstitutiveLaw<ConstitutiveLawsHandler_>::addElement(
ElementType type, UInt element, GhostType ghost_type) {
Array<UInt> & el_filter = this->element_filter(type, ghost_type);
el_filter.push_back(element);
return el_filter.size() - 1;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
inline UInt
ConstitutiveLaw<ConstitutiveLawsHandler_>::addElement(const Element & element) {
return this->addElement(element.type, element.element, element.ghost_type);
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
inline const Parameter &
ConstitutiveLaw<ConstitutiveLawsHandler_>::getParam(const ID & param) const {
try {
return get(param);
} catch (...) {
AKANTU_EXCEPTION("No parameter " << param << " in the constitutive law"
<< getID());
}
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
template <typename T>
inline void
ConstitutiveLaw<ConstitutiveLawsHandler_>::setParam(const ID & param, T value) {
try {
set<T>(param, value);
} catch (...) {
AKANTU_EXCEPTION("No parameter " << param << " in the constitutive law "
<< getID());
}
updateInternalParameters();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
template <typename T>
inline ElementTypeMap<UInt>
ConstitutiveLaw<ConstitutiveLawsHandler_>::getInternalDataPerElem(
const ID & field_id, ElementKind element_kind) const {
if (not this->template isInternal<T>(field_id, element_kind)) {
AKANTU_EXCEPTION("Cannot find internal field "
<< id << " in the constitutive law " << this->name);
}
const auto & internal_field = this->template getInternal<T>(field_id);
const FEEngine & fe_engine = internal_field.getFEEngine();
UInt nb_data_per_quad = internal_field.getNbComponent();
ElementTypeMap<UInt> res;
for (auto ghost_type : ghost_types) {
for (auto & type : internal_field.elementTypes(ghost_type)) {
UInt nb_quadrature_points =
fe_engine.getNbIntegrationPoints(type, ghost_type);
res(type, ghost_type) = nb_data_per_quad * nb_quadrature_points;
}
}
return res;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
template <typename T>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::flattenInternal(
const std::string & field_id, ElementTypeMapArray<T> & internal_flat,
const GhostType ghost_type, ElementKind element_kind) const {
if (!this->template isInternal<T>(field_id, element_kind)) {
AKANTU_EXCEPTION("Cannot find internal field "
<< id << " in the constitutive law " << this->name);
}
const InternalField<T> & internal_field =
this->template getInternal<T>(field_id);
const FEEngine & fe_engine = internal_field.getFEEngine();
const Mesh & mesh = fe_engine.getMesh();
for (auto && type : internal_field.filterTypes(ghost_type)) {
const Array<Real> & src_vect = internal_field(type, ghost_type);
const Array<UInt> & 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<Real> & dst_vect = internal_flat(type, ghost_type);
dst_vect.resize(nb_element_dst * nb_quad_per_elem);
auto it_dst = make_view(dst_vect, nb_data).begin();
for (auto && data : zip(filter, make_view(src_vect, nb_data))) {
it_dst[std::get<0>(data)] = std::get<1>(data);
}
}
}
} // namespace akantu
#endif // AKANTU_CONSTITUTIVE_LAW_TMPL_HH
diff --git a/src/model/common/constitutive_laws/constitutive_laws_handler.hh b/src/model/common/constitutive_laws/constitutive_laws_handler.hh
index aeeca8268..866fb882d 100644
--- a/src/model/common/constitutive_laws/constitutive_laws_handler.hh
+++ b/src/model/common/constitutive_laws/constitutive_laws_handler.hh
@@ -1,292 +1,294 @@
/**
* @file constitutive_laws_handler.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation ven mar 26 2021
*
* @brief A handler for multiple consistutive laws
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "constitutive_law.hh"
#include "constitutive_law_selector.hh"
#include "mesh.hh"
#include "mesh_events.hh"
#include "non_local_manager_callback.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_CONSTITUTIVE_LAWS_HANDLER_HH
#define AKANTU_CONSTITUTIVE_LAWS_HANDLER_HH
namespace akantu {
class NonLocalManager;
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
template <class ConstitutiveLawType>
class ConstitutiveLawsHandler : public DataAccessor<Element>,
public NonLocalManagerCallback {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ConstitutiveLawsHandler(const Mesh & mesh, UInt spatial_dimension,
const ID & parent_id)
: constitutive_law_index("constitutive_law index", parent_id),
constitutive_law_local_numbering("constitutive_law local numbering",
parent_id) {
constitutive_law_selector =
std::make_shared<DefaultConstitutiveLawSelector>(
constitutive_law_index);
constitutive_law_index.initialize(mesh, _element_kind = _ek_not_defined,
_default_value = UInt(-1),
_with_nb_element = true);
constitutive_law_local_numbering.initialize(
mesh, _element_kind = _ek_not_defined, _with_nb_element = true);
}
virtual ~ConstitutiveLawsHandler();
class NewConstitutiveLawElementsEvent : public NewElementsEvent {
public:
AKANTU_GET_MACRO_NOT_CONST(ConstitutiveLawList, constitutive_law,
Array<UInt> &);
AKANTU_GET_MACRO(ConstitutiveLawList, constitutive_law,
const Array<UInt> &);
protected:
Array<UInt> constitutive_law;
};
/* ------------------------------------------------------------------------ */
/* ConstitutiveLaws */
/* ------------------------------------------------------------------------ */
public:
/// register an empty constitutive_law of a given type
auto & registerNewConstitutiveLaw(const ID & cl_name, const ID & cl_type,
const ID & opt_param);
/// reassigns constitutive_laws depending on the constitutive_law selector
virtual void reassignConstitutiveLaw();
template <class Func> void for_each_constitutive_law(Func && func) {
for (auto && constitutive_law : constitutive_laws) {
std::forward<Func>(func)(
aka::as_type<ConstitutiveLawType>(*constitutive_laws));
}
}
protected:
/// register a constitutive_law in the dynamic database
auto & registerNewConstitutiveLaw(const ParserSection & cl_section);
/// read the constitutive_law files to instantiate all the constitutive_laws
void instantiateConstitutiveLaws();
/// set the element_id_by_constitutive_law and add the elements to the good
/// constitutive_laws
virtual void assignConstitutiveLawToElements(
const ElementTypeMapArray<UInt> * filter = nullptr);
protected:
void splitElementByConstitutiveLaw(
const Array<Element> & elements,
std::vector<Array<Element>> & elements_per_cl) const;
template <typename Operation>
void splitByConstitutiveLaw(const Array<Element> & elements,
Operation && op) const;
/* ------------------------------------------------------------------------ */
/* Dumpable interface (kept for convenience) and dumper relative functions */
/* ------------------------------------------------------------------------ */
public:
//! decide wether a field is a constitutive_law internal or not
bool isInternal(const std::string & field_name, ElementKind element_kind);
//! give the amount of data per element
virtual ElementTypeMap<UInt>
getInternalDataPerElem(const std::string & field_name, ElementKind kind);
//! flatten a given constitutive_law internal field
ElementTypeMapArray<Real> &
flattenInternal(const std::string & field_name, ElementKind kind,
GhostType ghost_type = _not_ghost);
//! flatten all the registered constitutive_law internals
void flattenAllRegisteredInternals(ElementKind kind);
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get an iterable on the constitutive_laws
inline decltype(auto) getConstitutiveLaws();
/// get an iterable on the constitutive_laws
inline decltype(auto) getConstitutiveLaws() const;
/// get a particular constitutive_law (by numerical constitutive_law index)
inline auto & getConstitutiveLaw(UInt cl_index);
/// get a particular constitutive_law (by numerical constitutive_law index)
inline const auto & getConstitutiveLaw(UInt cl_index) const;
/// get a particular constitutive_law (by constitutive_law name)
inline auto & getConstitutiveLaw(const std::string & name);
/// get a particular constitutive_law (by constitutive_law name)
inline const auto & getConstitutiveLaw(const std::string & name) const;
/// get a particular constitutive_law id from is name
inline UInt getConstitutiveLawIndex(const std::string & name) const;
/// give the number of constitutive_laws
inline UInt getNbConstitutiveLaws() const { return constitutive_laws.size(); }
/// give the constitutive_law internal index from its id
Int getInternalIndexFromID(const ID & id) const;
AKANTU_GET_MACRO(ConstitutiveLawByElement, constitutive_law_index,
const ElementTypeMapArray<UInt> &);
AKANTU_GET_MACRO(ConstitutiveLawLocalNumbering,
constitutive_law_local_numbering,
const ElementTypeMapArray<UInt> &);
/// vectors containing local constitutive_law element index for each global
/// element index
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConstitutiveLawByElement,
constitutive_law_index, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConstitutiveLawLocalNumbering,
constitutive_law_local_numbering,
UInt);
AKANTU_GET_MACRO_NOT_CONST(ConstitutiveLawSelector,
*constitutive_law_selector,
ConstitutiveLawSelector &);
/// Access the non_local_manager interface
AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &);
void setConstitutiveLawSelector(
std::shared_ptr<ConstitutiveLawSelector> constitutive_law_selector) {
this->constitutive_law_selector = std::move(constitutive_law_selector);
}
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) override;
/* ------------------------------------------------------------------------ */
/* Mesh Event Handler inherited members */
/* ------------------------------------------------------------------------ */
protected:
void onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) override;
void onElementsRemoved(const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) override;
/* ------------------------------------------------------------------------ */
/* NonLocalManager inherited members */
/* ------------------------------------------------------------------------ */
protected:
void initializeNonLocal() override;
void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override;
void computeNonLocalStresses(GhostType ghost_type) override;
void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override;
/// update the values of the non local internal
void updateLocalInternal(ElementTypeMapReal & internal_flat,
GhostType ghost_type, ElementKind kind) override;
/// copy the results of the averaging in the constitutive_laws
void updateNonLocalInternal(ElementTypeMapReal & internal_flat,
GhostType ghost_type, ElementKind kind) override;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// id of the derived class
const ID & parent_id;
/// underlying mesh
const Mesh & _mesh;
/// spatial_dimension of the derived model
UInt _spatial_dimension{_all_dimensions};
/// mapping between constitutive_law name and constitutive_law internal id
std::map<std::string, UInt> constitutive_laws_names_to_id;
/// Arrays containing the constitutive_law index for each element
ElementTypeMapArray<UInt> constitutive_law_index;
/// Arrays containing the position in the element filter of the
/// constitutive_law (constitutive_law's local numbering)
ElementTypeMapArray<UInt> constitutive_law_local_numbering;
/// list of used constitutive_laws
std::vector<std::unique_ptr<ConstitutiveLawType>> constitutive_laws;
/// class defining of to choose a constitutive_law
std::shared_ptr<ConstitutiveLawSelector> constitutive_law_selector;
using flatten_internal_map =
std::map<std::pair<std::string, ElementKind>,
std::unique_ptr<ElementTypeMapArray<Real>>>;
/// map a registered internals to be flattened for dump purposes
flatten_internal_map registered_internals;
/// tells if the constitutive_law are instantiated
bool are_constitutive_laws_instantiated{false};
/// non local manager
std::unique_ptr<NonLocalManager> non_local_manager;
template <class ConstitutiveLawsHandler> friend class ConstitutiveLaw;
};
-
+
} // namespace akantu
+
+#include "constitutive_laws_handler_tmpl.hh"
#endif /* AKANTU_CONSTITUTIVE_LAWS_HANDLER_HH */
diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc
index 0c605b236..3ae13adf4 100644
--- a/src/model/solid_mechanics/material.cc
+++ b/src/model/solid_mechanics/material.cc
@@ -1,1123 +1,1119 @@
/**
* @file material.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Implementation of the common part of the material class
*
*
* Copyright (©) 2010-2018 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material.hh"
#include "mesh_iterators.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
Material::Material(SolidMechanicsModel & model, const ID & id)
- : ConstitutiveLaw<SolidMechanicsModel>(model, id, _ek_regular),
- fem(model.getFEEngine()),
- model(model), spatial_dimension(this->model.getSpatialDimension()),
- 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),
- interpolation_inverse_coordinates("interpolation inverse coordinates",
- *this),
+ : ConstitutiveLaw<SolidMechanicsModel>(model, id, _ek_regular),
+ fem(model.getFEEngine()), model(model), 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),
+ interpolation_inverse_coordinates("interpolation inverse coordinates",
+ *this),
interpolation_points_matrices("interpolation points matrices", *this),
eigen_grad_u(model.getSpatialDimension(), model.getSpatialDimension(),
0.) {
AKANTU_DEBUG_IN();
this->registerParam("eigen_grad_u", eigen_grad_u, _pat_parsable,
"EigenGradU");
this->initialize();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Material::Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh,
FEEngine & fe_engine, const ID & id)
- : ConstitutiveLaw<SolidMechanicsModel>(model, id, _ek_regular),fem(fe_engine),
- model(model), spatial_dimension(dim),
+ : ConstitutiveLaw<SolidMechanicsModel>(model, id, dim, _ek_regular),
+ fem(fe_engine), model(model),
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),
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();
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();
}
Parent::initConstitutiveLaw();
AKANTU_DEBUG_OUT();
}
-
/* -------------------------------------------------------------------------- */
void Material::updateInternalParameters() {
auto dim = model.getSpatialDimension();
for (const auto & type :
element_filter.elementTypes(_element_kind = _ek_regular)) {
for (auto eigen_gradu : make_view(eigengradu(type), dim, dim)) {
eigen_gradu = eigen_grad_u;
}
}
Parent::updateInternalParameters();
}
-
+
/* -------------------------------------------------------------------------- */
/**
* Compute the internal forces by assembling @f$\int_{e} \sigma_e \frac{\partial
* \varphi}{\partial X} dX @f$
*
* @param[in] ghost_type compute the internal forces for _ghost or _not_ghost
* element
*/
void Material::assembleInternalForces(GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = model.getSpatialDimension();
if (!finite_deformation) {
auto & internal_force = const_cast<Array<Real> &>(model.getInternalForce());
// Mesh & mesh = fem.getMesh();
for (auto && type :
element_filter.elementTypes(spatial_dimension, ghost_type)) {
Array<UInt> & elem_filter = element_filter(type, ghost_type);
UInt nb_element = elem_filter.size();
if (nb_element == 0) {
continue;
}
const Array<Real> & 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$
auto * sigma_dphi_dx =
new Array<Real>(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);
/**
* 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$
*/
auto * int_sigma_dphi_dx =
new Array<Real>(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] 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<UInt> & elem_filter = element_filter(type, ghost_type);
if (elem_filter.empty()) {
continue;
}
Array<Real> & 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->StoCauchy<1>(type, ghost_type);
break;
case 2:
this->StoCauchy<2>(type, ghost_type);
break;
case 3:
this->StoCauchy<3>(type, ghost_type);
break;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void Material::StoCauchy(ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
auto gradu_it = this->gradu(el_type, ghost_type).begin(dim, dim);
auto gradu_end = this->gradu(el_type, ghost_type).end(dim, dim);
auto piola_it = this->piola_kirchhoff_2(el_type, ghost_type).begin(dim, dim);
auto stress_it = this->stress(el_type, ghost_type).begin(dim, dim);
for (; gradu_it != gradu_end; ++gradu_it, ++piola_it, ++stress_it) {
Matrix<Real> & grad_u = *gradu_it;
Matrix<Real> & piola = *piola_it;
Matrix<Real> & sigma = *stress_it;
auto F_tensor = gradUToF<dim>(grad_u);
this->StoCauchy<dim>(F_tensor, piola, sigma);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::setToSteadyState(GhostType ghost_type) {
AKANTU_DEBUG_IN();
const Array<Real> & displacement = model.getDisplacement();
// resizeInternalArray(gradu);
UInt spatial_dimension = model.getSpatialDimension();
for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) {
Array<UInt> & elem_filter = element_filter(type, ghost_type);
Array<Real> & 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] 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 <UInt dim>
void Material::assembleStiffnessMatrix(ElementType type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Array<UInt> & elem_filter = element_filter(type, ghost_type);
if (elem_filter.empty()) {
AKANTU_DEBUG_OUT();
return;
}
// const Array<Real> & shapes_derivatives =
// fem.getShapesDerivatives(type, ghost_type);
Array<Real> & 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);
auto * tangent_stiffness_matrix =
new Array<Real>(nb_element * nb_quadrature_points,
tangent_size * tangent_size, "tangent_stiffness_matrix");
tangent_stiffness_matrix->zero();
computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type);
/// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
UInt bt_d_b_size = dim * nb_nodes_per_element;
auto * bt_d_b = new Array<Real>(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);
delete tangent_stiffness_matrix;
/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
auto * K_e = new Array<Real>(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 <UInt dim>
void Material::assembleStiffnessMatrixNL(ElementType type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
const Array<Real> & shapes_derivatives =
fem.getShapesDerivatives(type, ghost_type);
Array<UInt> & elem_filter = element_filter(type, ghost_type);
// Array<Real> & 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);
auto * shapes_derivatives_filtered = new Array<Real>(
nb_element * nb_quadrature_points, dim * nb_nodes_per_element,
"shapes derivatives filtered");
FEEngine::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;
auto * bt_s_b = new Array<Real>(nb_element * nb_quadrature_points,
bt_s_b_size * bt_s_b_size, "B^t*D*B");
UInt piola_matrix_size = getCauchyStressMatrixSize(dim);
Matrix<Real> B(piola_matrix_size, bt_s_b_size);
Matrix<Real> Bt_S(bt_s_b_size, piola_matrix_size);
Matrix<Real> 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<dim>(Piola_kirchhoff_matrix, S);
VoigtHelper<dim>::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B,
nb_nodes_per_element);
Bt_S.template mul<true, false>(B, S);
Bt_S_B.template mul<false, false>(Bt_S, B);
}
delete shapes_derivatives_filtered;
/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
auto * K_e = new Array<Real>(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 <UInt dim>
void Material::assembleStiffnessMatrixL2(ElementType type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
const Array<Real> & shapes_derivatives =
fem.getShapesDerivatives(type, ghost_type);
Array<UInt> & elem_filter = element_filter(type, ghost_type);
Array<Real> & 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);
auto * tangent_stiffness_matrix =
new Array<Real>(nb_element * nb_quadrature_points,
tangent_size * tangent_size, "tangent_stiffness_matrix");
tangent_stiffness_matrix->zero();
computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type);
auto * shapes_derivatives_filtered = new Array<Real>(
nb_element * nb_quadrature_points, dim * nb_nodes_per_element,
"shapes derivatives filtered");
FEEngine::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;
auto * bt_d_b = new Array<Real>(nb_element * nb_quadrature_points,
bt_d_b_size * bt_d_b_size, "B^t*D*B");
Matrix<Real> B(tangent_size, dim * nb_nodes_per_element);
Matrix<Real> B2(tangent_size, dim * nb_nodes_per_element);
Matrix<Real> 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<dim > (*shapes_derivatives_filtered_it, B,
// nb_nodes_per_element);
VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(
*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
VoigtHelper<dim>::transferBMatrixToBL2(*shapes_derivatives_filtered_it,
grad_u, B2, nb_nodes_per_element);
B += B2;
Bt_D.template mul<true, false>(B, D);
Bt_D_B.template mul<false, false>(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$
auto * K_e = new Array<Real>(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 <UInt dim>
void Material::assembleInternalForces(GhostType ghost_type) {
AKANTU_DEBUG_IN();
Array<Real> & internal_force = model.getInternalForce();
Mesh & mesh = fem.getMesh();
for (auto type : element_filter.elementTypes(_ghost_type = ghost_type)) {
const Array<Real> & shapes_derivatives =
fem.getShapesDerivatives(type, ghost_type);
Array<UInt> & elem_filter = element_filter(type, ghost_type);
if (elem_filter.empty()) {
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);
auto * shapesd_filtered = new Array<Real>(
nb_element, size_of_shapes_derivatives, "filtered shapesd");
FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered,
type, ghost_type, elem_filter);
Array<Real>::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<Real>(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<Real>::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1);
Matrix<Real> B_tensor(stress_size, bt_s_size);
Matrix<Real> 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 = *bt_s_it;
auto & S = *stress_it;
VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(
*shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element);
VoigtHelper<dim>::transferBMatrixToBL2(*shapes_derivatives_filtered_it,
grad_u, B2_tensor,
nb_nodes_per_element);
B_tensor += B2_tensor;
auto S_vect = Material::stressToVoigt<dim>(S);
Matrix<Real> S_voigt(S_vect.storage(), stress_size, 1);
r.template mul<true, false>(B_tensor, S_voigt);
}
delete shapesd_filtered;
/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
auto * r_e = new Array<Real>(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::computePotentialEnergyByElements() {
AKANTU_DEBUG_IN();
for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) {
computePotentialEnergy(type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::computePotentialEnergy(ElementType /*unused*/) {
AKANTU_DEBUG_IN();
AKANTU_TO_IMPLEMENT();
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<Real> 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<Real> & 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<Real> & 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<Real> & result,
ElementTypeMapArray<Real> & 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<UInt> & elem_fil = element_filter(type, ghost_type);
Array<Real> & 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<Element> & 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<std::vector<Element>> * element_to_facet = nullptr;
GhostType current_ghost_type = _casper;
Array<Real> * result_vec = nullptr;
Array<Real>::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<Real> result_local(result_vec->storage() +
(global_facet * nb_quad_per_facet + q) *
result_vec->getNbComponent() +
static_cast<UInt>(is_second_element) *
stress_size,
stress_size);
const Matrix<Real> & result_tmp(result_it[global_el]);
result_local = result_tmp(f * nb_quad_per_facet + q);
}
}
}
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
const Array<T> & Material::getArray(const ID & vect_id, ElementType type,
GhostType ghost_type) const {
try {
return this->template getInternal<T>(vect_id)(type, ghost_type);
} catch (debug::Exception & e) {
AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
<< ") does not contain a vector "
<< vect_id << " [" << e << "]");
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
Array<T> & Material::getArray(const ID & vect_id, ElementType type,
GhostType ghost_type) {
try {
return this->template getInternal<T>(vect_id)(type, ghost_type);
} catch (debug::Exception & e) {
AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
<< ") does not contain a vector "
<< vect_id << " [" << e << "]");
}
}
template const Array<Real> & Material::getArray(const ID & vect_id,
ElementType type,
GhostType ghost_type) const;
/* -------------------------------------------------------------------------- */
template Array<Real> & Material::getArray(const ID & vect_id, ElementType type,
GhostType ghost_type);
/* -------------------------------------------------------------------------- */
template const Array<UInt> & Material::getArray(const ID & vect_id,
ElementType type,
GhostType ghost_type) const;
/* -------------------------------------------------------------------------- */
template Array<UInt> & Material::getArray(const ID & vect_id, ElementType type,
GhostType ghost_type);
-
/* -------------------------------------------------------------------------- */
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<Element> & /*unused*/,
const NewElementsEvent & /*unused*/) {
this->resizeInternals();
}
/* -------------------------------------------------------------------------- */
void Material::onElementsRemoved(
const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
[[gnu::unused]] const RemovedElementsEvent & event) {
UInt my_num = model.getInternalIndexFromID(getID());
ElementTypeMapArray<UInt> material_local_new_numbering(
"remove mat filter elem", getID());
auto el_begin = element_list.begin();
auto el_end = element_list.end();
for (auto && gt : ghost_types) {
for (auto && type :
new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) {
if (not element_filter.exists(type, gt) ||
element_filter(type, gt).empty()) {
continue;
}
auto & elem_filter = element_filter(type, gt);
auto & mat_indexes = this->model.material_index(type, gt);
auto & mat_loc_num = this->model.material_local_numbering(type, gt);
auto 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);
}
auto & mat_renumbering = material_local_new_numbering(type, gt);
const auto & renumbering = new_numbering(type, gt);
Array<UInt> 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(bool converged) {
if (not converged) {
this->restorePreviousState();
return;
}
for (const auto & type : element_filter.elementTypes(
_all_dimensions, _not_ghost, _ek_not_defined)) {
this->updateEnergies(type);
}
}
/* -------------------------------------------------------------------------- */
void Material::onDamageIteration() { this->savePreviousState(); }
/* -------------------------------------------------------------------------- */
void Material::onDamageUpdate() {
for (const auto & type : element_filter.elementTypes(
_all_dimensions, _not_ghost, _ek_not_defined)) {
this->updateEnergiesAfterDamage(type);
}
}
/* -------------------------------------------------------------------------- */
void Material::onDump() {
if (this->isFiniteDeformation()) {
this->computeAllCauchyStresses(_not_ghost);
}
}
/* -------------------------------------------------------------------------- */
void Material::printself(std::ostream & stream, int indent) const {
std::string space(indent, 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,
[[gnu::unused]] const Matrix<Real> & point,
Matrix<Real> & extrapolated) {
if (this->isInternal<Real>(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<Real> & internal =
this->getArray<Real>(id, element.type, element.ghost_type);
UInt nb_component = internal.getNbComponent();
Array<Real>::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<Real> & values = internal_it[local_element.element];
UInt index = 0;
Vector<Real> 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<Real> default_values(extrapolated.rows(), extrapolated.cols(), 0.);
extrapolated = default_values;
}
}
/* -------------------------------------------------------------------------- */
void Material::applyEigenGradU(const Matrix<Real> & 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).empty()) {
continue;
}
for (auto & eigengradu : make_view(this->eigengradu(type, ghost_type),
spatial_dimension, spatial_dimension)) {
eigengradu = prescribed_eigen_grad_u;
}
}
}
/* -------------------------------------------------------------------------- */
MaterialFactory & Material::getFactory() {
return MaterialFactory::getInstance();
}
} // namespace akantu
diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh
index 7754e2d08..afc1426f2 100644
--- a/src/model/solid_mechanics/material.hh
+++ b/src/model/solid_mechanics/material.hh
@@ -1,587 +1,585 @@
/**
* @file material.hh
*
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Mother class for all materials
*
*
* Copyright (©) 2010-2018 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_factory.hh"
#include "aka_voigthelper.hh"
#include "data_accessor.hh"
#include "integration_point.hh"
#include "parsable.hh"
#include "parser.hh"
#include "constitutive_law.hh"
/* -------------------------------------------------------------------------- */
#include "internal_field.hh"
#include "random_internal_field.hh"
/* -------------------------------------------------------------------------- */
#include "mesh_events.hh"
#include "solid_mechanics_model_event_handler.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_MATERIAL_HH_
#define AKANTU_MATERIAL_HH_
/* -------------------------------------------------------------------------- */
namespace akantu {
class Model;
class SolidMechanicsModel;
class Material;
} // namespace akantu
namespace akantu {
using MaterialFactory =
Factory<Material, ID, UInt, const ID &, SolidMechanicsModel &, const ID &>;
/**
* Interface of all materials
* Prerequisites for a new material
* - inherit from this class
* - implement the following methods:
* \code
* virtual Real getStableTimeStep(Real h, const Element & element =
* ElementNull);
*
* virtual void computeStress(ElementType el_type,
* GhostType ghost_type = _not_ghost);
*
* virtual void computeTangentStiffness(ElementType el_type,
* Array<Real> & tangent_matrix,
* GhostType ghost_type = _not_ghost);
* \endcode
*
*/
class Material : public ConstitutiveLaw<SolidMechanicsModel>,
protected SolidMechanicsModelEventHandler {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Material(const Material & mat) = delete;
Material & operator=(const Material & mat) = delete;
/// Initialize material with defaults
Material(SolidMechanicsModel & model, const ID & id = "");
/// Initialize material with custom mesh & fe_engine
Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh,
FEEngine & fe_engine, const ID & id = "");
/// Destructor
~Material() override;
protected:
void initialize();
/* ------------------------------------------------------------------------ */
/* Function that materials can/should reimplement */
/* ------------------------------------------------------------------------ */
protected:
/// constitutive law
virtual void computeStress(ElementType /* el_type */,
GhostType /* ghost_type */ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// compute the tangent stiffness matrix
virtual void computeTangentModuli(ElementType /*el_type*/,
Array<Real> & /*tangent_matrix*/,
GhostType /*ghost_type*/ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// compute the potential energy
virtual void computePotentialEnergy(ElementType el_type);
/// compute the potential energy for an element
virtual void
computePotentialEnergyByElement(ElementType /*type*/, UInt /*index*/,
Vector<Real> & /*epot_on_quad_points*/) {
AKANTU_TO_IMPLEMENT();
}
virtual void updateEnergies(ElementType /*el_type*/) {}
virtual void updateEnergiesAfterDamage(ElementType /*el_type*/) {}
/// set the material to steady state (to be implemented for materials that
/// need it)
virtual void setToSteadyState(ElementType /*el_type*/,
GhostType /*ghost_type*/ = _not_ghost) {}
/// function called to update the internal parameters when the
/// modifiable parameters are modified
virtual void updateInternalParameters();
public:
/// extrapolate internal values
virtual void extrapolateInternal(const ID & id, const Element & element,
const Matrix<Real> & points,
Matrix<Real> & extrapolated);
/// compute the p-wave speed in the material
virtual Real getPushWaveSpeed(const Element & /*element*/) const {
AKANTU_TO_IMPLEMENT();
}
/// compute the s-wave speed in the material
virtual Real getShearWaveSpeed(const Element & /*element*/) const {
AKANTU_TO_IMPLEMENT();
}
/// get a material celerity to compute the stable time step (default: is the
/// push wave speed)
virtual Real getCelerity(const Element & element) const {
return getPushWaveSpeed(element);
}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize the material computed parameter
virtual void initMaterial();
/// compute the residual for this material
// virtual void updateResidual(GhostType ghost_type = _not_ghost);
/// assemble the residual for this material
virtual void assembleInternalForces(GhostType ghost_type);
/// compute the stresses for this material
virtual void computeAllStresses(GhostType ghost_type = _not_ghost);
// virtual void
// computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost);
virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost);
/// set material to steady state
void setToSteadyState(GhostType ghost_type = _not_ghost);
/// compute the stiffness matrix
virtual void assembleStiffnessMatrix(GhostType ghost_type);
/// function to print the contain of the class
void printself(std::ostream & stream, int indent = 0) const override;
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points
*/
void interpolateStress(ElementTypeMapArray<Real> & result,
GhostType ghost_type = _not_ghost);
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points and store the
* results per facet
*/
void interpolateStressOnFacets(ElementTypeMapArray<Real> & result,
ElementTypeMapArray<Real> & by_elem_result,
GhostType ghost_type = _not_ghost);
/**
* function to initialize the elemental field interpolation
* function by inverting the quadrature points' coordinates
*/
void initElementalFieldInterpolation(
const ElementTypeMapArray<Real> & interpolation_points_coordinates);
/* ------------------------------------------------------------------------ */
/* Common part */
/* ------------------------------------------------------------------------ */
protected:
/* ------------------------------------------------------------------------ */
static inline UInt getTangentStiffnessVoigtSize(UInt dim);
/// compute the potential energy by element
void computePotentialEnergyByElements();
/* ------------------------------------------------------------------------ */
/* Finite deformation functions */
/* This functions area implementing what is described in the paper of Bathe */
/* et al, in IJNME, Finite Element Formulations For Large Deformation */
/* Dynamic Analysis, Vol 9, 353-386, 1975 */
/* ------------------------------------------------------------------------ */
protected:
/// assemble the residual
template <UInt dim> void assembleInternalForces(GhostType ghost_type);
template <UInt dim>
void computeAllStressesFromTangentModuli(ElementType type,
GhostType ghost_type);
template <UInt dim>
void assembleStiffnessMatrix(ElementType type, GhostType ghost_type);
/// assembling in finite deformation
template <UInt dim>
void assembleStiffnessMatrixNL(ElementType type, GhostType ghost_type);
template <UInt dim>
void assembleStiffnessMatrixL2(ElementType type, GhostType ghost_type);
/* ------------------------------------------------------------------------ */
/* Conversion functions */
/* ------------------------------------------------------------------------ */
public:
/// Size of the Stress matrix for the case of finite deformation see: Bathe et
/// al, IJNME, Vol 9, 353-386, 1975
static inline UInt getCauchyStressMatrixSize(UInt dim);
/// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386,
/// 1975
template <UInt dim>
static inline void setCauchyStressMatrix(const Matrix<Real> & S_t,
Matrix<Real> & sigma);
/// write the stress tensor in the Voigt notation.
template <UInt dim>
static inline decltype(auto) stressToVoigt(const Matrix<Real> & stress) {
return VoigtHelper<dim>::matrixToVoigt(stress);
}
/// write the strain tensor in the Voigt notation.
template <UInt dim>
static inline decltype(auto) strainToVoigt(const Matrix<Real> & strain) {
return VoigtHelper<dim>::matrixToVoigtWithFactors(strain);
}
/// write a voigt vector to stress
template <UInt dim>
static inline void voigtToStress(const Vector<Real> & voigt,
Matrix<Real> & stress) {
VoigtHelper<dim>::voigtToMatrix(voigt, stress);
}
/// Computation of Cauchy stress tensor in the case of finite deformation from
/// the 2nd Piola-Kirchhoff for a given element type
template <UInt dim>
void StoCauchy(ElementType el_type, GhostType ghost_type = _not_ghost);
/// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation
/// gradient
template <UInt dim>
inline void StoCauchy(const Matrix<Real> & F, const Matrix<Real> & S,
Matrix<Real> & sigma, const Real & C33 = 1.0) const;
template <UInt dim>
static inline void gradUToF(const Matrix<Real> & grad_u, Matrix<Real> & F);
template <UInt dim>
static inline decltype(auto) gradUToF(const Matrix<Real> & grad_u);
static inline void rightCauchy(const Matrix<Real> & F, Matrix<Real> & C);
static inline void leftCauchy(const Matrix<Real> & F, Matrix<Real> & B);
template <UInt dim>
static inline void gradUToEpsilon(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon);
template <UInt dim>
static inline decltype(auto) gradUToEpsilon(const Matrix<Real> & grad_u);
template <UInt dim>
static inline void gradUToE(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon);
template <UInt dim>
static inline decltype(auto) gradUToE(const Matrix<Real> & grad_u);
static inline Real stressToVonMises(const Matrix<Real> & stress);
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
inline UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void unpackData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) override;
/* ------------------------------------------------------------------------ */
/* SolidMechanicsModelEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
virtual void beforeSolveStep();
virtual void afterSolveStep(bool converged = true);
void onDamageIteration() override;
void onDamageUpdate() override;
void onDump() override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &)
AKANTU_GET_MACRO(Rho, rho, Real);
AKANTU_SET_MACRO(Rho, rho, Real);
- AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
-
/// return the potential energy for the subset of elements contained by the
/// material
Real getPotentialEnergy();
/// return the potential energy for the provided element
Real getPotentialEnergy(ElementType & type, UInt index);
/// return the energy (identified by id) for the subset of elements contained
/// by the material
virtual Real getEnergy(const std::string & type);
/// return the energy (identified by id) for the provided element
virtual Real getEnergy(const std::string & energy_id, ElementType type,
UInt index);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy,
Real);
AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(ElementFilter, element_filter,
const ElementTypeMapArray<UInt> &);
AKANTU_GET_MACRO(FEEngine, fem, FEEngine &);
bool isNonLocal() const { return is_non_local; }
bool isFiniteDeformation() const { return finite_deformation; }
bool isInelasticDeformation() const { return inelastic_deformation; }
/// apply a constant eigengrad_u everywhere in the material
virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u,
GhostType /*ghost_type*/ = _not_ghost);
bool hasMatrixChanged(const ID & id) {
if (id == "K") {
return hasStiffnessMatrixChanged() or finite_deformation;
}
return true;
}
MatrixType getMatrixType(const ID & id) {
if (id == "K") {
return getTangentType();
}
if (id == "M") {
return _symmetric;
}
return _mt_not_defined;
}
/// specify if the matrix need to be recomputed for this material
virtual bool hasStiffnessMatrixChanged() { return true; }
/// specify the type of matrix, if not overloaded the material is not valid
/// for static or implicit computations
virtual MatrixType getTangentType() { return _mt_not_defined; }
/// static method to reteive the material factory
static MaterialFactory & getFactory();
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// Link to the fem object in the model
FEEngine & fem;
/// Finite deformation
bool finite_deformation{false};
/// Finite deformation
bool inelastic_deformation{false};
/// The model to witch the material belong
SolidMechanicsModel & model;
/// density : rho
Real rho{0.};
/// spatial dimension
UInt spatial_dimension;
/// stresses arrays ordered by element types
InternalField<Real> stress;
/// eigengrad_u arrays ordered by element types
InternalField<Real> eigengradu;
/// grad_u arrays ordered by element types
InternalField<Real> gradu;
/// Green Lagrange strain (Finite deformation)
InternalField<Real> green_strain;
/// Second Piola-Kirchhoff stress tensor arrays ordered by element types
/// (Finite deformation)
InternalField<Real> piola_kirchhoff_2;
/// potential energy by element
InternalField<Real> potential_energy;
/// tell if using in non local mode or not
bool is_non_local{false};
/// tell if the material need the previous stress state
bool use_previous_stress{false};
/// tell if the material need the previous strain state
bool use_previous_gradu{false};
/// elemental field interpolation coordinates
InternalField<Real> interpolation_inverse_coordinates;
/// elemental field interpolation points
InternalField<Real> interpolation_points_matrices;
private:
/// eigen_grad_u for the parser
Matrix<Real> eigen_grad_u;
};
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream,
const Material & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#include "material_inline_impl.hh"
#include "internal_field_tmpl.hh"
#include "random_internal_field_tmpl.hh"
/* -------------------------------------------------------------------------- */
/* Auto loop */
/* -------------------------------------------------------------------------- */
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeStress. This macro in addition to write the loop
/// provides two tensors (matrices) sigma and grad_u
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \
auto && grad_u_view = \
make_view(this->gradu(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
auto stress_view = \
make_view(this->stress(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
if (this->isFiniteDeformation()) { \
stress_view = make_view(this->piola_kirchhoff_2(el_type, ghost_type), \
this->spatial_dimension, this->spatial_dimension); \
} \
\
for (auto && data : zip(grad_u_view, stress_view)) { \
[[gnu::unused]] Matrix<Real> & grad_u = std::get<0>(data); \
[[gnu::unused]] Matrix<Real> & sigma = std::get<1>(data)
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END }
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeTangentModuli. This macro in addition to write the
/// loop provides two tensors (matrices) sigma_tensor, grad_u, and a matrix
/// where the elemental tangent moduli should be stored in Voigt Notation
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \
auto && grad_u_view = \
make_view(this->gradu(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
auto && stress_view = \
make_view(this->stress(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
auto tangent_size = \
this->getTangentStiffnessVoigtSize(this->spatial_dimension); \
\
auto && tangent_view = make_view(tangent_mat, tangent_size, tangent_size); \
\
for (auto && data : zip(grad_u_view, stress_view, tangent_view)) { \
[[gnu::unused]] Matrix<Real> & grad_u = std::get<0>(data); \
[[gnu::unused]] Matrix<Real> & sigma = std::get<1>(data); \
Matrix<Real> & tangent = std::get<2>(data);
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END }
/* -------------------------------------------------------------------------- */
#define INSTANTIATE_MATERIAL_ONLY(mat_name) \
template class mat_name<1>; /* NOLINT */ \
template class mat_name<2>; /* NOLINT */ \
template class mat_name<3> /* NOLINT */
#define MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name) \
[](UInt dim, const ID &, SolidMechanicsModel & model, \
const ID & id) /* NOLINT */ \
-> std::unique_ptr< \
Material> { /* NOLINT */ \
switch (dim) { \
case 1: \
return std::make_unique<mat_name<1>>(/* NOLINT */ \
model, id); \
case 2: \
return std::make_unique<mat_name<2>>(/* NOLINT */ \
model, id); \
case 3: \
return std::make_unique<mat_name<3>>(/* NOLINT */ \
model, id); \
default: \
AKANTU_EXCEPTION( \
"The dimension " \
<< dim \
<< "is not a valid dimension for the material " \
<< #id); \
} \
}
#define INSTANTIATE_MATERIAL(id, mat_name) \
INSTANTIATE_MATERIAL_ONLY(mat_name); \
static bool material_is_alocated_##id [[gnu::unused]] = \
MaterialFactory::getInstance().registerAllocator( \
#id, MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name))
#endif /* AKANTU_MATERIAL_HH_ */