cohesive_element_inserter.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, May 29, 07:25

cohesive_element_inserter.cc
View Options

	/**
	* @file cohesive_element_inserter.cc
	*
	* @author Marco Vocialta <marco.vocialta@epfl.ch>
	*
	* @date creation: Wed Dec 04 2013
	* @date last modification: Sun Oct 04 2015
	*
	* @brief Cohesive element inserter functions
	*
	* @section LICENSE
	*
	* Copyright (©) 2014, 2015 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 "cohesive_element_inserter.hh"
	#include "element_group.hh"
	//#include "facet_synchronizer.hh"
	#include "global_ids_updater.hh"
	/* -------------------------------------------------------------------------- */
	#include <algorithm>
	#include <limits>
	/* -------------------------------------------------------------------------- */

	namespace akantu {

	CohesiveElementInserter::CohesiveElementInserter(Mesh & mesh, bool is_extrinsic,
	const ID & id)
	: id(id), mesh(mesh), mesh_facets(mesh.initMeshFacets()),
	insertion_facets("insertion_facets", id),
	insertion_limits(mesh.getSpatialDimension(), 2),
	check_facets("check_facets", id) {
	init(is_extrinsic);
	}

	/* -------------------------------------------------------------------------- */
	CohesiveElementInserter::~CohesiveElementInserter() = default;

	/* -------------------------------------------------------------------------- */
	void CohesiveElementInserter::init(bool is_extrinsic) {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh.getSpatialDimension();

	MeshUtils::resetFacetToDouble(mesh_facets);

	/// initialize facet insertion array
	insertion_facets.initialize(mesh_facets,
	_spatial_dimension = (spatial_dimension - 1),
	_with_nb_element = true);
	// mesh_facets.initElementTypeMapArray(
	// insertion_facets, 1, spatial_dimension - 1, false, _ek_regular, true);

	/// init insertion limits
	for (UInt dim = 0; dim < spatial_dimension; ++dim) {
	insertion_limits(dim, 0) = std::numeric_limits<Real>::max() * (-1.);
	insertion_limits(dim, 1) = std::numeric_limits<Real>::max();
	}

	if (is_extrinsic) {
	check_facets.initialize(mesh_facets,
	_spatial_dimension = spatial_dimension - 1);
	// mesh_facets.initElementTypeMapArray(check_facets, 1, spatial_dimension -
	// 1);
	initFacetsCheck();
	}

	#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
	facet_synchronizer = NULL;
	global_ids_updater = NULL;
	#endif

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void CohesiveElementInserter::initFacetsCheck() {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh.getSpatialDimension();

	for (ghost_type_t::iterator gt = ghost_type_t::begin();
	gt != ghost_type_t::end(); ++gt) {

	GhostType facet_gt = *gt;
	Mesh::type_iterator it =
	mesh_facets.firstType(spatial_dimension - 1, facet_gt);
	Mesh::type_iterator last =
	mesh_facets.lastType(spatial_dimension - 1, facet_gt);

	for (; it != last; ++it) {
	ElementType facet_type = *it;

	Array<bool> & f_check = check_facets(facet_type, facet_gt);

	const Array<std::vector<Element>> & element_to_facet =
	mesh_facets.getElementToSubelement(facet_type, facet_gt);

	UInt nb_facet = element_to_facet.size();
	f_check.resize(nb_facet);

	for (UInt f = 0; f < nb_facet; ++f) {
	if (element_to_facet(f)[1] == ElementNull \|\|
	(element_to_facet(f)[0].ghost_type == _ghost &&
	element_to_facet(f)[1].ghost_type == _ghost)) {
	f_check(f) = false;
	} else
	f_check(f) = true;
	}
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void CohesiveElementInserter::limitCheckFacets() {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh.getSpatialDimension();
	Vector<Real> bary_facet(spatial_dimension);

	for (ghost_type_t::iterator gt = ghost_type_t::begin();
	gt != ghost_type_t::end(); ++gt) {
	GhostType ghost_type = *gt;

	Mesh::type_iterator it =
	mesh_facets.firstType(spatial_dimension - 1, ghost_type);
	Mesh::type_iterator end =
	mesh_facets.lastType(spatial_dimension - 1, ghost_type);
	for (; it != end; ++it) {
	ElementType type = *it;
	Array<bool> & f_check = check_facets(type, ghost_type);
	UInt nb_facet = mesh_facets.getNbElement(type, ghost_type);

	for (UInt f = 0; f < nb_facet; ++f) {
	if (f_check(f)) {

	mesh_facets.getBarycenter(f, type, bary_facet.storage(), ghost_type);

	UInt coord_in_limit = 0;

	while (coord_in_limit < spatial_dimension &&
	bary_facet(coord_in_limit) >
	insertion_limits(coord_in_limit, 0) &&
	bary_facet(coord_in_limit) <
	insertion_limits(coord_in_limit, 1))
	++coord_in_limit;

	if (coord_in_limit != spatial_dimension)
	f_check(f) = false;
	}
	}
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void CohesiveElementInserter::setLimit(SpacialDirection axis, Real first_limit,
	Real second_limit) {

	AKANTU_DEBUG_ASSERT(
	axis < mesh.getSpatialDimension(),
	"You are trying to limit insertion in a direction that doesn't exist");

	insertion_limits(axis, 0) = std::min(first_limit, second_limit);
	insertion_limits(axis, 1) = std::max(first_limit, second_limit);
	}

	/* -------------------------------------------------------------------------- */
	UInt CohesiveElementInserter::insertIntrinsicElements() {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh.getSpatialDimension();

	Vector<Real> bary_facet(spatial_dimension);

	for (auto && ghost_type : ghost_types) {
	for (const auto & type_facet :
	mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) {
	auto & f_insertion = insertion_facets(type_facet, ghost_type);
	auto & element_to_facet =
	mesh_facets.getElementToSubelement(type_facet, ghost_type);

	UInt nb_facet = mesh_facets.getNbElement(type_facet, ghost_type);
	for (UInt f = 0; f < nb_facet; ++f) {

	if (element_to_facet(f)[1] == ElementNull)
	continue;

	mesh_facets.getBarycenter(f, type_facet, bary_facet.storage(),
	ghost_type);

	UInt coord_in_limit = 0;

	while (coord_in_limit < spatial_dimension &&
	bary_facet(coord_in_limit) >
	insertion_limits(coord_in_limit, 0) &&
	bary_facet(coord_in_limit) < insertion_limits(coord_in_limit, 1))
	++coord_in_limit;

	if (coord_in_limit == spatial_dimension)
	f_insertion(f) = true;
	}
	}
	}

	AKANTU_DEBUG_OUT();

	return insertElements();
	}

	/* -------------------------------------------------------------------------- */
	void CohesiveElementInserter::insertIntrinsicElements(
	const std::string & physname, UInt material_index) {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh.getSpatialDimension();
	ElementTypeMapArray<UInt> * phys_data;
	try {
	phys_data = &(mesh_facets.getData<UInt>("physical_names"));
	} catch (...) {
	phys_data = &(mesh_facets.registerData<UInt>("physical_names"));
	phys_data->initialize(mesh_facets,
	_spatial_dimension = spatial_dimension - 1,
	_with_nb_element = true);
	// mesh_facets.initElementTypeMapArray(*phys_data, 1, spatial_dimension - 1,
	// false, _ek_regular, true);
	}
	Vector<Real> bary_facet(spatial_dimension);
	mesh_facets.createElementGroup(physname);

	GhostType ghost_type = _not_ghost;

	Mesh::type_iterator it =
	mesh_facets.firstType(spatial_dimension - 1, ghost_type);
	Mesh::type_iterator end =
	mesh_facets.lastType(spatial_dimension - 1, ghost_type);

	for (; it != end; ++it) {
	const ElementType type_facet = *it;
	Array<bool> & f_insertion = insertion_facets(type_facet, ghost_type);
	Array<std::vector<Element>> & element_to_facet =
	mesh_facets.getElementToSubelement(type_facet, ghost_type);

	UInt nb_facet = mesh_facets.getNbElement(type_facet, ghost_type);
	UInt coord_in_limit = 0;

	ElementGroup & group = mesh.getElementGroup(physname);
	ElementGroup & group_facet = mesh_facets.getElementGroup(physname);

	Vector<Real> bary_physgroup(spatial_dimension);
	Real norm_bary;
	for (ElementGroup::const_element_iterator el_it(
	group.begin(type_facet, ghost_type));
	el_it != group.end(type_facet, ghost_type); ++el_it) {

	UInt e = *el_it;
	mesh.getBarycenter(e, type_facet, bary_physgroup.storage(), ghost_type);
	#ifndef AKANTU_NDEBUG
	bool find_a_partner = false;
	#endif
	norm_bary = bary_physgroup.norm();
	Array<UInt> & material_id = (*phys_data)(type_facet, ghost_type);

	for (UInt f = 0; f < nb_facet; ++f) {

	if (element_to_facet(f)[1] == ElementNull)
	continue;

	mesh_facets.getBarycenter(f, type_facet, bary_facet.storage(),
	ghost_type);

	coord_in_limit = 0;

	while (coord_in_limit < spatial_dimension &&
	(std::abs(bary_facet(coord_in_limit) -
	bary_physgroup(coord_in_limit)) /
	norm_bary <
	Math::getTolerance()))
	++coord_in_limit;

	if (coord_in_limit == spatial_dimension) {
	f_insertion(f) = true;
	#ifndef AKANTU_NDEBUG
	find_a_partner = true;
	#endif
	group_facet.add(type_facet, f, ghost_type, false);
	material_id(f) = material_index;
	break;
	}
	}
	AKANTU_DEBUG_ASSERT(find_a_partner,
	"The element nO "
	<< e << " of physical group " << physname
	<< " did not find its associated facet!"
	<< " Try to decrease math tolerance. "
	<< std::endl);
	}
	}
	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	UInt CohesiveElementInserter::insertElements(bool only_double_facets) {

	CohesiveNewNodesEvent node_event;
	NewElementsEvent element_event;

	Array<UInt> new_pairs(0, 2);

	UInt nb_new_elements = MeshUtils::insertCohesiveElements(
	mesh, mesh_facets, insertion_facets, new_pairs, element_event.getList(),
	only_double_facets);

	UInt nb_new_nodes = new_pairs.size();
	node_event.getList().reserve(nb_new_nodes);
	node_event.getOldNodesList().reserve(nb_new_nodes);
	for (UInt n = 0; n < nb_new_nodes; ++n) {
	node_event.getList().push_back(new_pairs(n, 1));
	node_event.getOldNodesList().push_back(new_pairs(n, 0));
	}

	#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
	if (mesh.isDistributed()) {

	/// update nodes type
	updateNodesType(mesh, node_event);
	updateNodesType(mesh_facets, node_event);

	/// update global ids
	nb_new_nodes = this->updateGlobalIDs(node_event);

	/// compute total number of new elements
	const auto & comm = mesh.getCommunicator();
	comm.allReduce(nb_new_elements, _so_sum);
	}
	#endif

	if (nb_new_nodes > 0)
	mesh.sendEvent(node_event);

	if (nb_new_elements > 0) {
	updateInsertionFacets();
	// mesh.updateTypesOffsets(_not_ghost);
	mesh.sendEvent(element_event);
	MeshUtils::resetFacetToDouble(mesh_facets);
	}

	#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
	if (mesh.isDistributed()) {
	/// update global ids
	this->synchronizeGlobalIDs(node_event);
	}
	#endif

	return nb_new_elements;
	}

	/* -------------------------------------------------------------------------- */
	void CohesiveElementInserter::updateInsertionFacets() {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh.getSpatialDimension();

	for (ghost_type_t::iterator gt = ghost_type_t::begin();
	gt != ghost_type_t::end(); ++gt) {

	GhostType facet_gt = *gt;
	Mesh::type_iterator it =
	mesh_facets.firstType(spatial_dimension - 1, facet_gt);
	Mesh::type_iterator last =
	mesh_facets.lastType(spatial_dimension - 1, facet_gt);

	for (; it != last; ++it) {
	ElementType facet_type = *it;

	Array<bool> & ins_facets = insertion_facets(facet_type, facet_gt);

	// this is the extrinsic case
	if (check_facets.exists(facet_type, facet_gt)) {
	Array<bool> & f_check = check_facets(facet_type, facet_gt);

	UInt nb_facets = f_check.size();

	for (UInt f = 0; f < ins_facets.size(); ++f) {
	if (ins_facets(f)) {
	++nb_facets;
	ins_facets(f) = false;
	f_check(f) = false;
	}
	}

	f_check.resize(nb_facets);
	}
	// and this the intrinsic one
	else {
	ins_facets.resize(mesh_facets.getNbElement(facet_type, facet_gt));
	ins_facets.set(false);
	}
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void CohesiveElementInserter::printself(std::ostream & stream,
	int indent) const {
	std::string space;
	for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
	;

	stream << space << "CohesiveElementInserter [" << std::endl;

	stream << space << " + mesh [" << std::endl;
	mesh.printself(stream, indent + 2);
	stream << space << AKANTU_INDENT << "]" << std::endl;

	stream << space << " + mesh_facets [" << std::endl;
	mesh_facets.printself(stream, indent + 2);
	stream << space << AKANTU_INDENT << "]" << std::endl;

	stream << space << "]" << std::endl;
	}

	} // namespace akantu

cohesive_element_inserter.ccNo OneTemporaryActions

File Metadata

cohesive_element_inserter.ccView Options

Event Timeline

cohesive_element_inserter.cc
No OneTemporary
Actions

cohesive_element_inserter.cc
View Options