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 * @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 "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 inline void registerInternal(std::shared_ptr> & 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 const InternalField & getInternal(const ID & id) const; template InternalField & getInternal(const ID & id); template inline bool isInternal(const ID & id, const ElementKind & element_kind) const; 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); 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 & getElementFilter() const = 0; AKANTU_GET_MACRO(Name, name, const std::string &); AKANTU_GET_MACRO(ID, id, const ID &); + + private: std::map> internal_vectors; protected: ID id; + // spatial dimension for constitutive law + UInt spatial_dimension; + /// constitutive law name std::string name; }; template class ConstitutiveLaw : public ConstitutiveLawInternalHandler, public DataAccessor, 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 & 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: /// 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 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: /* ------------------------------------------------------------------------ */ virtual void onNodesAdded(const Array & /*unused*/, const NewNodesEvent & /*unused*/) override{}; virtual void onNodesRemoved(const Array & /*unused*/, const Array & /*unused*/, const RemovedNodesEvent & /*unused*/) override{}; virtual void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; virtual void onElementsAdded(const Array & /*unused*/, const NewElementsEvent & /*unused*/); virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); public: 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_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 & getElementFilter() const override {return element_filter;} + 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}; protected: + /// 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_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 * @author Nicolas Richart * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "constitutive_law.hh" #include "constitutive_law_non_local_interface.hh" #include "non_local_neighborhood.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -template -ConstitutiveLawNonLocal::ConstitutiveLawNonLocal( +template +ConstitutiveLawNonLocal::ConstitutiveLawNonLocal( typename ConstitutiveLawParent::ConstitutiveLawsHandler & handler, const ID & id) : ConstitutiveLawParent(handler, id) { } /* -------------------------------------------------------------------------- */ -template -void ConstitutiveLawNonLocal:: +template +void ConstitutiveLawNonLocal:: 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 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 -void ConstitutiveLawNonLocal::updateNonLocalInternals( +template +void ConstitutiveLawNonLocal::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 & internal = this->template getInternal(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 -void ConstitutiveLawNonLocal::registerNeighborhood() { +template +void ConstitutiveLawNonLocal::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 * @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" #include "fe_engine.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_CONSTITUTIVE_LAW_TMPL_HH #define AKANTU_CONSTITUTIVE_LAW_TMPL_HH namespace akantu { /* -------------------------------------------------------------------------- */ template inline void ConstitutiveLawInternalHandler::registerInternal( std::shared_ptr> & 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 const InternalField & ConstitutiveLawInternalHandler::getInternal(const ID & id) const { auto it = internal_vectors.find(this->id + ":" + id); if (it != internal_vectors.end() and aka::is_of_type>(*it->second)) { return aka::as_type>(*it->second); } AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << id << " (" << (getID() + ":" + id) << ")"); } /* -------------------------------------------------------------------------- */ template InternalField & ConstitutiveLawInternalHandler::getInternal(const ID & id) { auto it = internal_vectors.find(getID() + ":" + id); if (it != internal_vectors.end() and aka::is_of_type>(*it->second)) { return aka::as_type>(*it->second); } AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << id << " (" << (getID() + ":" + id) << ")"); } /* -------------------------------------------------------------------------- */ template 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>(*it->second) and aka::as_type>(*it->second).getElementKind() != element_kind); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template ConstitutiveLaw::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 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::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 && pair : internal_vectors) { pair.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 pair : internal_vectors) { pair.second->removeIntegrationPoints(constitutive_law_local_new_numbering); } } /* -------------------------------------------------------------------------- */ template template inline void ConstitutiveLaw::packElementDataHelper( const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) const { DataAccessor::packElementalDataHelper(data_to_pack, buffer, elements, true, handler.getFEEngine(fem_id)); } /* -------------------------------------------------------------------------- */ template template inline void ConstitutiveLaw::unpackElementDataHelper( ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) { DataAccessor::unpackElementalDataHelper(data_to_unpack, buffer, elements, true, handler.getFEEngine(fem_id)); } /* -------------------------------------------------------------------------- */ template inline Element ConstitutiveLaw::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 inline Element ConstitutiveLaw::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 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 ElementTypeMap ConstitutiveLaw::getInternalDataPerElem( const ID & field_id, ElementKind element_kind) const { if (not this->template isInternal(field_id, element_kind)) { AKANTU_EXCEPTION("Cannot find internal field " << id << " in the constitutive law " << this->name); } const auto & 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 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/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 * @author Mohit Pundir * * @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 . * */ /* -------------------------------------------------------------------------- */ #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 ConstitutiveLawsHandler : public DataAccessor, 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( 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 &); AKANTU_GET_MACRO(ConstitutiveLawList, constitutive_law, const Array &); protected: Array 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 void for_each_constitutive_law(Func && func) { for (auto && constitutive_law : constitutive_laws) { std::forward(func)( aka::as_type(*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 * filter = nullptr); protected: void splitElementByConstitutiveLaw( const Array & elements, std::vector> & elements_per_cl) const; template void splitByConstitutiveLaw(const Array & 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 getInternalDataPerElem(const std::string & field_name, ElementKind kind); //! flatten a given constitutive_law internal field ElementTypeMapArray & 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 &); AKANTU_GET_MACRO(ConstitutiveLawLocalNumbering, constitutive_law_local_numbering, const ElementTypeMapArray &); /// 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 constitutive_law_selector) { this->constitutive_law_selector = std::move(constitutive_law_selector); } /* ------------------------------------------------------------------------ */ /* 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; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & 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 constitutive_laws_names_to_id; /// Arrays containing the constitutive_law index for each element ElementTypeMapArray constitutive_law_index; /// Arrays containing the position in the element filter of the /// constitutive_law (constitutive_law's local numbering) ElementTypeMapArray constitutive_law_local_numbering; /// list of used constitutive_laws std::vector> constitutive_laws; /// class defining of to choose a constitutive_law std::shared_ptr constitutive_law_selector; using flatten_internal_map = std::map, std::unique_ptr>>; /// 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 non_local_manager; template 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 * @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 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 . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "mesh_iterators.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) - : ConstitutiveLaw(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(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(model, id, _ek_regular),fem(fe_engine), - model(model), spatial_dimension(dim), + : ConstitutiveLaw(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 &>(model.getInternalForce()); // Mesh & mesh = fem.getMesh(); for (auto && type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by /// @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ auto * sigma_dphi_dx = new Array(nb_element * nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); fem.computeBtD(stress(type, ghost_type), *sigma_dphi_dx, type, ghost_type, elem_filter); /** * 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(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 & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { continue; } Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, spatial_dimension, type, ghost_type, elem_filter); gradu_vect -= eigengradu(type, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllCauchyStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(finite_deformation, "The Cauchy stress can only be " "computed if you are working in " "finite deformation."); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { switch (spatial_dimension) { case 1: this->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 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 & grad_u = *gradu_it; Matrix & piola = *piola_it; Matrix & sigma = *stress_it; auto F_tensor = gradUToF(grad_u); this->StoCauchy(F_tensor, piola, sigma); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & displacement = model.getDisplacement(); // resizeInternalArray(gradu); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(displacement, gradu_vect, spatial_dimension, type, ghost_type, elem_filter); setToSteadyState(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { if (finite_deformation) { switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL<1>(type, ghost_type); assembleStiffnessMatrixL2<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrixNL<2>(type, ghost_type); assembleStiffnessMatrixL2<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrixNL<3>(type, ghost_type); assembleStiffnessMatrixL2<3>(type, ghost_type); break; } } } else { switch (spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(type, ghost_type); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { AKANTU_DEBUG_OUT(); return; } // const Array & shapes_derivatives = // fem.getShapesDerivatives(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); auto * tangent_stiffness_matrix = new Array(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(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(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); // Array & gradu_vect = delta_gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto * shapes_derivatives_filtered = new Array( 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(nb_element * nb_quadrature_points, bt_s_b_size * bt_s_b_size, "B^t*D*B"); UInt piola_matrix_size = getCauchyStressMatrixSize(dim); Matrix B(piola_matrix_size, bt_s_b_size); Matrix Bt_S(bt_s_b_size, piola_matrix_size); Matrix S(piola_matrix_size, piola_matrix_size); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); auto Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); auto piola_it = piola_kirchhoff_2(type, ghost_type).begin(dim, dim); for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) { auto & Bt_S_B = *Bt_S_B_it; const auto & Piola_kirchhoff_matrix = *piola_it; setCauchyStressMatrix(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.template mul(B, S); Bt_S_B.template mul(Bt_S, B); } delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto * K_e = new Array(nb_element, bt_s_b_size * bt_s_b_size, "K_e"); fem.integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, elem_filter); delete bt_s_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); auto * tangent_stiffness_matrix = new Array(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( 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(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_D_B_it = bt_d_b->begin(bt_d_b_size, bt_d_b_size); auto grad_u_it = gradu_vect.begin(dim, dim); auto D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); auto D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) { const auto & grad_u = *grad_u_it; const auto & D = *D_it; auto & Bt_D_B = *Bt_D_B_it; // transferBMatrixToBL1 (*shapes_derivatives_filtered_it, B, // nb_nodes_per_element); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element); B += B2; Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); Mesh & mesh = fem.getMesh(); for (auto type : element_filter.elementTypes(_ghost_type = ghost_type)) { const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & 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( nb_element, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); // Set stress vectors UInt stress_size = getTangentStiffnessVoigtSize(dim); // Set matrices B and BNL* UInt bt_s_size = dim * nb_nodes_per_element; auto * bt_s = new Array(nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); auto grad_u_it = this->gradu(type, ghost_type).begin(dim, dim); auto grad_u_end = this->gradu(type, ghost_type).end(dim, dim); auto stress_it = this->piola_kirchhoff_2(type, ghost_type).begin(dim, dim); shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1); Matrix B_tensor(stress_size, bt_s_size); Matrix B2_tensor(stress_size, bt_s_size); for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) { auto & grad_u = *grad_u_it; auto & r = *bt_s_it; auto & S = *stress_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element); B_tensor += B2_tensor; auto S_vect = Material::stressToVoigt(S); Matrix S_voigt(S_vect.storage(), stress_size, 1); r.template mul(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(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 epot_on_quad_points(fem.getNbIntegrationPoints(type)); computePotentialEnergyByElement(type, index, epot_on_quad_points); epot = fem.integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(const std::string & type) { AKANTU_DEBUG_IN(); if (type == "potential") { return getPotentialEnergy(); } AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(const std::string & energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "potential") { return getPotentialEnergy(type, index); } AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation( const ElementTypeMapArray & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); this->fem.initElementalFieldInterpolationFromIntegrationPoints( interpolation_points_coordinates, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, &(this->element_filter)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type) { this->fem.interpolateElementalFieldFromIntegrationPoints( this->stress, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, result, ghost_type, &(this->element_filter)); } /* -------------------------------------------------------------------------- */ void Material::interpolateStressOnFacets( ElementTypeMapArray & result, ElementTypeMapArray & by_elem_result, const GhostType ghost_type) { interpolateStress(by_elem_result, ghost_type); UInt stress_size = this->stress.getNbComponent(); const Mesh & mesh = this->model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_fil = element_filter(type, ghost_type); Array & by_elem_res = by_elem_result(type, ghost_type); UInt nb_element = elem_fil.size(); UInt nb_element_full = this->model.getMesh().getNbElement(type, ghost_type); UInt nb_interpolation_points_per_elem = by_elem_res.size() / nb_element_full; const Array & facet_to_element = mesh_facets.getSubelementToElement(type, ghost_type); ElementType type_facet = Mesh::getFacetType(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); UInt nb_quad_per_facet = nb_interpolation_points_per_elem / nb_facet_per_elem; Element element_for_comparison{type, 0, ghost_type}; const Array> * element_to_facet = nullptr; GhostType current_ghost_type = _casper; Array * result_vec = nullptr; Array::const_matrix_iterator result_it = by_elem_res.begin_reinterpret( stress_size, nb_interpolation_points_per_elem, nb_element_full); for (UInt el = 0; el < nb_element; ++el) { UInt global_el = elem_fil(el); element_for_comparison.element = global_el; for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element facet_elem = facet_to_element(global_el, f); UInt global_facet = facet_elem.element; if (facet_elem.ghost_type != current_ghost_type) { current_ghost_type = facet_elem.ghost_type; element_to_facet = &mesh_facets.getElementToSubelement( type_facet, current_ghost_type); result_vec = &result(type_facet, current_ghost_type); } bool is_second_element = (*element_to_facet)(global_facet)[0] != element_for_comparison; for (UInt q = 0; q < nb_quad_per_facet; ++q) { Vector result_local(result_vec->storage() + (global_facet * nb_quad_per_facet + q) * result_vec->getNbComponent() + static_cast(is_second_element) * stress_size, stress_size); const Matrix & result_tmp(result_it[global_el]); result_local = result_tmp(f * nb_quad_per_facet + q); } } } } } /* -------------------------------------------------------------------------- */ template const Array & Material::getArray(const ID & vect_id, ElementType type, GhostType ghost_type) const { try { return this->template getInternal(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 Array & Material::getArray(const ID & vect_id, ElementType type, GhostType ghost_type) { try { return this->template getInternal(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 & Material::getArray(const ID & vect_id, ElementType type, GhostType ghost_type) const; /* -------------------------------------------------------------------------- */ template Array & Material::getArray(const ID & vect_id, ElementType type, GhostType ghost_type); /* -------------------------------------------------------------------------- */ template const Array & Material::getArray(const ID & vect_id, ElementType type, GhostType ghost_type) const; /* -------------------------------------------------------------------------- */ template Array & 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 & /*unused*/, const NewElementsEvent & /*unused*/) { this->resizeInternals(); } /* -------------------------------------------------------------------------- */ void Material::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, [[gnu::unused]] const RemovedElementsEvent & event) { UInt my_num = model.getInternalIndexFromID(getID()); ElementTypeMapArray 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 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 & point, Matrix & extrapolated) { if (this->isInternal(id, element.kind())) { UInt nb_element = this->element_filter(element.type, element.ghost_type).size(); const ID name = this->getID() + ":" + id; UInt nb_quads = this->internal_vectors_real[name]->getFEEngine().getNbIntegrationPoints( element.type, element.ghost_type); const Array & internal = this->getArray(id, element.type, element.ghost_type); UInt nb_component = internal.getNbComponent(); Array::const_matrix_iterator internal_it = internal.begin_reinterpret(nb_component, nb_quads, nb_element); Element local_element = this->convertToLocalElement(element); /// instead of really extrapolating, here the value of the first GP /// is copied into the result vector. This works only for linear /// elements /// @todo extrapolate!!!! AKANTU_DEBUG_WARNING("This is a fix, values are not truly extrapolated"); const Matrix & values = internal_it[local_element.element]; UInt index = 0; Vector tmp(nb_component); for (UInt j = 0; j < values.cols(); ++j) { tmp = values(j); if (tmp.norm() > 0) { index = j; break; } } for (UInt i = 0; i < extrapolated.size(); ++i) { extrapolated(i) = values(index); } } else { Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); extrapolated = default_values; } } /* -------------------------------------------------------------------------- */ void Material::applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const GhostType ghost_type) { for (auto && type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { if (element_filter(type, ghost_type).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 * @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, 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; /* ------------------------------------------------------------------------ */ /* 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; } 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_ */