constitutive_laws_tmpl.hh
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Nov 30, 02:10

constitutive_laws_tmpl.hh
View Options

	/**
	* @file constitutive_law_tmpl.hh *
	* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
	* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
	* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
	* @author Nicolas Richart <nicolas.richart@epfl.ch>
	* @author Marco Vocialta <marco.vocialta@epfl.ch>
	*
	* @date creation: Tue Jul 27 2010
	* @date last modification: Wed Feb 21 2018
	*
	* @brief Implementation of the templated part of the constitutive law class
	*
	*
	* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
	* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
	*
	* Akantu is free software: you can redistribute it and/or modify it under the
	* terms of the GNU Lesser General Public License as published by the Free
	* Software Foundation, either version 3 of the License, or (at your option) any
	* later version.
	*
	* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
	* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
	* details.
	*
	* You should have received a copy of the GNU Lesser General Public License
	* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
	*
	*/

	/* -------------------------------------------------------------------------- */
	#include "constitutive_laws_handler.hh"
	/* -------------------------------------------------------------------------- */

	#ifndef AKANTU_CONSTITUTIVE_TMPL_HH__
	#define AKANTU_CONSTITUTIVE_TMPL_HH__

	namespace akantu {
	template<class ConstitutiveLawsHandler_>
	ConstitutiveLaw<ConstitutiveLawsHandler_>::ConstitutiveLaw(ConstitutiveLawsHandler_ & handler,
	const ID & id,
	ElementKind element_kind)
	: Parsable(ParserType::_constitutive_law, id), handler(handler),
	id(id), element_filter("element_filter", id) {

	AKANTU_DEBUG_IN();

	/// 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();

	AKANU_DEBUG_OUT();

	}

	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::initialize(){
	registerParam("name", name, std::string(), _pat_parsable \| _pat_readable);
	}


	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::initConstitutiveLaw() {

	this->resizeInternals();

	updateInternalParmaters();

	is_init = true;

	}


	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::savePreviousState() {
	AKANTU_DEBUG_IN();

	for (auto pair : internal_vectors_real) {
	if (pair.second->hasHistory()) {
	pair.second->saveCurrentValues();
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::restorePreviousState() {
	AKANTU_DEBUG_IN();

	for (auto pair : internal_vectors_real) {
	if (pair.second->hasHistory()) {
	pair.second->restorePreviousValues();
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::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();
	}



	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::addElements(const Array<Element> & elements_to_add) {
	AKANTU_DEBUG_IN();

	UInt law_id = handler.getConstitutiveLawIndex(name);
	for (const auto & element : elements_to_add) {
	auto index = this->addElement(element);
	handler.constitutive_law_index(element) = law_id;
	handler.constitutive_law_local_numbering(element) = index;
	}

	this->resizeInternals();

	AKANTU_DEBUG_OUT();
	}


	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	void ConstitutiveLaw::removeElements(const Array<Element> & elements_to_remove) {
	AKANTU_DEBUG_IN();

	auto el_begin = elements_to_remove.begin();
	auto el_end = elements_to_remove.end();

	if (elements_to_remove.empty()) {
	return;
	}

	auto & mesh = handler.getMesh();

	ElementTypeMapArray<UInt> constitutive_law_local_new_numbering(
	"remove constitutive law filter elem", id);

	constitutive_law_local_new_numbering.initialize(
	mesh, _element_filter = &element_filter, _element_kind = _ek_not_defined,
	_with_nb_element = true);

	ElementTypeMapArray<UInt> element_filter_tmp("element_filter_tmp", id);

	element_filter_tmp.initialize(mesh, _element_filter = &element_filter,
	_element_kind = _ek_not_defined);

	ElementTypeMap<UInt> new_ids, element_ids;

	for_each_element(
	mesh,
	[&](auto && el) {
	if (not new_ids(el.type, el.ghost_type)) {
	element_ids(el.type, el.ghost_type) = 0;
	}

	auto & element_id = element_ids(el.type, el.ghost_type);
	auto l_el = Element{el.type, element_id, el.ghost_type};
	if (std::find(el_begin, el_end, el) != el_end) {
	constitutive_law_local_new_numbering(l_el) = UInt(-1);
	return;
	}

	element_filter_tmp(el.type, el.ghost_type).push_back(el.element);
	if (not new_ids(el.type, el.ghost_type)) {
	new_ids(el.type, el.ghost_type) = 0;
	}

	auto & new_id = new_ids(el.type, el.ghost_type);

	constitutive_law_local_new_numbering(l_el) = new_id;
	handler.constitutive_law_local_numbering(el) = new_id;

	++new_id;
	++element_id;
	},
	_element_filter = &element_filter, _element_kind = _ek_not_defined);

	for (auto ghost_type : ghost_types) {
	for (const auto & type : element_filter.elementTypes(
	_ghost_type = ghost_type, _element_kind = _ek_not_defined)) {
	element_filter(type, ghost_type)
	.copy(element_filter_tmp(type, ghost_type));
	}
	}

	for (auto it = internal_vectors_real.begin();
	it != internal_vectors_real.end(); ++it) {
	it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
	}

	for (auto it = internal_vectors_uint.begin();
	it != internal_vectors_uint.end(); ++it) {
	it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
	}

	for (auto it = internal_vectors_bool.begin();
	it != internal_vectors_bool.end(); ++it) {
	it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
	}

	AKANTU_DEBUG_OUT();
	}


	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::onElementsAdded(const Array<Element> & /unused/,
	const NewElementsEvent & /unused/) {
	this->resizeInternals();
	}

	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::onElementsRemoved(
	const Array<Element> & element_list,
	const ElementTypeMapArray<UInt> & new_numbering,
	[[gnu::unused]] const RemovedElementsEvent & event) {

	UInt my_num = handler.getInternalIndexFromID(getID());

	ElementTypeMapArray<UInt> constitutive_law_local_new_numbering(
	"remove constitutive law filter elem", getID());

	auto el_begin = element_list.begin();
	auto el_end = element_list.end();

	for (auto && gt : ghost_types) {
	for (auto && type :
	new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) {

	if (not element_filter.exists(type, gt) \|\|
	element_filter(type, gt).empty()) {
	continue;
	}

	auto & elem_filter = element_filter(type, gt);
	auto & law_indexes = this->handler.constitutive_law_index(type, gt);
	auto & law_loc_num = this->handler.constitutive_law_local_numbering(type, gt);
	auto nb_element = this->handler.getMesh().getNbElement(type, gt);

	// all constitutive laws will resized to the same size...
	law_indexes.resize(nb_element);
	law_loc_num.resize(nb_element);

	if (!constitutive_law_local_new_numbering.exists(type, gt)) {
	constitutive_law_local_new_numbering.alloc(elem_filter.size(), 1, type, gt);
	}

	auto & law_renumbering = constitutive_law_local_new_numbering(type, gt);
	const auto & renumbering = new_numbering(type, gt);
	Array<UInt> elem_filter_tmp;
	UInt ni = 0;
	Element el{type, 0, gt};

	for (UInt i = 0; i < elem_filter.size(); ++i) {
	el.element = elem_filter(i);
	if (std::find(el_begin, el_end, el) == el_end) {
	UInt new_el = renumbering(el.element);
	AKANTU_DEBUG_ASSERT(new_el != UInt(-1),
	"A not removed element as been badly renumbered");
	elem_filter_tmp.push_back(new_el);
	law_renumbering(i) = ni;

	law_indexes(new_el) = my_num;
	law_loc_num(new_el) = ni;
	++ni;
	} else {
	law_renumbering(i) = UInt(-1);
	}
	}

	elem_filter.resize(elem_filter_tmp.size());
	elem_filter.copy(elem_filter_tmp);
	}
	}

	for (auto it = internal_vectors_real.begin();
	it != internal_vectors_real.end(); ++it) {
	it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
	}

	for (auto it = internal_vectors_uint.begin();
	it != internal_vectors_uint.end(); ++it) {
	it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
	}

	for (auto it = internal_vectors_bool.begin();
	it != internal_vectors_bool.end(); ++it) {
	it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
	}
	}



	inline UInt ConstitutiveLaw<ConstitutiveLawsHandler_>::addElement(ElementType type, UInt element,
	GhostType ghost_type) {
	Array<UInt> & el_filter = this->element_filter(type, ghost_type);
	el_filter.push_back(element);
	return el_filter.size() - 1;
	}

	/* -------------------------------------------------------------------------- */
	inline UInt ConstitutiveLaw<ConstitutiveLawsHandler_>::addElement(const Element & element) {
	return this->addElement(element.type, element.element, element.ghost_type);
	}


	/* -------------------------------------------------------------------------- */
	inline const Parameter & ConstitutiveLaw<ConstitutiveLawsHandler_>::getParam(const ID & param) const {
	try {
	return get(param);
	} catch (...) {
	AKANTU_EXCEPTION("No parameter " << param << " in the constitutive law"
	<< getID());
	}
	}

	/* -------------------------------------------------------------------------- */
	template <typename T>
	inline void ConstitutiveLaw<ConstitutiveLawsHandler_>::setParam(const ID & param, T value) {
	try {
	set<T>(param, value);
	} catch (...) {
	AKANTU_EXCEPTION("No parameter " << param << " in the constitutive law "
	<< getID());
	}
	updateInternalParameters();
	}



	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	template <>
	inline void
	ConstitutiveLaw<ConstitutiveLawsHandler_>::registerInternal<Real>(InternalField<Real> & vect) {
	internal_vectors_real[vect.getID()] = &vect;
	}

	template<class ConstitutiveLawsHandler_>
	template <>
	inline void
	ConstitutiveLaw<ConstitutiveLawsHandler_>::registerInternal<UInt>(InternalField<UInt> & vect) {
	internal_vectors_uint[vect.getID()] = &vect;
	}

	template<class ConstitutiveLawsHandler_>
	template <>
	inline void
	ConstitutiveLaw<ConstitutiveLawsHandler_>::registerInternal<bool>(InternalField<bool> & vect) {
	internal_vectors_bool[vect.getID()] = &vect;
	}

	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	template <>
	inline void
	ConstitutiveLaw<ConstitutiveLawsHandler_>::unregisterInternal<Real>(InternalField<Real> & vect) {
	internal_vectors_real.erase(vect.getID());
	}

	template<class ConstitutiveLawsHandler_>
	template <>
	inline void
	ConstitutiveLaw<ConstitutiveLawsHandler_>::unregisterInternal<UInt>(InternalField<UInt> & vect) {
	internal_vectors_uint.erase(vect.getID());
	}

	template<class ConstitutiveLawsHandler_>
	template <>
	inline void
	ConstitutiveLaw<ConstitutiveLawsHandler_>::unregisterInternal<bool>(InternalField<bool> & vect) {
	internal_vectors_bool.erase(vect.getID());
	}



	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	template<typename T>
	const InternalField<T> &
	ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal([[gnu::unused]] const ID & int_id) const {
	AKANTU_TO_IMPLEMENT();
	return NULL;
	}

	/* -------------------------------------------------------------------------- */
	template <typename T>
	InternalField<T> & ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal([[gnu::unused]] const ID & int_id) {
	AKANTU_TO_IMPLEMENT();
	return NULL;
	}

	/* -------------------------------------------------------------------------- */
	template <>
	const InternalField<Real> & ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal(const ID & int_id) const {
	auto it = internal_vectors_real.find(getID() + ":" + int_id);
	if (it == internal_vectors_real.end()) {
	AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
	<< ") does not contain an internal "
	<< int_id << " ("
	<< (getID() + ":" + int_id) << ")");
	}
	return *it->second;
	}

	/* -------------------------------------------------------------------------- */
	template <> InternalField<Real> & ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal(const ID & int_id) {
	auto it = internal_vectors_real.find(getID() + ":" + int_id);
	if (it == internal_vectors_real.end()) {
	AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
	<< ") does not contain an internal "
	<< int_id << " ("
	<< (getID() + ":" + int_id) << ")");
	}
	return *it->second;
	}

	/* -------------------------------------------------------------------------- */
	template <>
	const InternalField<UInt> & ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal(const ID & int_id) const {
	auto it = internal_vectors_uint.find(getID() + ":" + int_id);
	if (it == internal_vectors_uint.end()) {
	AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
	<< ") does not contain an internal "
	<< int_id << " ("
	<< (getID() + ":" + int_id) << ")");
	}
	return *it->second;
	}

	/* -------------------------------------------------------------------------- */
	template <> InternalField<UInt> & ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal(const ID & int_id) {
	auto it = internal_vectors_uint.find(getID() + ":" + int_id);
	if (it == internal_vectors_uint.end()) {
	AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
	<< ") does not contain an internal "
	<< int_id << " ("
	<< (getID() + ":" + int_id) << ")");
	}
	return *it->second;
	}


	/* -------------------------------------------------------------------------- */
	template<class ConstitutiveLawsHandler_>
	inline bool ConstitutiveLaw<ConstitutiveLawsHandler_>::isInternal(const ID & /id/,
	ElementKind /element_kind/) const {
	AKANTU_TO_IMPLEMENT();
	}

	template <class ConstitutiveLawsHandler_>
	template <>
	inline bool
	ConstitutiveLaw<ConstitutiveLawsHandler_>::isInternal<Real>(const ID & id,
	ElementKind element_kind) const {
	auto internal_array = internal_vectors_real.find(this->getID() + ":" + id);

	return not (internal_array == internal_vectors_real.end() \|\|
	internal_array->second->getElementKind() != element_kind);
	}


	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	inline ElementTypeMap<UInt>
	ConstitutiveLaw<ConstitutiveLawsHandler_>::getInternalDataPerElem(const ID & field_id,
	ElementKind element_kind) const {

	if (!this->template isInternal<T>(field_id, element_kind)) {
	AKANTU_EXCEPTION("Cannot find internal field " << id << " in the constitutive law "
	<< this->name);
	}

	const InternalField<T> & internal_field =
	this->template getInternal<T>(field_id);
	const FEEngine & fe_engine = internal_field.getFEEngine();
	UInt nb_data_per_quad = internal_field.getNbComponent();

	ElementTypeMap<UInt> res;
	for (auto ghost_type : ghost_types) {
	for (auto & type : internal_field.elementTypes(ghost_type)) {
	UInt nb_quadrature_points =
	fe_engine.getNbIntegrationPoints(type, ghost_type);
	res(type, ghost_type) = nb_data_per_quad * nb_quadrature_points;
	}
	}

	return res;
	}


	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::flattenInternal(const std::string & field_id,
	ElementTypeMapArray<T> & internal_flat,
	const GhostType ghost_type,
	ElementKind element_kind) const {

	if (!this->template isInternal<T>(field_id, element_kind)) {
	AKANTU_EXCEPTION("Cannot find internal field " << id << " in the constitutive law "
	<< this->name);
	}

	const InternalField<T> & internal_field =
	this->template getInternal<T>(field_id);

	const FEEngine & fe_engine = internal_field.getFEEngine();
	const Mesh & mesh = fe_engine.getMesh();

	for (auto && type : internal_field.filterTypes(ghost_type)) {
	const Array<Real> & src_vect = internal_field(type, ghost_type);
	const Array<UInt> & filter = internal_field.getFilter(type, ghost_type);

	// total number of elements in the corresponding mesh
	UInt nb_element_dst = mesh.getNbElement(type, ghost_type);
	// number of element in the internal field
	UInt nb_element_src = filter.size();
	// number of quadrature points per elem
	UInt nb_quad_per_elem = fe_engine.getNbIntegrationPoints(type);
	// number of data per quadrature point
	UInt nb_data_per_quad = internal_field.getNbComponent();

	if (!internal_flat.exists(type, ghost_type)) {
	internal_flat.alloc(nb_element_dst * nb_quad_per_elem, nb_data_per_quad,
	type, ghost_type);
	}

	if (nb_element_src == 0) {
	continue;
	}

	// number of data per element
	UInt nb_data = nb_quad_per_elem * nb_data_per_quad;

	Array<Real> & dst_vect = internal_flat(type, ghost_type);
	dst_vect.resize(nb_element_dst * nb_quad_per_elem);

	auto it_dst = make_view(dst_vect, nb_data).begin();

	for (auto && data : zip(filter, make_view(src_vect, nb_data))) {
	it_dst[std::get<0>(data)] = std::get<1>(data);
	}
	}
	}

	}

	#endif

constitutive_laws_tmpl.hhNo OneTemporaryActions

File Metadata

constitutive_laws_tmpl.hhView Options

Event Timeline

constitutive_laws_tmpl.hh
No OneTemporary
Actions

constitutive_laws_tmpl.hh
View Options