diff --git a/src/model/common/constitutive_laws/constitutive_law.hh b/src/model/common/constitutive_laws/constitutive_law.hh index fef4dde4e..61fb1ea44 100644 --- a/src/model/common/constitutive_laws/constitutive_law.hh +++ b/src/model/common/constitutive_laws/constitutive_law.hh @@ -1,182 +1,196 @@ /** * @file constitutive_law.hh * * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "data_accessor.hh" +#include "mesh_events.hh" #include "parsable.hh" #include "parser.hh" -#include "mesh_events.hh" /* -------------------------------------------------------------------------- */ #include "internal_field.hh" #include "random_internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONSTITUTIVE_LAW_HH__ #define __AKANTU_CONSTITUTIVE_LAW_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { template class ConstitutiveLaw : public DataAccessor, 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); /// Destructor ~ConstitutiveLaw() override; protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - template void registerInternal(InternalField & /*vect*/) { + template + void registerInternal(InternalField & /*vect*/) { AKANTU_TO_IMPLEMENT(); } - template void unregisterInternal(InternalField & /*vect*/) { + template + void unregisterInternal(InternalField & /*vect*/) { AKANTU_TO_IMPLEMENT(); } /// initialize the constitutive law computed parameter virtual void initConstitutiveLaw(); /// save the internals in the previous_state if needed virtual void savePreviousState(); /// restore the internals from previous_state if needed virtual void restorePreviousState(); /// 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 & elements_to_add); /// remove many element at once void removeElements(const Array & elements_to_remove); /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; protected: /// resize the internals arrrays virtual void resizeInternals(); /// function called to update the internal parameters when the /// modifiable parameters are modified virtual void updateInternalParameters() {} public: virtual void onElementsAdded(const Array & /*unused*/, const NewElementsEvent & /*unused*/); - virtual void onElementsRemoved(const Array & element_list, - const ElementTypeMapArray & new_numbering, - const RemovedElementsEvent & event); + virtual void + onElementsRemoved(const Array & element_list, + const ElementTypeMapArray & new_numbering, + const RemovedElementsEvent & event); public: - template - const InternalField & getInternal(const ID & id) const; - template InternalField & getInternal(const ID & id); + template + const InternalField & getInternal(const ID & /*id*/) const { + AKANTU_TO_IMPLEMENT(); + return NULL; + } - template - inline bool isInternal(const ID & id, const ElementKind & element_kind) const; + template + InternalField & getInternal(const ID & /*id*/) { + AKANTU_TO_IMPLEMENT(); + return NULL; + } + + template + inline bool isInternal(const ID & /*id*/, + const ElementKind & /*element_kind*/) const { + return false; + } template inline void setParam(const ID & param, T value); inline const Parameter & getParam(const ID & param) const; template void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Name, name, const std::string &); AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); template ElementTypeMap 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}; std::map *> internal_vectors_real; std::map *> internal_vectors_uint; std::map *> internal_vectors_bool; protected: ID id; /// constitutive law name std::string name; /// list of element handled by the constitutive law ElementTypeMapArray 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_tmpl.hh b/src/model/common/constitutive_laws/constitutive_law_tmpl.hh index aa1593808..5ab6e034f 100644 --- a/src/model/common/constitutive_laws/constitutive_law_tmpl.hh +++ b/src/model/common/constitutive_laws/constitutive_law_tmpl.hh @@ -1,565 +1,607 @@ /** * @file constitutive_law_tmpl.hh * * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: 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 . * */ /* -------------------------------------------------------------------------- */ #include "constitutive_law.hh" #include "constitutive_laws_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_CONSTITUTIVE_LAW_TMPL_HH #define AKANTU_CONSTITUTIVE_LAW_TMPL_HH namespace akantu { template ConstitutiveLaw::ConstitutiveLaw( ConstitutiveLawsHandler_ & handler, const ID & id, ElementKind element_kind) : Parsable(ParserType::_constitutive_law, id), handler(handler), id(id), 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 void ConstitutiveLaw::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); } /* -------------------------------------------------------------------------- */ template void ConstitutiveLaw::initConstitutiveLaw() { this->resizeInternals(); this->updateInternalParmaters(); is_init = true; } /* -------------------------------------------------------------------------- */ template void ConstitutiveLaw::savePreviousState() { for (auto pair : internal_vectors_real) { if (pair.second->hasHistory()) { pair.second->saveCurrentValues(); } } } /* -------------------------------------------------------------------------- */ template void ConstitutiveLaw::restorePreviousState() { for (auto pair : internal_vectors_real) { if (pair.second->hasHistory()) { pair.second->restorePreviousValues(); } } } /* -------------------------------------------------------------------------- */ template void ConstitutiveLaw::resizeInternals() { 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(); } } /* -------------------------------------------------------------------------- */ template void ConstitutiveLaw::addElements( const Array & 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 void ConstitutiveLaw::removeElements( const Array & 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 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 element_filter_tmp("element_filter_tmp", id); element_filter_tmp.initialize(mesh, _element_filter = &element_filter, _element_kind = _ek_not_defined); ElementTypeMap 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 it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { it->second->removeIntegrationPoints(constitutive_law_local_new_numbering); } for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { it->second->removeIntegrationPoints(constitutive_law_local_new_numbering); } for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) { it->second->removeIntegrationPoints(constitutive_law_local_new_numbering); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ConstitutiveLaw::onElementsAdded( const Array & /*unused*/, const NewElementsEvent & /*unused*/) { this->resizeInternals(); } /* -------------------------------------------------------------------------- */ template void ConstitutiveLaw::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & /*event*/) { UInt my_num = handler.getInternalIndexFromID(getID()); ElementTypeMapArray 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 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 it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { it->second->removeIntegrationPoints(constitutive_law_local_new_numbering); } for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { it->second->removeIntegrationPoints(constitutive_law_local_new_numbering); } for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) { it->second->removeIntegrationPoints(constitutive_law_local_new_numbering); } } /* -------------------------------------------------------------------------- */ template inline UInt ConstitutiveLaw::addElement( ElementType type, UInt element, GhostType ghost_type) { Array & el_filter = this->element_filter(type, ghost_type); el_filter.push_back(element); return el_filter.size() - 1; } /* -------------------------------------------------------------------------- */ template inline UInt ConstitutiveLaw::addElement(const Element & element) { return this->addElement(element.type, element.element, element.ghost_type); } /* -------------------------------------------------------------------------- */ template inline const Parameter & ConstitutiveLaw::getParam(const ID & param) const { try { return get(param); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the constitutive law" << getID()); } } /* -------------------------------------------------------------------------- */ template template inline void ConstitutiveLaw::setParam(const ID & param, T value) { try { set(param, value); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the constitutive law " << getID()); } updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template -template <> -inline void ConstitutiveLaw::registerInternal( +template +inline void +ConstitutiveLaw::registerInternal( InternalField & vect) { internal_vectors_real[vect.getID()] = &vect; } template -template <> -inline void ConstitutiveLaw::registerInternal( +template +inline void +ConstitutiveLaw::registerInternal( InternalField & vect) { internal_vectors_uint[vect.getID()] = &vect; } template -template <> -inline void ConstitutiveLaw::registerInternal( +template +inline void +ConstitutiveLaw::registerInternal( InternalField & vect) { internal_vectors_bool[vect.getID()] = &vect; } /* -------------------------------------------------------------------------- */ template -template <> -inline void ConstitutiveLaw::unregisterInternal( +template +inline void +ConstitutiveLaw::unregisterInternal( InternalField & vect) { internal_vectors_real.erase(vect.getID()); } template -template <> -inline void ConstitutiveLaw::unregisterInternal( +template +inline void +ConstitutiveLaw::unregisterInternal( InternalField & vect) { internal_vectors_uint.erase(vect.getID()); } template -template <> -inline void ConstitutiveLaw::unregisterInternal( +template +inline void +ConstitutiveLaw::unregisterInternal( InternalField & vect) { internal_vectors_bool.erase(vect.getID()); } /* -------------------------------------------------------------------------- */ template -template -const InternalField & ConstituitveLaw::getInternal( - const ID & /*int_id*/) const { - AKANTU_TO_IMPLEMENT(); - return NULL; -} - -/* -------------------------------------------------------------------------- */ -template -InternalField & -ConstituitveLaw::getInternal(const ID & /*int_id*/) { - AKANTU_TO_IMPLEMENT(); - return NULL; -} - -/* -------------------------------------------------------------------------- */ -template <> +template const InternalField & -ConstituitveLaw::getInternal( +ConstituitveLaw::getInternal( const ID & int_id) const { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ +template template <> -InternalField & -ConstituitveLaw::getInternal(const ID & int_id) { +InternalField & +ConstituitveLaw::getInternal( + const ID & int_id) { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ -template <> +template +template const InternalField & -ConstituitveLaw::getInternal( +ConstituitveLaw::getInternal( const ID & int_id) const { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ -template <> +template +template InternalField & -ConstituitveLaw::getInternal(const ID & int_id) { +ConstituitveLaw::getInternal( + const ID & int_id) { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template -inline bool ConstitutiveLaw::isInternal( - const ID & /*id*/, ElementKind /*element_kind*/) const { - AKANTU_TO_IMPLEMENT(); +template +const InternalField & +ConstituitveLaw::getInternal( + const ID & int_id) const { + auto it = internal_vectors_bool.find(getID() + ":" + int_id); + if (it == internal_vectors_bool.end()) { + AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() + << ") does not contain an internal " + << int_id << " (" + << (getID() + ":" + int_id) << ")"); + } + return *it->second; +} + +/* -------------------------------------------------------------------------- */ +template +template +InternalField & +ConstituitveLaw::getInternal( + const ID & int_id) { + auto it = internal_vectors_bool.find(getID() + ":" + int_id); + if (it == internal_vectors_bool.end()) { + AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() + << ") does not contain an internal " + << int_id << " (" + << (getID() + ":" + int_id) << ")"); + } + return *it->second; } /* -------------------------------------------------------------------------- */ template -template <> -inline bool ConstitutiveLaw::isInternal( +template +inline bool ConstitutiveLaw::isInternal( const ID & id, ElementKind element_kind) const { auto internal_array = internal_vectors_real.find(this->getID() + ":" + id); return not(internal_array == internal_vectors_real.end() || internal_array->second->getElementKind() != element_kind); } +/* -------------------------------------------------------------------------- */ +template +template +inline bool ConstitutiveLaw::isInternal( + const ID & id, ElementKind element_kind) const { + auto internal_array = internal_vectors_uint.find(this->getID() + ":" + id); + + return not(internal_array == internal_vectors_uint.end() || + internal_array->second->getElementKind() != element_kind); +} + +/* -------------------------------------------------------------------------- */ +template +template +inline bool ConstitutiveLaw::isInternal( + const ID & id, ElementKind element_kind) const { + auto internal_array = internal_vectors_bool.find(this->getID() + ":" + id); + + return not(internal_array == internal_vectors_bool.end() || + internal_array->second->getElementKind() != element_kind); +} + /* -------------------------------------------------------------------------- */ template inline ElementTypeMap ConstitutiveLaw::getInternalDataPerElem( const ID & field_id, ElementKind element_kind) const { if (!this->template isInternal(field_id, element_kind)) { AKANTU_EXCEPTION("Cannot find internal field " << id << " in the constitutive law " << this->name); } const InternalField & internal_field = this->template getInternal(field_id); const FEEngine & fe_engine = internal_field.getFEEngine(); UInt nb_data_per_quad = internal_field.getNbComponent(); ElementTypeMap res; for (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 void ConstitutiveLaw::flattenInternal( const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type, ElementKind element_kind) const { if (!this->template isInternal(field_id, element_kind)) { AKANTU_EXCEPTION("Cannot find internal field " << id << " in the constitutive law " << this->name); } const InternalField & internal_field = this->template getInternal(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 & src_vect = internal_field(type, ghost_type); const Array & filter = internal_field.getFilter(type, ghost_type); // total number of elements in the corresponding mesh UInt nb_element_dst = mesh.getNbElement(type, ghost_type); // number of element in the internal field UInt nb_element_src = filter.size(); // number of quadrature points per elem UInt nb_quad_per_elem = fe_engine.getNbIntegrationPoints(type); // number of data per quadrature point UInt nb_data_per_quad = internal_field.getNbComponent(); if (!internal_flat.exists(type, ghost_type)) { internal_flat.alloc(nb_element_dst * nb_quad_per_elem, nb_data_per_quad, type, ghost_type); } if (nb_element_src == 0) { continue; } // number of data per element UInt nb_data = nb_quad_per_elem * nb_data_per_quad; Array & dst_vect = internal_flat(type, ghost_type); dst_vect.resize(nb_element_dst * nb_quad_per_elem); 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/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 4653e13df..339591dba 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,630 +1,621 @@ /** * @file material.hh * * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @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 . * */ /* -------------------------------------------------------------------------- */ #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; /** * 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 & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : public ConstitutiveLaw public DataAccessor, public Parsable, public MeshEventHandler, 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 & /*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 & /*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 & points, Matrix & 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 & 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 & result, ElementTypeMapArray & 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 & 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 void assembleInternalForces(GhostType ghost_type); template void computeAllStressesFromTangentModuli(ElementType type, GhostType ghost_type); template void assembleStiffnessMatrix(ElementType type, GhostType ghost_type); /// assembling in finite deformation template void assembleStiffnessMatrixNL(ElementType type, GhostType ghost_type); template 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 static inline void setCauchyStressMatrix(const Matrix & S_t, Matrix & sigma); /// write the stress tensor in the Voigt notation. template static inline decltype(auto) stressToVoigt(const Matrix & stress) { return VoigtHelper::matrixToVoigt(stress); } /// write the strain tensor in the Voigt notation. template static inline decltype(auto) strainToVoigt(const Matrix & strain) { return VoigtHelper::matrixToVoigtWithFactors(strain); } /// write a voigt vector to stress template static inline void voigtToStress(const Vector & voigt, Matrix & stress) { VoigtHelper::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 void StoCauchy(ElementType el_type, GhostType ghost_type = _not_ghost); /// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation /// gradient template inline void StoCauchy(const Matrix & F, const Matrix & S, Matrix & sigma, const Real & C33 = 1.0) const; template static inline void gradUToF(const Matrix & grad_u, Matrix & F); template static inline decltype(auto) gradUToF(const Matrix & grad_u); static inline void rightCauchy(const Matrix & F, Matrix & C); static inline void leftCauchy(const Matrix & F, Matrix & B); template static inline void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); template static inline decltype(auto) gradUToEpsilon(const Matrix & grad_u); template static inline void gradUToE(const Matrix & grad_u, Matrix & epsilon); template static inline decltype(auto) gradUToE(const Matrix & grad_u); static inline Real stressToVonMises(const Matrix & stress); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; template inline void packElementDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); /* ------------------------------------------------------------------------ */ /* MeshEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ void onNodesAdded(const Array & /*unused*/, const NewNodesEvent & /*unused*/) override{}; void onNodesRemoved(const Array & /*unused*/, const Array & /*unused*/, const RemovedNodesEvent & /*unused*/) override{}; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) 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 &); AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray &); AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); AKANTU_GET_MACRO(FEEngine, fem, FEEngine &); bool isNonLocal() const { return is_non_local; } - template - const Array & getArray(const ID & id, ElementType type, - GhostType ghost_type = _not_ghost) const; - template - Array & getArray(const ID & id, ElementType type, - GhostType ghost_type = _not_ghost); - - - 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 & 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 stress; /// eigengrad_u arrays ordered by element types InternalField eigengradu; /// grad_u arrays ordered by element types InternalField gradu; /// Green Lagrange strain (Finite deformation) InternalField green_strain; /// Second Piola-Kirchhoff stress tensor arrays ordered by element types /// (Finite deformation) InternalField piola_kirchhoff_2; /// potential energy by element InternalField 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 interpolation_inverse_coordinates; /// elemental field interpolation points InternalField interpolation_points_matrices; private: /// eigen_grad_u for the parser Matrix 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 & grad_u = std::get<0>(data); \ [[gnu::unused]] Matrix & 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 & grad_u = std::get<0>(data); \ [[gnu::unused]] Matrix & sigma = std::get<1>(data); \ Matrix & 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>(/* NOLINT */ \ model, id); \ case 2: \ return std::make_unique>(/* NOLINT */ \ model, id); \ case 3: \ return std::make_unique>(/* 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_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index ff8c82b57..f30c14649 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,495 +1,495 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * * 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 . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "constitutive_laws_handler.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_HH_ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public ConstitutiveLawsHandler, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "solid_mechanics_model", ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; public: /// initialize all internal arrays for materials virtual void initMaterials(); protected: /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(ElementType type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; void dump(const std::string & dumper_name) override; void dump(const std::string & dumper_name, UInt step) override; void dump(const std::string & dumper_name, Real time, UInt step) override; void dump() override; void dump(UInt step) override; void dump(Real time, UInt step) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement); /// get the SolidMechanicsModel::previous_displacement array AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement); /// get the SolidMechanicsModel::current_position array const Array & getCurrentPosition(); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment); /// get the lumped SolidMechanicsModel::mass array AKANTU_GET_MACRO_DEREF_PTR(Mass, mass); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force); /// get the SolidMechanicsModel::force array (external forces) [[deprecated("Use getExternalForce instead of this function")]] Array & getForce() { return getExternalForce(); } /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force); /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs); /// get an iterable on the materials - inline decltype(auto) getMaterials() { - this->getConstitutiveLaws(); - } + inline decltype(auto) getMaterials() { return this->getConstitutiveLaws(); } /// get an iterable on the materials inline decltype(auto) getMaterials() const { - this->getConstitutiveLaws(); + return this->getConstitutiveLaws(); } /// get a particular material (by numerical material index) inline Material & getMaterial(UInt mat_index) { - this->getConstitutiveLaw(mat_index); + return this->getConstitutiveLaw(mat_index); } /// get a particular material (by numerical material index) inline const Material & getMaterial(UInt mat_index) const { - this->getConstitutiveLaw(mat_index); + return this->getConstitutiveLaw(mat_index); } /// get a particular material (by material name) - inline Material & getMaterial(const std::string & name) { - this->getConstitutiveLaw(name); + inline Material & getMaterial(const std::string & name) { + return this->getConstitutiveLaw(name); } /// get a particular material (by material name) - inline const Material & getMaterial(const std::string & name) const { - this->getConstitutiveLaw(name); + inline const Material & getMaterial(const std::string & name) const { + return this->getConstitutiveLaw(name); } /// get a particular material id from is name - inline UInt getMaterialIndex(const std::string & name) const; + inline UInt getMaterialIndex(const std::string & name) const { + return this->getConstitutiveLawIndex(name); + } /// give the number of materials - inline UInt getNbMaterials() const { return getNbConstitutiveLaws(); } + inline UInt getNbMaterials() const { return this->getNbConstitutiveLaws(); } /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for an element Real getEnergy(const std::string & energy_id, ElementType type, UInt index); /// compute the energy for an element Real getEnergy(const std::string & energy_id, const Element & element) { return getEnergy(energy_id, element.type, element.element); } /// compute the energy for an element group Real getEnergy(const ID & energy_id, const ID & group_id); // this function is kept for backward compatinility decltype(auto) getMaterialByElement() const { return this->getConstitutiveLawByElement(); } // this function is kept for backward compatinility decltype(auto) getMaterialLocalNumbering() const { return this->getConstitutiveLawLocalNumbering(); } // this function is kept for backward compatinility decltype(auto) getMaterialByElement(ElementType type, GhostType ghost_type = _not_ghost) const { return this->getConstitutiveLawByElement(type, ghost_type); } // this function is kept for backward compatinility decltype(auto) getMaterialLocalNumbering(ElementType type, GhostType ghost_type = _not_ghost) const { return this->getConstitutiveLawLocalNumbering(type, ghost_type); } // this function is kept for backward compatinility decltype(auto) getMaterialSelector() { return this->getConstitutiveLawSelector(); } // this function is kept for backward compatinility void setMaterialSelector(std::shared_ptr material_selector) { this->setConstitutiveLawSelector(std::move(material_selector)); } /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: /// compute the stable time step Real getStableTimeStep(GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// release version of the displacement array UInt displacement_release{0}; /// release version of the current_position array UInt current_position_release{0}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; protected: /// conversion coefficient form force/mass to acceleration Real f_m2a{1.0}; /// displacements array std::unique_ptr> displacement; /// displacements array at the previous time step (used in finite deformation) std::unique_ptr> previous_displacement; /// increment of displacement std::unique_ptr> displacement_increment; /// lumped mass array std::unique_ptr> mass; /// velocities array std::unique_ptr> velocity; /// accelerations array std::unique_ptr> acceleration; /// external forces array std::unique_ptr> external_force; /// internal forces array std::unique_ptr> internal_force; /// array specifing if a degree of freedom is blocked or not std::unique_ptr> blocked_dofs; /// array of current position used during update residual std::unique_ptr> current_position; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */