constitutive_law_tmpl.hh
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, Oct 28, 21:17

constitutive_law_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_law.hh" // NOLINT
	#include "constitutive_laws_handler.hh"
	#include "fe_engine.hh"
	#include "internal_field.hh"
	/* -------------------------------------------------------------------------- */

	#ifndef AKANTU_CONSTITUTIVE_LAW_TMPL_HH
	#define AKANTU_CONSTITUTIVE_LAW_TMPL_HH

	namespace akantu {

	/* -------------------------------------------------------------------------- */
	template <typename T, template <typename Type> class InternalFieldType>
	inline InternalFieldType<T> &
	ConstitutiveLawInternalHandler::registerInternal(const ID & id,
	Int nb_component) {
	return this->registerInternal<T, InternalFieldType>(
	id, nb_component, this->default_fe_engine_id);
	}

	/* -------------------------------------------------------------------------- */
	template <typename T, template <typename Type> class InternalFieldType>
	inline InternalFieldType<T> & ConstitutiveLawInternalHandler::registerInternal(
	const ID & id, Int nb_component, const ID & fe_engine_id) {
	return this->registerInternal<T, InternalFieldType>(
	id, nb_component, fe_engine_id, this->element_filter);
	}
	/* -------------------------------------------------------------------------- */
	template <typename T, template <typename Type> class InternalFieldType>
	inline InternalFieldType<T> & ConstitutiveLawInternalHandler::registerInternal(
	const ID & id, Int nb_component, const ID & fe_engine_id,
	const ElementTypeMapArray<Idx> & element_filter) {
	auto && internal =
	std::shared_ptr<InternalFieldType<T>>(new InternalFieldType<T>(
	id, *this, this->spatial_dimension, fe_engine_id, element_filter));
	internal->initialize(nb_component);
	internal_vectors[internal->getRegisterID()] = internal;
	return *internal;
	}

	/* -------------------------------------------------------------------------- */
	inline void ConstitutiveLawInternalHandler::unregisterInternal(const ID & id) {
	internal_vectors.erase(id);
	}

	/* -------------------------------------------------------------------------- */
	inline void ConstitutiveLawInternalHandler::savePreviousState() {
	for (auto && [_, internal] : internal_vectors) {
	if (internal->hasHistory()) {
	internal->saveCurrentValues();
	}
	}
	}

	/* -------------------------------------------------------------------------- */
	inline void ConstitutiveLawInternalHandler::restorePreviousState() {
	for (auto && [_, internal] : internal_vectors) {
	if (internal->hasHistory()) {
	internal->restorePreviousValues();
	}
	}
	}

	/* -------------------------------------------------------------------------- */
	inline void ConstitutiveLawInternalHandler::resizeInternals() {
	for (auto && [_, internal] : internal_vectors) {
	internal->resize();
	}
	}

	/* -------------------------------------------------------------------------- */
	template <typename T>
	const InternalField<T> &
	ConstitutiveLawInternalHandler::getInternal(const ID & id) const {
	auto it = internal_vectors.find(id);
	if (it != internal_vectors.end() and
	aka::is_of_type<InternalField<T>>(*it->second)) {
	return aka::as_type<InternalField<T>>(*it->second);
	}

	AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
	<< ") does not contain an internal "
	<< id << " (" << (getID() + ":" + id)
	<< ")");
	}

	/* -------------------------------------------------------------------------- */
	template <typename T>
	InternalField<T> & ConstitutiveLawInternalHandler::getInternal(const ID & id) {

	if (auto it = internal_vectors.find(getID() + ":" + id);
	it != internal_vectors.end() and
	aka::is_of_type<InternalField<T>>(*it->second)) {
	return aka::as_type<InternalField<T>>(*it->second);
	}

	AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
	<< ") does not contain an internal "
	<< id << " (" << (getID() + ":" + id)
	<< ")");
	}

	/* -------------------------------------------------------------------------- */
	template <typename T>
	inline bool ConstitutiveLawInternalHandler::isInternal(
	const ID & id, const ElementKind & element_kind) const {
	auto it = internal_vectors.find(id);

	return ((it != internal_vectors.end()) and
	aka::is_of_type<InternalField<T>>(*it->second) and
	aka::as_type<InternalField<T>>(*it->second).getElementKind() ==
	element_kind);
	}

	/* -------------------------------------------------------------------------- */
	/* -------------------------------------------------------------------------- */
	template <typename T>
	const Array<T> &
	ConstitutiveLawInternalHandler::getArray(const ID & vect_id, ElementType type,
	GhostType ghost_type) const {
	if (isInternal<T>(vect_id, Mesh::getKind(type))) {
	auto && internal = this->template getInternal<T>(vect_id);
	if (internal.exists(type, ghost_type)) {
	return internal(type, ghost_type);
	}

	AKANTU_SILENT_EXCEPTION(
	"The internal " << vect_id << " in the constitutive law " << name
	<< " (" << getID() << ") does not contain the type ["
	<< type << ":" << ghost_type << "]");
	}

	AKANTU_SILENT_EXCEPTION("The constitutive law "
	<< name << " (" << getID()
	<< ") does not contain an internal field "
	<< vect_id);
	}

	/* -------------------------------------------------------------------------- */
	template <typename T>
	Array<T> & ConstitutiveLawInternalHandler::getArray(const ID & vect_id,
	ElementType type,
	GhostType ghost_type) {
	if (isInternal<T>(vect_id, Mesh::getKind(type))) {
	auto && internal = this->template getInternal<T>(vect_id);
	if (internal.exists(type, ghost_type)) {
	return internal(type, ghost_type);
	}
	AKANTU_SILENT_EXCEPTION(
	"The internal " << vect_id << " in the constitutive law " << name
	<< " (" << getID() << ") does not contain the type ["
	<< type << ":" << ghost_type << "]");
	}

	AKANTU_SILENT_EXCEPTION("The constitutive law "
	<< name << " (" << getID()
	<< ") does not contain an internal field "
	<< vect_id);
	}

	/* -------------------------------------------------------------------------- */
	inline void ConstitutiveLawInternalHandler::removeIntegrationPoints(
	ElementTypeMapArray<Idx> & new_numbering) {
	for (auto && [_, internal] : internal_vectors) {
	internal->removeIntegrationPoints(new_numbering);
	}
	}
	/* -------------------------------------------------------------------------- */
	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	ConstitutiveLaw<ConstitutiveLawsHandler_>::ConstitutiveLaw(
	ConstitutiveLawsHandler_ & handler, const ID & id, Int spatial_dimension,
	ElementKind element_kind, const ID & fe_engine_id)
	: ConstitutiveLawInternalHandler(id, spatial_dimension, fe_engine_id),
	Parsable(handler.getConstitutiveLawParserType(), id), handler(handler) {

	/// for each connectivity types allocate the element filer array of
	/// the constitutive law
	this->element_filter.initialize(
	handler.getMesh(), _spatial_dimension = handler.getSpatialDimension(),
	_element_kind = element_kind);

	this->initialize();
	}

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

	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::initConstitutiveLaw() {
	this->resizeInternals();
	this->updateInternalParameters();
	is_init = true;
	}

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

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

	this->resizeInternals();

	AKANTU_DEBUG_OUT();
	}

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

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

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

	auto & mesh = handler.getMesh();

	ElementTypeMapArray<Idx> 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<Idx> element_filter_tmp("element_filter_tmp", id);

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

	ElementTypeMap<Idx> 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) = Idx(-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));
	}
	}

	this->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<Idx> & new_numbering,
	const RemovedElementsEvent & /event/) {

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

	ElementTypeMapArray<Idx> 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) or
	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 (not 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<Idx> elem_filter_tmp;
	Int ni = 0;
	Element el{type, 0, gt};

	for (auto && [i, el_id] : enumerate(elem_filter)) {
	el.element = el_id;

	if (std::find(el_begin, el_end, el) == el_end) {
	auto new_el = renumbering(el_id);
	AKANTU_DEBUG_ASSERT(new_el != -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) = -1;
	}
	}

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

	this->removeIntegrationPoints(constitutive_law_local_new_numbering);
	}

	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	template <typename T>
	inline void ConstitutiveLaw<ConstitutiveLawsHandler_>::packInternalFieldHelper(
	const InternalField<T> & data_to_pack, CommunicationBuffer & buffer,
	const Array<Element> & elements) const {
	DataAccessor::packElementalDataHelper<T>(data_to_pack, buffer, elements,
	data_to_pack.getFEEngine());
	}

	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	template <typename T>
	inline void
	ConstitutiveLaw<ConstitutiveLawsHandler_>::unpackInternalFieldHelper(
	InternalField<T> & data_to_unpack, CommunicationBuffer & buffer,
	const Array<Element> & elements) {
	DataAccessor::unpackElementalDataHelper<T>(data_to_unpack, buffer, elements,
	data_to_unpack.getFEEngine());
	}

	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	inline Element ConstitutiveLaw<ConstitutiveLawsHandler_>::convertToLocalElement(
	const Element & global_element) const {
	#ifndef AKANTU_NDEBUG
	auto model_law_index = handler.getConstitutiveLawByElement()(global_element);
	auto 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
	auto local_element = global_element;
	local_element.element =
	handler.getConstitutiveLawLocalNumbering()(global_element);

	return local_element;
	}

	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	inline Element
	ConstitutiveLaw<ConstitutiveLawsHandler_>::convertToGlobalElement(
	const Element & local_element) const {
	auto global_element = local_element;
	global_element.element = this->element_filter(local_element);
	return global_element;
	}

	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	inline Idx
	ConstitutiveLaw<ConstitutiveLawsHandler_>::addElement(const Element & element) {
	auto & el_filter = this->element_filter(element.type, element.ghost_type);
	el_filter.push_back(element.element);
	return el_filter.size() - 1;
	}

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

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

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

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

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

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

	return res;
	}

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

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

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

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

	for (auto && type : internal_field.filterTypes(ghost_type)) {
	const auto & src_vect = internal_field(type, ghost_type);
	const auto & filter = internal_field.getFilter(type, ghost_type);

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

	if (not 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
	auto nb_data = nb_quad_per_elem * nb_data_per_quad;

	auto & 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);
	}
	}
	}

	/* -------------------------------------------------------------------------- */
	template <class ConstitutiveLawsHandler_>
	template <typename T>
	void ConstitutiveLaw<ConstitutiveLawsHandler_>::inflateInternal(
	const std::string & field_id, const ElementTypeMapArray<T> & field,
	GhostType ghost_type, ElementKind element_kind) {
	if (not this->template isInternal<T>(field_id, element_kind)) {
	AKANTU_EXCEPTION("Cannot find internal field " << id << " in material "
	<< this->name);
	}

	auto & internal_field = this->template getInternal<T>(field_id);
	const auto & fe_engine = internal_field.getFEEngine();

	for (auto && type : field.elementTypes(ghost_type)) {
	if (not internal_field.exists(type, ghost_type)) {
	continue;
	}
	const auto & filter = internal_field.getFilter(type, ghost_type);

	const auto & src_array = field(type, ghost_type);
	auto & dest_array = internal_field(type, ghost_type);

	auto nb_quad_per_elem = fe_engine.getNbIntegrationPoints(type);
	auto nb_component = src_array.getNbComponent();

	AKANTU_DEBUG_ASSERT(
	field.size() == fe_engine.getMesh().getNbElement(type, ghost_type) *
	nb_quad_per_elem,
	"The ElementTypeMapArray to inflate is not of the proper size");
	AKANTU_DEBUG_ASSERT(dest_array.getNbComponent() == nb_component,
	"The ElementTypeMapArray has not the proper "
	"number of components");

	auto src =
	make_view(field(type, ghost_type), nb_component, nb_quad_per_elem)
	.begin();
	for (auto && data :
	zip(filter, make_view(dest_array, nb_component, nb_quad_per_elem))) {
	std::get<1>(data) = src[std::get<0>(data)];
	}
	}
	}

	} // namespace akantu

	#endif // AKANTU_CONSTITUTIVE_LAW_TMPL_HH

constitutive_law_tmpl.hhNo OneTemporaryActions

File Metadata

constitutive_law_tmpl.hhView Options

Event Timeline

constitutive_law_tmpl.hh
No OneTemporary
Actions

constitutive_law_tmpl.hh
View Options