mesh_utils.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, Nov 3, 23:06

mesh_utils.cc
View Options

	/**
	* @file mesh_utils.cc
	*
	* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
	* @author Dana Christen <dana.christen@epfl.ch>
	* @author David Simon Kammer <david.kammer@epfl.ch>
	* @author Nicolas Richart <nicolas.richart@epfl.ch>
	* @author Leonardo Snozzi <leonardo.snozzi@epfl.ch>
	* @author Marco Vocialta <marco.vocialta@epfl.ch>
	*
	* @date creation: Fri Aug 20 2010
	* @date last modification: Fri Jan 22 2016
	*
	* @brief All mesh utils necessary for various tasks
	*
	* @section LICENSE
	*
	* Copyright (©) 2010-2012, 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 "mesh_utils.hh"
	#include "element_synchronizer.hh"
	#include "fe_engine.hh"
	#include "mesh_accessor.hh"
	/* -------------------------------------------------------------------------- */
	#include <limits>
	#include <numeric>
	#include <queue>
	#include <set>
	/* -------------------------------------------------------------------------- */

	namespace akantu {

	/* -------------------------------------------------------------------------- */
	void MeshUtils::buildNode2Elements(const Mesh & mesh,
	CSR<Element> & node_to_elem,
	UInt spatial_dimension) {
	AKANTU_DEBUG_IN();
	if (spatial_dimension == _all_dimensions)
	spatial_dimension = mesh.getSpatialDimension();

	/// count number of occurrence of each node
	UInt nb_nodes = mesh.getNbNodes();

	/// array for the node-element list
	node_to_elem.resizeRows(nb_nodes);
	node_to_elem.clearRows();

	AKANTU_DEBUG_ASSERT(
	mesh.firstType(spatial_dimension) != mesh.lastType(spatial_dimension),
	"Some elements must be found in right dimension to compute facets!");

	for (ghost_type_t::iterator gt = ghost_type_t::begin();
	gt != ghost_type_t::end(); ++gt) {
	Mesh::type_iterator first =
	mesh.firstType(spatial_dimension, *gt, _ek_not_defined);
	Mesh::type_iterator last =
	mesh.lastType(spatial_dimension, *gt, _ek_not_defined);

	for (; first != last; ++first) {
	ElementType type = *first;
	UInt nb_element = mesh.getNbElement(type, *gt);
	Array<UInt>::const_iterator<Vector<UInt>> conn_it =
	mesh.getConnectivity(type, *gt).begin(
	Mesh::getNbNodesPerElement(type));

	for (UInt el = 0; el < nb_element; ++el, ++conn_it)
	for (UInt n = 0; n < conn_it->size(); ++n)
	++node_to_elem.rowOffset((*conn_it)(n));
	}
	}

	node_to_elem.countToCSR();
	node_to_elem.resizeCols();

	/// rearrange element to get the node-element list
	Element e;
	node_to_elem.beginInsertions();

	for (ghost_type_t::iterator gt = ghost_type_t::begin();
	gt != ghost_type_t::end(); ++gt) {
	Mesh::type_iterator first =
	mesh.firstType(spatial_dimension, *gt, _ek_not_defined);
	Mesh::type_iterator last =
	mesh.lastType(spatial_dimension, *gt, _ek_not_defined);
	e.ghost_type = *gt;
	for (; first != last; ++first) {
	ElementType type = *first;
	e.type = type;
	e.kind = Mesh::getKind(type);
	UInt nb_element = mesh.getNbElement(type, *gt);
	Array<UInt>::const_iterator<Vector<UInt>> conn_it =
	mesh.getConnectivity(type, *gt).begin(
	Mesh::getNbNodesPerElement(type));

	for (UInt el = 0; el < nb_element; ++el, ++conn_it) {
	e.element = el;
	for (UInt n = 0; n < conn_it->size(); ++n)
	node_to_elem.insertInRow((*conn_it)(n), e);
	}
	}
	}

	node_to_elem.endInsertions();

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	/**
	* This function should disappear in the future (used in mesh partitioning)
	*/
	void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR<UInt> & node_to_elem,
	UInt spatial_dimension) {
	AKANTU_DEBUG_IN();
	if (spatial_dimension == _all_dimensions)
	spatial_dimension = mesh.getSpatialDimension();
	UInt nb_nodes = mesh.getNbNodes();

	const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList();
	Mesh::ConnectivityTypeList::const_iterator it;

	UInt nb_types = type_list.size();
	UInt nb_good_types = 0;

	Vector<UInt> nb_nodes_per_element(nb_types);
	UInt ** conn_val = new UInt *[nb_types];
	Vector<UInt> nb_element(nb_types);

	for (it = type_list.begin(); it != type_list.end(); ++it) {
	ElementType type = *it;
	if (Mesh::getSpatialDimension(type) != spatial_dimension)
	continue;

	nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type);
	conn_val[nb_good_types] = mesh.getConnectivity(type, _not_ghost).storage();
	nb_element[nb_good_types] =
	mesh.getConnectivity(type, _not_ghost).getSize();
	nb_good_types++;
	}

	AKANTU_DEBUG_ASSERT(
	nb_good_types != 0,
	"Some elements must be found in right dimension to compute facets!");

	/// array for the node-element list
	node_to_elem.resizeRows(nb_nodes);
	node_to_elem.clearRows();

	/// count number of occurrence of each node
	for (UInt t = 0; t < nb_good_types; ++t) {
	for (UInt el = 0; el < nb_element[t]; ++el) {
	UInt el_offset = el * nb_nodes_per_element[t];
	for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) {
	++node_to_elem.rowOffset(conn_val[t][el_offset + n]);
	}
	}
	}

	node_to_elem.countToCSR();
	node_to_elem.resizeCols();
	node_to_elem.beginInsertions();

	/// rearrange element to get the node-element list
	for (UInt t = 0, linearized_el = 0; t < nb_good_types; ++t)
	for (UInt el = 0; el < nb_element[t]; ++el, ++linearized_el) {
	UInt el_offset = el * nb_nodes_per_element[t];
	for (UInt n = 0; n < nb_nodes_per_element[t]; ++n)
	node_to_elem.insertInRow(conn_val[t][el_offset + n], linearized_el);
	}

	node_to_elem.endInsertions();
	delete[] conn_val;
	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::buildNode2ElementsElementTypeMap(const Mesh & mesh,
	CSR<UInt> & node_to_elem,
	const ElementType & type,
	const GhostType & ghost_type) {
	AKANTU_DEBUG_IN();
	UInt nb_nodes = mesh.getNbNodes();

	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
	UInt nb_elements = mesh.getConnectivity(type, ghost_type).getSize();

	UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage();

	/// array for the node-element list
	node_to_elem.resizeRows(nb_nodes);
	node_to_elem.clearRows();

	/// count number of occurrence of each node
	for (UInt el = 0; el < nb_elements; ++el) {
	UInt el_offset = el * nb_nodes_per_element;
	for (UInt n = 0; n < nb_nodes_per_element; ++n)
	++node_to_elem.rowOffset(conn_val[el_offset + n]);
	}

	/// convert the occurrence array in a csr one
	node_to_elem.countToCSR();

	node_to_elem.resizeCols();
	node_to_elem.beginInsertions();

	/// save the element index in the node-element list
	for (UInt el = 0; el < nb_elements; ++el) {
	UInt el_offset = el * nb_nodes_per_element;
	for (UInt n = 0; n < nb_nodes_per_element; ++n) {
	node_to_elem.insertInRow(conn_val[el_offset + n], el);
	}
	}

	node_to_elem.endInsertions();

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::buildFacets(Mesh & mesh) {
	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 gt_facet = *gt;
	Mesh::type_iterator it = mesh.firstType(spatial_dimension - 1, gt_facet);
	Mesh::type_iterator end = mesh.lastType(spatial_dimension - 1, gt_facet);
	for (; it != end; ++it) {
	mesh.getConnectivity(it, gt).resize(0);
	// \todo inform the mesh event handler
	}
	}

	buildFacetsDimension(mesh, mesh, true, spatial_dimension);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets,
	UInt to_dimension) {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh.getSpatialDimension();

	buildAllFacets(mesh, mesh_facets, spatial_dimension, to_dimension);

	AKANTU_DEBUG_OUT();
	}
	/* -------------------------------------------------------------------------- */
	void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets,
	UInt from_dimension, UInt to_dimension) {
	AKANTU_DEBUG_IN();

	AKANTU_DEBUG_ASSERT(
	mesh_facets.isMeshFacets(),
	"The mesh_facets should be initialized with initMeshFacets");

	const ElementTypeMapArray<UInt> * prank_to_element = nullptr;

	if (mesh.isDistributed()) {
	prank_to_element = &(mesh.getElementSynchronizer().getPrankToElement());
	}

	/// generate facets
	buildFacetsDimension(mesh, mesh_facets, false, from_dimension,
	prank_to_element);

	/// copy nodes type
	MeshAccessor mesh_accessor_facets(mesh_facets);
	mesh_accessor_facets.getNodesType().copy(mesh.getNodesType());

	/// sort facets and generate sub-facets
	for (UInt i = from_dimension - 1; i > to_dimension; --i) {
	buildFacetsDimension(mesh_facets, mesh_facets, false, i, prank_to_element);
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::buildFacetsDimension(
	const Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension,
	const ElementTypeMapArray<UInt> * prank_to_element) {
	AKANTU_DEBUG_IN();

	// save the current parent of mesh_facets and set it temporarly to mesh since
	// mesh is the one containing the elements for which mesh_facets has the
	// sub-elements
	// example: if the function is called with mesh = mesh_facets
	const Mesh * mesh_facets_parent = nullptr;
	try {
	mesh_facets_parent = &mesh_facets.getMeshParent();
	} catch (...) {
	}

	mesh_facets.defineMeshParent(mesh);
	MeshAccessor mesh_accessor(mesh_facets);

	UInt spatial_dimension = mesh.getSpatialDimension();

	const Array<Real> & mesh_facets_nodes = mesh_facets.getNodes();
	const Array<Real>::const_vector_iterator mesh_facets_nodes_it =
	mesh_facets_nodes.begin(spatial_dimension);

	CSR<Element> node_to_elem;
	buildNode2Elements(mesh, node_to_elem, dimension);

	Array<UInt> counter;
	std::vector<Element> connected_elements;

	// init the SubelementToElement data to improve performance
	for (ghost_type_t::iterator gt = ghost_type_t::begin();
	gt != ghost_type_t::end(); ++gt) {
	GhostType ghost_type = *gt;
	Mesh::type_iterator first = mesh.firstType(dimension, ghost_type);
	Mesh::type_iterator last = mesh.lastType(dimension, ghost_type);

	for (; first != last; ++first) {
	ElementType type = *first;

	mesh_accessor.getSubelementToElement(type, ghost_type);

	Vector<ElementType> facet_types = mesh.getAllFacetTypes(type);

	for (UInt ft = 0; ft < facet_types.size(); ++ft) {
	ElementType facet_type = facet_types(ft);
	mesh_accessor.getElementToSubelement(facet_type, ghost_type);
	mesh_accessor.getConnectivity(facet_type, ghost_type);
	}
	}
	}

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

	current_element.ghost_type = ghost_type;
	Mesh::type_iterator first = mesh.firstType(dimension, ghost_type);
	Mesh::type_iterator last = mesh.lastType(dimension, ghost_type);

	for (; first != last; ++first) {
	ElementType type = *first;
	Vector<ElementType> facet_types = mesh.getAllFacetTypes(type);

	current_element.type = type;

	for (UInt ft = 0; ft < facet_types.size(); ++ft) {
	ElementType facet_type = facet_types(ft);
	UInt nb_element = mesh.getNbElement(type, ghost_type);

	Array<std::vector<Element>> * element_to_subelement =
	&mesh_facets.getElementToSubelement(facet_type, ghost_type);
	Array<UInt> * connectivity_facets =
	&mesh_facets.getConnectivity(facet_type, ghost_type);

	UInt nb_facet_per_element = mesh.getNbFacetsPerElement(type, ft);
	const Array<UInt> & element_connectivity =
	mesh.getConnectivity(type, ghost_type);

	const Matrix<UInt> facet_local_connectivity =
	mesh.getFacetLocalConnectivity(type, ft);
	UInt nb_nodes_per_facet = connectivity_facets->getNbComponent();
	Vector<UInt> facet(nb_nodes_per_facet);

	for (UInt el = 0; el < nb_element; ++el) {
	current_element.element = el;

	for (UInt f = 0; f < nb_facet_per_element; ++f) {
	for (UInt n = 0; n < nb_nodes_per_facet; ++n)
	facet(n) =
	element_connectivity(el, facet_local_connectivity(f, n));

	UInt first_node_nb_elements = node_to_elem.getNbCols(facet(0));
	counter.resize(first_node_nb_elements);
	counter.clear();

	// loop over the other nodes to search intersecting elements,
	// which are the elements that share another node with the
	// starting element after first_node
	CSR<Element>::iterator first_node_elements =
	node_to_elem.begin(facet(0));
	CSR<Element>::iterator first_node_elements_end =
	node_to_elem.end(facet(0));
	UInt local_el = 0;
	for (; first_node_elements != first_node_elements_end;
	++first_node_elements, ++local_el) {
	for (UInt n = 1; n < nb_nodes_per_facet; ++n) {
	CSR<Element>::iterator node_elements_begin =
	node_to_elem.begin(facet(n));
	CSR<Element>::iterator node_elements_end =
	node_to_elem.end(facet(n));
	counter(local_el) +=
	std::count(node_elements_begin, node_elements_end,
	*first_node_elements);
	}
	}

	// counting the number of elements connected to the facets and
	// taking the minimum element number, because the facet should
	// be inserted just once
	UInt nb_element_connected_to_facet = 0;
	Element minimum_el = ElementNull;
	connected_elements.clear();
	for (UInt el_f = 0; el_f < first_node_nb_elements; el_f++) {
	Element real_el = node_to_elem(facet(0), el_f);
	if (counter(el_f) == nb_nodes_per_facet - 1) {
	++nb_element_connected_to_facet;
	minimum_el = std::min(minimum_el, real_el);
	connected_elements.push_back(real_el);
	}
	}

	if (minimum_el == current_element) {
	bool full_ghost_facet = false;

	UInt n = 0;
	while (n < nb_nodes_per_facet && mesh.isPureGhostNode(facet(n)))
	++n;
	if (n == nb_nodes_per_facet)
	full_ghost_facet = true;

	if (!full_ghost_facet) {
	if (!boundary_only \|\|
	(boundary_only && nb_element_connected_to_facet == 1)) {

	std::vector<Element> elements;

	// build elements_on_facets: linearized_el must come first
	// in order to store the facet in the correct direction
	// and avoid to invert the sign in the normal computation
	elements.push_back(current_element);

	/// boundary facet
	if (nb_element_connected_to_facet == 1)
	elements.push_back(ElementNull);
	/// internal facet
	else if (nb_element_connected_to_facet == 2) {
	elements.push_back(connected_elements[1]);

	/// check if facet is in between ghost and normal
	/// elements: if it's the case, the facet is either
	/// ghost or not ghost. The criterion to decide this
	/// is arbitrary. It was chosen to check the processor
	/// id (prank) of the two neighboring elements. If
	/// prank of the ghost element is lower than prank of
	/// the normal one, the facet is not ghost, otherwise
	/// it's ghost
	GhostType gt[2] = {_not_ghost, _not_ghost};
	for (UInt el = 0; el < connected_elements.size(); ++el)
	gt[el] = connected_elements[el].ghost_type;

	if (gt[0] + gt[1] == 1) {
	if (prank_to_element) {
	UInt prank[2];
	for (UInt el = 0; el < 2; ++el) {
	UInt current_el = connected_elements[el].element;
	ElementType current_type =
	connected_elements[el].type;
	GhostType current_gt =
	connected_elements[el].ghost_type;

	const Array<UInt> & prank_to_el =
	(*prank_to_element)(current_type, current_gt);

	prank[el] = prank_to_el(current_el);
	}

	bool ghost_one = (gt[0] != _ghost);

	if (prank[ghost_one] > prank[!ghost_one])
	facet_ghost_type = _not_ghost;
	else
	facet_ghost_type = _ghost;

	connectivity_facets = &mesh_facets.getConnectivity(
	facet_type, facet_ghost_type);
	element_to_subelement =
	&mesh_facets.getElementToSubelement(
	facet_type, facet_ghost_type);
	}
	}
	}
	/// facet of facet
	else {
	for (UInt i = 1; i < nb_element_connected_to_facet; ++i) {
	elements.push_back(connected_elements[i]);
	}
	}

	element_to_subelement->push_back(elements);
	connectivity_facets->push_back(facet);

	/// current facet index
	UInt current_facet = connectivity_facets->getSize() - 1;

	/// loop on every element connected to current facet and
	/// insert current facet in the first free spot of the
	/// subelement_to_element vector
	for (UInt elem = 0; elem < elements.size(); ++elem) {
	Element loc_el = elements[elem];

	if (loc_el.type != _not_defined) {

	Array<Element> & subelement_to_element =
	mesh_facets.getSubelementToElement(loc_el.type,
	loc_el.ghost_type);

	UInt nb_facet_per_loc_element =
	subelement_to_element.getNbComponent();

	for (UInt f_in = 0; f_in < nb_facet_per_loc_element;
	++f_in) {
	if (subelement_to_element(loc_el.element, f_in).type ==
	_not_defined) {
	subelement_to_element(loc_el.element, f_in).type =
	facet_type;
	subelement_to_element(loc_el.element, f_in).element =
	current_facet;
	subelement_to_element(loc_el.element, f_in)
	.ghost_type = facet_ghost_type;
	break;
	}
	}
	}
	}

	/// reset connectivity in case a facet was found in
	/// between ghost and normal elements
	if (facet_ghost_type != ghost_type) {
	facet_ghost_type = ghost_type;
	connectivity_facets = &mesh_accessor.getConnectivity(
	facet_type, facet_ghost_type);
	element_to_subelement =
	&mesh_accessor.getElementToSubelement(facet_type,
	facet_ghost_type);
	}
	}
	}
	}
	}
	}
	}
	}
	}

	// restore the parent of mesh_facet
	if (mesh_facets_parent)
	mesh_facets.defineMeshParent(*mesh_facets_parent);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::renumberMeshNodes(Mesh & mesh,
	Array<UInt> & local_connectivities,
	UInt nb_local_element, UInt nb_ghost_element,
	ElementType type,
	Array<UInt> & old_nodes_numbers) {
	AKANTU_DEBUG_IN();

	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);

	std::map<UInt, UInt> renumbering_map;
	for (UInt i = 0; i < old_nodes_numbers.getSize(); ++i) {
	renumbering_map[old_nodes_numbers(i)] = i;
	}

	/// renumber the nodes
	renumberNodesInConnectivity(local_connectivities,
	(nb_local_element + nb_ghost_element) *
	nb_nodes_per_element,
	renumbering_map);

	std::map<UInt, UInt>::iterator it = renumbering_map.begin();
	std::map<UInt, UInt>::iterator end = renumbering_map.end();
	old_nodes_numbers.resize(renumbering_map.size());
	for (; it != end; ++it) {
	old_nodes_numbers(it->second) = it->first;
	}
	renumbering_map.clear();

	MeshAccessor mesh_accessor(mesh);

	/// copy the renumbered connectivity to the right place
	Array<UInt> & local_conn = mesh_accessor.getConnectivity(type);
	local_conn.resize(nb_local_element);
	memcpy(local_conn.storage(), local_connectivities.storage(),
	nb_local_element * nb_nodes_per_element * sizeof(UInt));

	Array<UInt> & ghost_conn = mesh_accessor.getConnectivity(type, _ghost);
	ghost_conn.resize(nb_ghost_element);
	memcpy(ghost_conn.storage(),
	local_connectivities.storage() +
	nb_local_element * nb_nodes_per_element,
	nb_ghost_element * nb_nodes_per_element * sizeof(UInt));

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::renumberNodesInConnectivity(
	Array<UInt> & list_nodes, UInt nb_nodes,
	std::map<UInt, UInt> & renumbering_map) {
	AKANTU_DEBUG_IN();

	UInt * connectivity = list_nodes.storage();
	UInt new_node_num = renumbering_map.size();
	for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) {
	UInt & node = *connectivity;
	std::map<UInt, UInt>::iterator it = renumbering_map.find(node);
	if (it == renumbering_map.end()) {
	UInt old_node = node;
	renumbering_map[old_node] = new_node_num;
	node = new_node_num;
	++new_node_num;
	} else {
	node = it->second;
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::purifyMesh(Mesh & mesh) {
	AKANTU_DEBUG_IN();

	std::map<UInt, UInt> renumbering_map;

	RemovedNodesEvent remove_nodes(mesh);
	Array<UInt> & nodes_removed = remove_nodes.getList();

	for (UInt gt = _not_ghost; gt <= _ghost; ++gt) {
	GhostType ghost_type = (GhostType)gt;

	Mesh::type_iterator it =
	mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined);
	Mesh::type_iterator end =
	mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined);
	for (; it != end; ++it) {

	ElementType type(*it);
	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);

	Array<UInt> & connectivity = mesh.getConnectivity(type, ghost_type);
	UInt nb_element(connectivity.getSize());

	renumberNodesInConnectivity(
	connectivity, nb_element * nb_nodes_per_element, renumbering_map);
	}
	}

	Array<UInt> & new_numbering = remove_nodes.getNewNumbering();
	std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1));

	std::map<UInt, UInt>::iterator it = renumbering_map.begin();
	std::map<UInt, UInt>::iterator end = renumbering_map.end();
	for (; it != end; ++it) {
	new_numbering(it->first) = it->second;
	}

	for (UInt i = 0; i < new_numbering.getSize(); ++i) {
	if (new_numbering(i) == UInt(-1))
	nodes_removed.push_back(i);
	}

	mesh.sendEvent(remove_nodes);

	AKANTU_DEBUG_OUT();
	}

	#if defined(AKANTU_COHESIVE_ELEMENT)
	/* -------------------------------------------------------------------------- */
	UInt MeshUtils::insertCohesiveElements(
	Mesh & mesh, Mesh & mesh_facets,
	const ElementTypeMapArray<bool> & facet_insertion,
	Array<UInt> & doubled_nodes, Array<Element> & new_elements,
	bool only_double_facets) {
	UInt spatial_dimension = mesh.getSpatialDimension();
	UInt elements_to_insert = updateFacetToDouble(mesh_facets, facet_insertion);

	if (elements_to_insert > 0) {

	if (spatial_dimension == 1) {
	doublePointFacet(mesh, mesh_facets, doubled_nodes);
	} else {
	doubleFacet(mesh, mesh_facets, spatial_dimension - 1, doubled_nodes,
	true);
	findSubfacetToDouble<false>(mesh, mesh_facets);

	if (spatial_dimension == 2) {
	doubleSubfacet<2>(mesh, mesh_facets, doubled_nodes);
	} else if (spatial_dimension == 3) {
	doubleFacet(mesh, mesh_facets, 1, doubled_nodes, false);
	findSubfacetToDouble<true>(mesh, mesh_facets);
	doubleSubfacet<3>(mesh, mesh_facets, doubled_nodes);
	}
	}

	if (!only_double_facets)
	updateCohesiveData(mesh, mesh_facets, new_elements);
	}

	return elements_to_insert;
	}
	#endif

	/* -------------------------------------------------------------------------- */
	void MeshUtils::doubleNodes(Mesh & mesh, const std::vector<UInt> & old_nodes,
	Array<UInt> & doubled_nodes) {
	AKANTU_DEBUG_IN();

	Array<Real> & position = mesh.getNodes();
	UInt spatial_dimension = mesh.getSpatialDimension();

	UInt old_nb_nodes = position.getSize();
	UInt new_nb_nodes = old_nb_nodes + old_nodes.size();

	UInt old_nb_doubled_nodes = doubled_nodes.getSize();
	UInt new_nb_doubled_nodes = old_nb_doubled_nodes + old_nodes.size();

	position.resize(new_nb_nodes);
	doubled_nodes.resize(new_nb_doubled_nodes);

	Array<Real>::iterator<Vector<Real>> position_begin =
	position.begin(spatial_dimension);

	for (UInt n = 0; n < old_nodes.size(); ++n) {
	UInt new_node = old_nb_nodes + n;

	/// store doubled nodes
	doubled_nodes(old_nb_doubled_nodes + n, 0) = old_nodes[n];
	doubled_nodes(old_nb_doubled_nodes + n, 1) = new_node;

	/// update position
	std::copy(position_begin + old_nodes[n], position_begin + old_nodes[n] + 1,
	position_begin + new_node);
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::doubleFacet(Mesh & mesh, Mesh & mesh_facets,
	UInt facet_dimension, Array<UInt> & doubled_nodes,
	bool facet_mode) {
	AKANTU_DEBUG_IN();

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

	GhostType gt_facet = *gt;

	Mesh::type_iterator it = mesh_facets.firstType(facet_dimension, gt_facet);
	Mesh::type_iterator end = mesh_facets.lastType(facet_dimension, gt_facet);

	for (; it != end; ++it) {

	ElementType type_facet = *it;

	Array<UInt> & f_to_double =
	mesh_facets.getData<UInt>("facet_to_double", type_facet, gt_facet);
	UInt nb_facet_to_double = f_to_double.getSize();

	if (nb_facet_to_double == 0)
	continue;

	ElementType type_subfacet = Mesh::getFacetType(type_facet);
	const UInt nb_subfacet_per_facet =
	Mesh::getNbFacetsPerElement(type_facet);
	GhostType gt_subfacet = _casper;
	Array<std::vector<Element>> * f_to_subfacet = NULL;

	Array<Element> & subfacet_to_facet =
	mesh_facets.getSubelementToElement(type_facet, gt_facet);

	Array<UInt> & conn_facet =
	mesh_facets.getConnectivity(type_facet, gt_facet);
	UInt nb_nodes_per_facet = conn_facet.getNbComponent();

	UInt old_nb_facet = conn_facet.getSize();
	UInt new_nb_facet = old_nb_facet + nb_facet_to_double;

	conn_facet.resize(new_nb_facet);
	subfacet_to_facet.resize(new_nb_facet);

	UInt new_facet = old_nb_facet - 1;
	Element new_facet_el(type_facet, 0, gt_facet);

	Array<Element>::iterator<Vector<Element>> subfacet_to_facet_begin =
	subfacet_to_facet.begin(nb_subfacet_per_facet);
	Array<UInt>::iterator<Vector<UInt>> conn_facet_begin =
	conn_facet.begin(nb_nodes_per_facet);

	for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {
	UInt old_facet = f_to_double(facet);
	++new_facet;

	/// adding a new facet by copying original one

	/// copy connectivity in new facet
	std::copy(conn_facet_begin + old_facet,
	conn_facet_begin + old_facet + 1,
	conn_facet_begin + new_facet);

	/// update subfacet_to_facet
	std::copy(subfacet_to_facet_begin + old_facet,
	subfacet_to_facet_begin + old_facet + 1,
	subfacet_to_facet_begin + new_facet);

	new_facet_el.element = new_facet;

	/// loop on every subfacet
	for (UInt sf = 0; sf < nb_subfacet_per_facet; ++sf) {
	Element & subfacet = subfacet_to_facet(old_facet, sf);
	if (subfacet == ElementNull)
	continue;

	if (gt_subfacet != subfacet.ghost_type) {
	gt_subfacet = subfacet.ghost_type;
	f_to_subfacet = &mesh_facets.getElementToSubelement(
	type_subfacet, subfacet.ghost_type);
	}

	/// update facet_to_subfacet array
	(*f_to_subfacet)(subfacet.element).push_back(new_facet_el);
	}
	}

	/// update facet_to_subfacet and _segment_3 facets if any
	if (!facet_mode) {
	updateSubfacetToFacet(mesh_facets, type_facet, gt_facet, true);
	updateFacetToSubfacet(mesh_facets, type_facet, gt_facet, true);
	updateQuadraticSegments<true>(mesh, mesh_facets, type_facet, gt_facet,
	doubled_nodes);
	} else
	updateQuadraticSegments<false>(mesh, mesh_facets, type_facet, gt_facet,
	doubled_nodes);
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	UInt MeshUtils::updateFacetToDouble(
	Mesh & mesh_facets, const ElementTypeMapArray<bool> & facet_insertion) {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh_facets.getSpatialDimension();
	UInt nb_facets_to_double = 0.;

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

	GhostType gt_facet = *gt;

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

	for (; it != end; ++it) {

	ElementType type_facet = *it;

	const Array<bool> & f_insertion = facet_insertion(type_facet, gt_facet);
	Array<UInt> & f_to_double =
	mesh_facets.getData<UInt>("facet_to_double", type_facet, gt_facet);

	Array<std::vector<Element>> & element_to_facet =
	mesh_facets.getElementToSubelement(type_facet, gt_facet);

	ElementType el_type = _not_defined;
	GhostType el_gt = _casper;
	UInt nb_facet_per_element = 0;
	Element old_facet_el(type_facet, 0, gt_facet);

	Array<Element> * facet_to_element = NULL;

	for (UInt f = 0; f < f_insertion.getSize(); ++f) {

	if (f_insertion(f) == false)
	continue;

	++nb_facets_to_double;

	if (element_to_facet(f)[1].type == _not_defined
	#if defined(AKANTU_COHESIVE_ELEMENT)
	\|\| element_to_facet(f)[1].kind == _ek_cohesive
	#endif
	) {
	AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary");
	continue;
	}

	f_to_double.push_back(f);

	UInt new_facet = mesh_facets.getNbElement(type_facet, gt_facet) +
	f_to_double.getSize() - 1;
	old_facet_el.element = f;

	/// update facet_to_element vector
	Element & elem_to_update = element_to_facet(f)[1];
	UInt el = elem_to_update.element;

	if (elem_to_update.ghost_type != el_gt \|\|
	elem_to_update.type != el_type) {
	el_type = elem_to_update.type;
	el_gt = elem_to_update.ghost_type;
	facet_to_element =
	&mesh_facets.getSubelementToElement(el_type, el_gt);
	nb_facet_per_element = facet_to_element->getNbComponent();
	}

	Element * f_update =
	std::find(facet_to_element->storage() + el * nb_facet_per_element,
	facet_to_element->storage() + el * nb_facet_per_element +
	nb_facet_per_element,
	old_facet_el);

	AKANTU_DEBUG_ASSERT(
	facet_to_element->storage() + el * nb_facet_per_element !=
	facet_to_element->storage() + el * nb_facet_per_element +
	nb_facet_per_element,
	"Facet not found");

	f_update->element = new_facet;

	/// update elements connected to facet
	std::vector<Element> first_facet_list = element_to_facet(f);
	element_to_facet.push_back(first_facet_list);

	/// set new and original facets as boundary facets
	element_to_facet(new_facet)[0] = element_to_facet(new_facet)[1];

	element_to_facet(f)[1] = ElementNull;
	element_to_facet(new_facet)[1] = ElementNull;
	}
	}
	}

	AKANTU_DEBUG_OUT();
	return nb_facets_to_double;
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::resetFacetToDouble(Mesh & mesh_facets) {
	AKANTU_DEBUG_IN();

	for (UInt g = _not_ghost; g <= _ghost; ++g) {
	GhostType gt = (GhostType)g;

	Mesh::type_iterator it = mesh_facets.firstType(_all_dimensions, gt);
	Mesh::type_iterator end = mesh_facets.lastType(_all_dimensions, gt);
	for (; it != end; ++it) {
	ElementType type = *it;
	mesh_facets.getDataPointer<UInt>("facet_to_double", type, gt, 1, false);

	mesh_facets.getDataPointer<std::vector<Element>>(
	"facets_to_subfacet_double", type, gt, 1, false);

	mesh_facets.getDataPointer<std::vector<Element>>(
	"elements_to_subfacet_double", type, gt, 1, false);

	mesh_facets.getDataPointer<std::vector<Element>>(
	"subfacets_to_subsubfacet_double", type, gt, 1, false);
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	template <bool subsubfacet_mode>
	void MeshUtils::findSubfacetToDouble(Mesh & mesh, Mesh & mesh_facets) {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh_facets.getSpatialDimension();
	if (spatial_dimension == 1)
	return;

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

	GhostType gt_facet = *gt;

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

	for (; it != end; ++it) {

	ElementType type_facet = *it;

	Array<UInt> & f_to_double =
	mesh_facets.getData<UInt>("facet_to_double", type_facet, gt_facet);
	UInt nb_facet_to_double = f_to_double.getSize();
	if (nb_facet_to_double == 0)
	continue;

	ElementType type_subfacet = Mesh::getFacetType(type_facet);
	GhostType gt_subfacet = _casper;

	ElementType type_subsubfacet = Mesh::getFacetType(type_subfacet);
	GhostType gt_subsubfacet = _casper;

	Array<UInt> * conn_subfacet = NULL;
	Array<UInt> * sf_to_double = NULL;
	Array<std::vector<Element>> * sf_to_subfacet_double = NULL;
	Array<std::vector<Element>> * f_to_subfacet_double = NULL;
	Array<std::vector<Element>> * el_to_subfacet_double = NULL;

	UInt nb_subfacet = Mesh::getNbFacetsPerElement(type_facet);

	UInt nb_subsubfacet;
	UInt nb_nodes_per_sf_el;

	if (subsubfacet_mode) {
	nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subsubfacet);
	nb_subsubfacet = Mesh::getNbFacetsPerElement(type_subfacet);
	} else
	nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subfacet);

	Array<Element> & subfacet_to_facet =
	mesh_facets.getSubelementToElement(type_facet, gt_facet);

	Array<std::vector<Element>> & element_to_facet =
	mesh_facets.getElementToSubelement(type_facet, gt_facet);

	Array<Element> * subsubfacet_to_subfacet = NULL;

	UInt old_nb_facet = subfacet_to_facet.getSize() - nb_facet_to_double;

	Element current_facet(type_facet, 0, gt_facet);
	std::vector<Element> element_list;
	std::vector<Element> facet_list;
	std::vector<Element> * subfacet_list;
	if (subsubfacet_mode)
	subfacet_list = new std::vector<Element>;

	/// map to filter subfacets
	Array<std::vector<Element>> * facet_to_subfacet = NULL;

	/// this is used to make sure that both new and old facets are
	/// checked
	UInt facets[2];

	/// loop on every facet
	for (UInt f_index = 0; f_index < 2; ++f_index) {
	for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {
	facets[bool(f_index)] = f_to_double(facet);
	facets[!bool(f_index)] = old_nb_facet + facet;

	UInt old_facet = facets[0];
	UInt new_facet = facets[1];

	Element & starting_element = element_to_facet(new_facet)[0];
	current_facet.element = old_facet;

	/// loop on every subfacet
	for (UInt sf = 0; sf < nb_subfacet; ++sf) {

	Element & subfacet = subfacet_to_facet(old_facet, sf);
	if (subfacet == ElementNull)
	continue;

	if (gt_subfacet != subfacet.ghost_type) {
	gt_subfacet = subfacet.ghost_type;

	if (subsubfacet_mode) {
	subsubfacet_to_subfacet = &mesh_facets.getSubelementToElement(
	type_subfacet, gt_subfacet);
	} else {
	conn_subfacet =
	&mesh_facets.getConnectivity(type_subfacet, gt_subfacet);
	sf_to_double = &mesh_facets.getData<UInt>(
	"facet_to_double", type_subfacet, gt_subfacet);

	f_to_subfacet_double =
	&mesh_facets.getData<std::vector<Element>>(
	"facets_to_subfacet_double", type_subfacet,
	gt_subfacet);

	el_to_subfacet_double =
	&mesh_facets.getData<std::vector<Element>>(
	"elements_to_subfacet_double", type_subfacet,
	gt_subfacet);

	facet_to_subfacet = &mesh_facets.getElementToSubelement(
	type_subfacet, gt_subfacet);
	}
	}

	if (subsubfacet_mode) {
	/// loop on every subsubfacet
	for (UInt ssf = 0; ssf < nb_subsubfacet; ++ssf) {
	Element & subsubfacet =
	(*subsubfacet_to_subfacet)(subfacet.element, ssf);

	if (subsubfacet == ElementNull)
	continue;

	if (gt_subsubfacet != subsubfacet.ghost_type) {
	gt_subsubfacet = subsubfacet.ghost_type;
	conn_subfacet = &mesh_facets.getConnectivity(type_subsubfacet,
	gt_subsubfacet);
	sf_to_double = &mesh_facets.getData<UInt>(
	"facet_to_double", type_subsubfacet, gt_subsubfacet);

	sf_to_subfacet_double =
	&mesh_facets.getData<std::vector<Element>>(
	"subfacets_to_subsubfacet_double", type_subsubfacet,
	gt_subsubfacet);

	f_to_subfacet_double =
	&mesh_facets.getData<std::vector<Element>>(
	"facets_to_subfacet_double", type_subsubfacet,
	gt_subsubfacet);

	el_to_subfacet_double =
	&mesh_facets.getData<std::vector<Element>>(
	"elements_to_subfacet_double", type_subsubfacet,
	gt_subsubfacet);

	facet_to_subfacet = &mesh_facets.getElementToSubelement(
	type_subsubfacet, gt_subsubfacet);
	}

	UInt global_ssf = subsubfacet.element;

	Vector<UInt> subsubfacet_connectivity(
	conn_subfacet->storage() + global_ssf * nb_nodes_per_sf_el,
	nb_nodes_per_sf_el);

	/// check if subsubfacet is to be doubled
	if (findElementsAroundSubfacet<true>(
	mesh, mesh_facets, starting_element, current_facet,
	subsubfacet_connectivity, element_list, facet_list,
	subfacet_list) == false &&
	removeElementsInVector(*subfacet_list,
	(*facet_to_subfacet)(global_ssf)) ==
	false) {

	sf_to_double->push_back(global_ssf);
	sf_to_subfacet_double->push_back(*subfacet_list);
	f_to_subfacet_double->push_back(facet_list);
	el_to_subfacet_double->push_back(element_list);
	}
	}
	} else {
	const UInt global_sf = subfacet.element;

	Vector<UInt> subfacet_connectivity(
	conn_subfacet->storage() + global_sf * nb_nodes_per_sf_el,
	nb_nodes_per_sf_el);

	/// check if subfacet is to be doubled
	if (findElementsAroundSubfacet<false>(
	mesh, mesh_facets, starting_element, current_facet,
	subfacet_connectivity, element_list,
	facet_list) == false &&
	removeElementsInVector(
	facet_list, (*facet_to_subfacet)(global_sf)) == false) {

	sf_to_double->push_back(global_sf);
	f_to_subfacet_double->push_back(facet_list);
	el_to_subfacet_double->push_back(element_list);
	}
	}
	}
	}
	}

	if (subsubfacet_mode)
	delete subfacet_list;
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	#if defined(AKANTU_COHESIVE_ELEMENT)
	void MeshUtils::updateCohesiveData(Mesh & mesh, Mesh & mesh_facets,
	Array<Element> & new_elements) {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh.getSpatialDimension();
	bool third_dimension = spatial_dimension == 3;

	MeshAccessor mesh_facets_accessor(mesh_facets);

	for (auto gt_facet : ghost_types) {
	for (auto type_facet :
	mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {

	Array<UInt> & f_to_double =
	mesh_facets.getData<UInt>("facet_to_double", type_facet, gt_facet);

	UInt nb_facet_to_double = f_to_double.getSize();
	if (nb_facet_to_double == 0)
	continue;

	ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet);

	auto & facet_to_coh_element =
	mesh_facets_accessor.getSubelementToElement(type_cohesive, gt_facet);

	auto & conn_facet =
	mesh_facets.getConnectivity(type_facet, gt_facet);
	auto & conn_cohesive =
	mesh.getConnectivity(type_cohesive, gt_facet);
	UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(type_facet);

	Array<std::vector<Element>> & element_to_facet =
	mesh_facets.getElementToSubelement(type_facet, gt_facet);

	UInt old_nb_cohesive_elements = conn_cohesive.getSize();
	UInt new_nb_cohesive_elements =
	conn_cohesive.getSize() + nb_facet_to_double;

	UInt old_nb_facet = element_to_facet.getSize() - nb_facet_to_double;
	facet_to_coh_element.resize(new_nb_cohesive_elements);
	conn_cohesive.resize(new_nb_cohesive_elements);

	UInt new_elements_old_size = new_elements.getSize();
	new_elements.resize(new_elements_old_size + nb_facet_to_double);

	Element c_element(type_cohesive, 0, gt_facet, _ek_cohesive);
	Element f_element(type_facet, 0, gt_facet);

	UInt facets[2];

	for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {

	/// (in 3D cohesive elements connectivity is inverted)
	facets[third_dimension] = f_to_double(facet);
	facets[!third_dimension] = old_nb_facet + facet;

	UInt cohesive_element = old_nb_cohesive_elements + facet;

	/// store doubled facets
	f_element.element = facets[0];
	facet_to_coh_element(cohesive_element, 0) = f_element;
	f_element.element = facets[1];
	facet_to_coh_element(cohesive_element, 1) = f_element;

	/// modify cohesive elements' connectivity
	for (UInt n = 0; n < nb_nodes_per_facet; ++n) {
	conn_cohesive(cohesive_element, n) = conn_facet(facets[0], n);
	conn_cohesive(cohesive_element, n + nb_nodes_per_facet) =
	conn_facet(facets[1], n);
	}

	/// update element_to_facet vectors
	c_element.element = cohesive_element;
	element_to_facet(facets[0])[1] = c_element;
	element_to_facet(facets[1])[1] = c_element;

	/// add cohesive element to the element event list
	new_elements(new_elements_old_size + facet) = c_element;
	}
	}
	}

	AKANTU_DEBUG_OUT();
	}
	#endif

	/* -------------------------------------------------------------------------- */
	void MeshUtils::doublePointFacet(Mesh & mesh, Mesh & mesh_facets,
	Array<UInt> & doubled_nodes) {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh.getSpatialDimension();
	if (spatial_dimension != 1)
	return;

	Array<Real> & position = mesh.getNodes();

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

	GhostType gt_facet = *gt;

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

	for (; it != end; ++it) {

	ElementType type_facet = *it;

	Array<UInt> & conn_facet =
	mesh_facets.getConnectivity(type_facet, gt_facet);
	Array<std::vector<Element>> & element_to_facet =
	mesh_facets.getElementToSubelement(type_facet, gt_facet);

	const Array<UInt> & f_to_double =
	mesh_facets.getData<UInt>("facet_to_double", type_facet, gt_facet);

	UInt nb_facet_to_double = f_to_double.getSize();

	UInt old_nb_facet = element_to_facet.getSize() - nb_facet_to_double;
	UInt new_nb_facet = element_to_facet.getSize();

	UInt old_nb_nodes = position.getSize();
	UInt new_nb_nodes = old_nb_nodes + nb_facet_to_double;

	position.resize(new_nb_nodes);
	conn_facet.resize(new_nb_facet);

	UInt old_nb_doubled_nodes = doubled_nodes.getSize();
	doubled_nodes.resize(old_nb_doubled_nodes + nb_facet_to_double);

	for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {
	UInt old_facet = f_to_double(facet);
	UInt new_facet = old_nb_facet + facet;

	ElementType type = element_to_facet(new_facet)[0].type;
	UInt el = element_to_facet(new_facet)[0].element;
	GhostType gt = element_to_facet(new_facet)[0].ghost_type;

	UInt old_node = conn_facet(old_facet);
	UInt new_node = old_nb_nodes + facet;

	/// update position
	position(new_node) = position(old_node);

	conn_facet(new_facet) = new_node;
	Array<UInt> & conn_segment = mesh.getConnectivity(type, gt);
	UInt nb_nodes_per_segment = conn_segment.getNbComponent();

	/// update facet connectivity
	UInt i;
	for (i = 0;
	conn_segment(el, i) != old_node && i <= nb_nodes_per_segment; ++i)
	;
	conn_segment(el, i) = new_node;

	doubled_nodes(old_nb_doubled_nodes + facet, 0) = old_node;
	doubled_nodes(old_nb_doubled_nodes + facet, 1) = new_node;
	}
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	template <bool third_dim_segments>
	void MeshUtils::updateQuadraticSegments(Mesh & mesh, Mesh & mesh_facets,
	ElementType type_facet,
	GhostType gt_facet,
	Array<UInt> & doubled_nodes) {
	AKANTU_DEBUG_IN();

	if (type_facet != _segment_3)
	return;

	Array<UInt> & f_to_double =
	mesh_facets.getData<UInt>("facet_to_double", type_facet, gt_facet);
	UInt nb_facet_to_double = f_to_double.getSize();

	UInt old_nb_facet =
	mesh_facets.getNbElement(type_facet, gt_facet) - nb_facet_to_double;

	Array<UInt> & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet);

	Array<std::vector<Element>> & element_to_facet =
	mesh_facets.getElementToSubelement(type_facet, gt_facet);

	/// this ones matter only for segments in 3D
	Array<std::vector<Element>> * el_to_subfacet_double = NULL;
	Array<std::vector<Element>> * f_to_subfacet_double = NULL;

	if (third_dim_segments) {
	el_to_subfacet_double = &mesh_facets.getData<std::vector<Element>>(
	"elements_to_subfacet_double", type_facet, gt_facet);

	f_to_subfacet_double = &mesh_facets.getData<std::vector<Element>>(
	"facets_to_subfacet_double", type_facet, gt_facet);
	}

	std::vector<UInt> middle_nodes;

	for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {
	UInt old_facet = f_to_double(facet);
	UInt node = conn_facet(old_facet, 2);
	if (!mesh.isPureGhostNode(node))
	middle_nodes.push_back(node);
	}

	UInt n = doubled_nodes.getSize();

	doubleNodes(mesh, middle_nodes, doubled_nodes);

	for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {
	UInt old_facet = f_to_double(facet);
	UInt old_node = conn_facet(old_facet, 2);
	if (mesh.isPureGhostNode(old_node))
	continue;

	UInt new_node = doubled_nodes(n, 1);
	UInt new_facet = old_nb_facet + facet;

	conn_facet(new_facet, 2) = new_node;

	if (third_dim_segments) {
	updateElementalConnectivity(mesh_facets, old_node, new_node,
	element_to_facet(new_facet));

	updateElementalConnectivity(mesh, old_node, new_node,
	(*el_to_subfacet_double)(facet),
	&(*f_to_subfacet_double)(facet));
	} else {
	updateElementalConnectivity(mesh, old_node, new_node,
	element_to_facet(new_facet));
	}
	++n;
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::updateSubfacetToFacet(Mesh & mesh_facets,
	ElementType type_subfacet,
	GhostType gt_subfacet, bool facet_mode) {
	AKANTU_DEBUG_IN();

	Array<UInt> & sf_to_double =
	mesh_facets.getData<UInt>("facet_to_double", type_subfacet, gt_subfacet);
	UInt nb_subfacet_to_double = sf_to_double.getSize();

	/// update subfacet_to_facet vector
	ElementType type_facet = _not_defined;
	GhostType gt_facet = _casper;
	Array<Element> * subfacet_to_facet = NULL;
	UInt nb_subfacet_per_facet = 0;
	UInt old_nb_subfacet = mesh_facets.getNbElement(type_subfacet, gt_subfacet) -
	nb_subfacet_to_double;

	Array<std::vector<Element>> * facet_list = NULL;
	if (facet_mode)
	facet_list = &mesh_facets.getData<std::vector<Element>>(
	"facets_to_subfacet_double", type_subfacet, gt_subfacet);
	else
	facet_list = &mesh_facets.getData<std::vector<Element>>(
	"subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet);

	Element old_subfacet_el(type_subfacet, 0, gt_subfacet);
	Element new_subfacet_el(type_subfacet, 0, gt_subfacet);

	for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) {
	old_subfacet_el.element = sf_to_double(sf);
	new_subfacet_el.element = old_nb_subfacet + sf;

	for (UInt f = 0; f < (*facet_list)(sf).size(); ++f) {
	Element & facet = (*facet_list)(sf)[f];

	if (facet.type != type_facet \|\| facet.ghost_type != gt_facet) {
	type_facet = facet.type;
	gt_facet = facet.ghost_type;

	subfacet_to_facet =
	&mesh_facets.getSubelementToElement(type_facet, gt_facet);
	nb_subfacet_per_facet = subfacet_to_facet->getNbComponent();
	}

	Element * sf_update = std::find(
	subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet,
	subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet +
	nb_subfacet_per_facet,
	old_subfacet_el);

	AKANTU_DEBUG_ASSERT(subfacet_to_facet->storage() +
	facet.element * nb_subfacet_per_facet !=
	subfacet_to_facet->storage() +
	facet.element * nb_subfacet_per_facet +
	nb_subfacet_per_facet,
	"Subfacet not found");

	*sf_update = new_subfacet_el;
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::updateFacetToSubfacet(Mesh & mesh_facets,
	ElementType type_subfacet,
	GhostType gt_subfacet, bool facet_mode) {
	AKANTU_DEBUG_IN();

	Array<UInt> & sf_to_double =
	mesh_facets.getData<UInt>("facet_to_double", type_subfacet, gt_subfacet);
	UInt nb_subfacet_to_double = sf_to_double.getSize();

	Array<std::vector<Element>> & facet_to_subfacet =
	mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet);

	Array<std::vector<Element>> * facet_to_subfacet_double = NULL;

	if (facet_mode) {
	facet_to_subfacet_double = &mesh_facets.getData<std::vector<Element>>(
	"facets_to_subfacet_double", type_subfacet, gt_subfacet);
	} else {
	facet_to_subfacet_double = &mesh_facets.getData<std::vector<Element>>(
	"subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet);
	}

	UInt old_nb_subfacet = facet_to_subfacet.getSize();
	facet_to_subfacet.resize(old_nb_subfacet + nb_subfacet_to_double);

	for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf)
	facet_to_subfacet(old_nb_subfacet + sf) = (*facet_to_subfacet_double)(sf);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	template <UInt spatial_dimension>
	void MeshUtils::doubleSubfacet(Mesh & mesh, Mesh & mesh_facets,
	Array<UInt> & doubled_nodes) {
	AKANTU_DEBUG_IN();

	if (spatial_dimension == 1)
	return;

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

	GhostType gt_subfacet = *gt;

	Mesh::type_iterator it = mesh_facets.firstType(0, gt_subfacet);
	Mesh::type_iterator end = mesh_facets.lastType(0, gt_subfacet);

	for (; it != end; ++it) {

	ElementType type_subfacet = *it;

	Array<UInt> & sf_to_double = mesh_facets.getData<UInt>(
	"facet_to_double", type_subfacet, gt_subfacet);
	UInt nb_subfacet_to_double = sf_to_double.getSize();

	if (nb_subfacet_to_double == 0)
	continue;

	AKANTU_DEBUG_ASSERT(
	type_subfacet == _point_1,
	"Only _point_1 subfacet doubling is supported at the moment");

	Array<UInt> & conn_subfacet =
	mesh_facets.getConnectivity(type_subfacet, gt_subfacet);

	UInt old_nb_subfacet = conn_subfacet.getSize();
	UInt new_nb_subfacet = old_nb_subfacet + nb_subfacet_to_double;

	conn_subfacet.resize(new_nb_subfacet);

	std::vector<UInt> nodes_to_double;
	UInt old_nb_doubled_nodes = doubled_nodes.getSize();

	/// double nodes
	for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) {
	UInt old_subfacet = sf_to_double(sf);
	nodes_to_double.push_back(conn_subfacet(old_subfacet));
	}

	doubleNodes(mesh, nodes_to_double, doubled_nodes);

	/// add new nodes in connectivity
	for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) {
	UInt new_subfacet = old_nb_subfacet + sf;
	UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1);

	conn_subfacet(new_subfacet) = new_node;
	}

	/// update facet and element connectivity
	Array<std::vector<Element>> & f_to_subfacet_double =
	mesh_facets.getData<std::vector<Element>>("facets_to_subfacet_double",
	type_subfacet, gt_subfacet);

	Array<std::vector<Element>> & el_to_subfacet_double =
	mesh_facets.getData<std::vector<Element>>(
	"elements_to_subfacet_double", type_subfacet, gt_subfacet);

	Array<std::vector<Element>> * sf_to_subfacet_double = NULL;

	if (spatial_dimension == 3)
	sf_to_subfacet_double = &mesh_facets.getData<std::vector<Element>>(
	"subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet);

	for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) {
	UInt old_node = doubled_nodes(old_nb_doubled_nodes + sf, 0);
	UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1);

	updateElementalConnectivity(mesh, old_node, new_node,
	el_to_subfacet_double(sf),
	&f_to_subfacet_double(sf));

	updateElementalConnectivity(mesh_facets, old_node, new_node,
	f_to_subfacet_double(sf));

	if (spatial_dimension == 3)
	updateElementalConnectivity(mesh_facets, old_node, new_node,
	(*sf_to_subfacet_double)(sf));
	}

	if (spatial_dimension == 2) {
	updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, true);
	updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, true);
	} else if (spatial_dimension == 3) {
	updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, false);
	updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, false);
	}
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::flipFacets(
	Mesh & mesh_facets, const ElementTypeMapArray<UInt> & global_connectivity,
	GhostType gt_facet) {
	AKANTU_DEBUG_IN();

	UInt spatial_dimension = mesh_facets.getSpatialDimension();

	/// get global connectivity for local mesh
	ElementTypeMapArray<UInt> global_connectivity_tmp("global_connectivity_tmp",
	mesh_facets.getID(),
	mesh_facets.getMemoryID());

	mesh_facets.initElementTypeMapArray(global_connectivity_tmp, 1,
	spatial_dimension - 1, gt_facet, true,
	_ek_regular, true);

	mesh_facets.getGlobalConnectivity(global_connectivity_tmp,
	spatial_dimension - 1, gt_facet);

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

	/// loop on every facet
	for (; it != end; ++it) {
	ElementType type_facet = *it;

	Array<UInt> & connectivity =
	mesh_facets.getConnectivity(type_facet, gt_facet);
	const Array<UInt> & g_connectivity =
	global_connectivity(type_facet, gt_facet);

	Array<std::vector<Element>> & el_to_f =
	mesh_facets.getElementToSubelement(type_facet, gt_facet);
	Array<Element> & subfacet_to_facet =
	mesh_facets.getSubelementToElement(type_facet, gt_facet);

	UInt nb_subfacet_per_facet = subfacet_to_facet.getNbComponent();
	UInt nb_nodes_per_facet = connectivity.getNbComponent();
	UInt nb_facet = connectivity.getSize();

	UInt nb_nodes_per_P1_facet =
	Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type_facet));

	Array<UInt> & global_conn_tmp =
	global_connectivity_tmp(type_facet, gt_facet);

	Array<UInt>::iterator<Vector<UInt>> conn_it =
	connectivity.begin(nb_nodes_per_facet);
	Array<UInt>::iterator<Vector<UInt>> gconn_tmp_it =
	global_conn_tmp.begin(nb_nodes_per_facet);
	Array<UInt>::const_iterator<Vector<UInt>> conn_glob_it =
	g_connectivity.begin(nb_nodes_per_facet);
	Array<Element>::iterator<Vector<Element>> subf_to_f =
	subfacet_to_facet.begin(nb_subfacet_per_facet);

	UInt * conn_tmp_pointer = new UInt[nb_nodes_per_facet];
	Vector<UInt> conn_tmp(conn_tmp_pointer, nb_nodes_per_facet);

	Element el_tmp;
	Element * subf_tmp_pointer = new Element[nb_subfacet_per_facet];
	Vector<Element> subf_tmp(subf_tmp_pointer, nb_subfacet_per_facet);

	for (UInt f = 0; f < nb_facet;
	++f, ++conn_it, ++gconn_tmp_it, ++subf_to_f, ++conn_glob_it) {

	Vector<UInt> & gconn_tmp = *gconn_tmp_it;
	const Vector<UInt> & conn_glob = *conn_glob_it;
	Vector<UInt> & conn_local = *conn_it;

	/// skip facet if connectivities are the same
	if (gconn_tmp == conn_glob)
	continue;

	/// re-arrange connectivity
	conn_tmp = conn_local;

	UInt * begin = conn_glob.storage();
	UInt * end = conn_glob.storage() + nb_nodes_per_facet;

	for (UInt n = 0; n < nb_nodes_per_facet; ++n) {

	UInt * new_node = std::find(begin, end, gconn_tmp(n));
	AKANTU_DEBUG_ASSERT(new_node != end, "Node not found");

	UInt new_position = new_node - begin;

	conn_local(new_position) = conn_tmp(n);
	}

	/// if 3D, check if facets are just rotated
	if (spatial_dimension == 3) {
	/// find first node
	UInt * new_node = std::find(begin, end, gconn_tmp(0));
	AKANTU_DEBUG_ASSERT(new_node != end, "Node not found");

	UInt new_position = new_node - begin;
	UInt n = 1;

	/// count how many nodes in the received connectivity follow
	/// the same order of those in the local connectivity
	for (; n < nb_nodes_per_P1_facet &&
	gconn_tmp(n) ==
	conn_glob((new_position + n) % nb_nodes_per_P1_facet);
	++n)
	;

	/// skip the facet inversion if facet is just rotated
	if (n == nb_nodes_per_P1_facet)
	continue;
	}

	/// update data to invert facet
	el_tmp = el_to_f(f)[0];
	el_to_f(f)[0] = el_to_f(f)[1];
	el_to_f(f)[1] = el_tmp;

	subf_tmp = (*subf_to_f);
	(*subf_to_f)(0) = subf_tmp(1);
	(*subf_to_f)(1) = subf_tmp(0);
	}

	delete[] subf_tmp_pointer;
	delete[] conn_tmp_pointer;
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::fillElementToSubElementsData(Mesh & mesh) {
	AKANTU_DEBUG_IN();

	if (mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) {
	AKANTU_DEBUG_INFO("There are not facets, add them in the mesh file or call "
	"the buildFacet method.");
	return;
	}

	UInt spatial_dimension = mesh.getSpatialDimension();
	ElementTypeMapArray<Real> barycenters("barycenter_tmp", mesh.getID(),
	mesh.getMemoryID());
	mesh.initElementTypeMapArray(barycenters, spatial_dimension, _all_dimensions);

	for (ghost_type_t::iterator gt = ghost_type_t::begin();
	gt != ghost_type_t::end(); ++gt) {
	Mesh::type_iterator it = mesh.firstType(_all_dimensions, *gt);
	Mesh::type_iterator end = mesh.lastType(_all_dimensions, *gt);
	for (; it != end; ++it) {
	UInt nb_element = mesh.getNbElement(it, gt);
	Array<Real> & barycenters_arr = barycenters(it, gt);
	barycenters_arr.resize(nb_element);
	Array<Real>::vector_iterator bary =
	barycenters_arr.begin(spatial_dimension);
	Array<Real>::vector_iterator bary_end =
	barycenters_arr.end(spatial_dimension);

	for (UInt el = 0; bary != bary_end; ++bary, ++el) {
	mesh.getBarycenter(el, it, bary->storage(), gt);
	}
	}
	}

	MeshAccessor mesh_accessor(mesh);
	for (Int sp(spatial_dimension); sp >= 1; --sp) {
	if (mesh.getNbElement(sp) == 0)
	continue;

	for (ghost_type_t::iterator git = ghost_type_t::begin();
	git != ghost_type_t::end(); ++git) {
	Mesh::type_iterator tit = mesh.firstType(sp, *git);
	Mesh::type_iterator tend = mesh.lastType(sp, *git);
	for (; tit != tend; ++tit) {
	mesh_accessor.getSubelementToElement(tit, git)
	.resize(mesh.getNbElement(tit, git));
	mesh_accessor.getSubelementToElement(tit, git).clear();
	}

	tit = mesh.firstType(sp - 1, *git);
	tend = mesh.lastType(sp - 1, *git);
	for (; tit != tend; ++tit) {
	mesh_accessor.getElementToSubelement(tit, git)
	.resize(mesh.getNbElement(tit, git));
	mesh.getElementToSubelement(tit, git).clear();
	}
	}

	CSR<Element> nodes_to_elements;
	buildNode2Elements(mesh, nodes_to_elements, sp);

	Element facet_element;

	for (ghost_type_t::iterator git = ghost_type_t::begin();
	git != ghost_type_t::end(); ++git) {
	Mesh::type_iterator tit = mesh.firstType(sp - 1, *git);
	Mesh::type_iterator tend = mesh.lastType(sp - 1, *git);

	facet_element.ghost_type = *git;
	for (; tit != tend; ++tit) {
	facet_element.type = *tit;

	Array<std::vector<Element>> & element_to_subelement =
	mesh.getElementToSubelement(tit, git);

	const Array<UInt> & connectivity = mesh.getConnectivity(tit, git);

	Array<UInt>::const_iterator<Vector<UInt>> fit =
	connectivity.begin(mesh.getNbNodesPerElement(*tit));
	Array<UInt>::const_iterator<Vector<UInt>> fend =
	connectivity.end(mesh.getNbNodesPerElement(*tit));

	UInt fid = 0;
	for (; fit != fend; ++fit, ++fid) {
	const Vector<UInt> & facet = *fit;
	facet_element.element = fid;
	std::map<Element, UInt> element_seen_counter;
	UInt nb_nodes_per_facet =
	mesh.getNbNodesPerElement(Mesh::getP1ElementType(*tit));
	for (UInt n(0); n < nb_nodes_per_facet; ++n) {
	CSR<Element>::iterator eit = nodes_to_elements.begin(facet(n));
	CSR<Element>::iterator eend = nodes_to_elements.end(facet(n));
	for (; eit != eend; ++eit) {
	Element & elem = *eit;
	std::map<Element, UInt>::iterator cit =
	element_seen_counter.find(elem);
	if (cit != element_seen_counter.end()) {
	cit->second++;
	} else {
	element_seen_counter[elem] = 1;
	}
	}
	}

	std::vector<Element> connected_elements;
	std::map<Element, UInt>::iterator cit = element_seen_counter.begin();
	std::map<Element, UInt>::iterator cend = element_seen_counter.end();
	for (; cit != cend; ++cit) {
	if (cit->second == nb_nodes_per_facet)
	connected_elements.push_back(cit->first);
	}

	std::vector<Element>::iterator ceit = connected_elements.begin();
	std::vector<Element>::iterator ceend = connected_elements.end();
	for (; ceit != ceend; ++ceit)
	element_to_subelement(fid).push_back(*ceit);

	for (UInt ce = 0; ce < connected_elements.size(); ++ce) {
	Element & elem = connected_elements[ce];
	Array<Element> & subelement_to_element =
	mesh_accessor.getSubelementToElement(elem.type,
	elem.ghost_type);

	UInt f(0);
	for (; f < mesh.getNbFacetsPerElement(elem.type) &&
	subelement_to_element(elem.element, f) != ElementNull;
	++f)
	;

	AKANTU_DEBUG_ASSERT(
	f < mesh.getNbFacetsPerElement(elem.type),
	"The element " << elem << " seems to have too many facets!! ("
	<< f << " < "
	<< mesh.getNbFacetsPerElement(elem.type) << ")");

	subelement_to_element(elem.element, f) = facet_element;
	}
	}
	}
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	template <bool third_dim_points>
	bool MeshUtils::findElementsAroundSubfacet(
	const Mesh & mesh, const Mesh & mesh_facets,
	const Element & starting_element, const Element & end_facet,
	const Vector<UInt> & subfacet_connectivity,
	std::vector<Element> & elem_list, std::vector<Element> & facet_list,
	std::vector<Element> * subfacet_list) {
	AKANTU_DEBUG_IN();

	/// preallocated stuff before starting
	bool facet_matched = false;

	elem_list.clear();
	facet_list.clear();

	if (third_dim_points)
	subfacet_list->clear();

	elem_list.push_back(starting_element);

	const Array<UInt> * facet_connectivity = NULL;
	const Array<UInt> * sf_connectivity = NULL;
	const Array<Element> * facet_to_element = NULL;
	const Array<Element> * subfacet_to_facet = NULL;

	ElementType current_type = _not_defined;
	GhostType current_ghost_type = _casper;

	ElementType current_facet_type = _not_defined;
	GhostType current_facet_ghost_type = _casper;

	ElementType current_subfacet_type = _not_defined;
	GhostType current_subfacet_ghost_type = _casper;

	const Array<std::vector<Element>> * element_to_facet = NULL;

	const Element * opposing_el = NULL;

	std::queue<Element> elements_to_check;
	elements_to_check.push(starting_element);

	/// keep going until there are elements to check
	while (!elements_to_check.empty()) {

	/// check current element
	Element & current_el = elements_to_check.front();

	if (current_el.type != current_type \|\|
	current_el.ghost_type != current_ghost_type) {

	current_type = current_el.type;
	current_ghost_type = current_el.ghost_type;

	facet_to_element =
	&mesh_facets.getSubelementToElement(current_type, current_ghost_type);
	}

	/// loop over each facet of the element
	for (UInt f = 0; f < facet_to_element->getNbComponent(); ++f) {

	const Element & current_facet =
	(*facet_to_element)(current_el.element, f);

	if (current_facet == ElementNull)
	continue;

	if (current_facet_type != current_facet.type \|\|
	current_facet_ghost_type != current_facet.ghost_type) {

	current_facet_type = current_facet.type;
	current_facet_ghost_type = current_facet.ghost_type;

	element_to_facet = &mesh_facets.getElementToSubelement(
	current_facet_type, current_facet_ghost_type);
	facet_connectivity = &mesh_facets.getConnectivity(
	current_facet_type, current_facet_ghost_type);

	if (third_dim_points)
	subfacet_to_facet = &mesh_facets.getSubelementToElement(
	current_facet_type, current_facet_ghost_type);
	}

	/// check if end facet is reached
	if (current_facet == end_facet)
	facet_matched = true;

	/// add this facet if not already passed
	if (std::find(facet_list.begin(), facet_list.end(), current_facet) ==
	facet_list.end() &&
	hasElement(*facet_connectivity, current_facet,
	subfacet_connectivity)) {
	facet_list.push_back(current_facet);

	if (third_dim_points) {
	/// check subfacets
	for (UInt sf = 0; sf < subfacet_to_facet->getNbComponent(); ++sf) {
	const Element & current_subfacet =
	(*subfacet_to_facet)(current_facet.element, sf);

	if (current_subfacet == ElementNull)
	continue;

	if (current_subfacet_type != current_subfacet.type \|\|
	current_subfacet_ghost_type != current_subfacet.ghost_type) {
	current_subfacet_type = current_subfacet.type;
	current_subfacet_ghost_type = current_subfacet.ghost_type;

	sf_connectivity = &mesh_facets.getConnectivity(
	current_subfacet_type, current_subfacet_ghost_type);
	}

	if (std::find(subfacet_list->begin(), subfacet_list->end(),
	current_subfacet) == subfacet_list->end() &&
	hasElement(*sf_connectivity, current_subfacet,
	subfacet_connectivity))
	subfacet_list->push_back(current_subfacet);
	}
	}
	} else
	continue;

	/// consider opposing element
	if ((*element_to_facet)(current_facet.element)[0] == current_el)
	opposing_el = &(*element_to_facet)(current_facet.element)[1];
	else
	opposing_el = &(*element_to_facet)(current_facet.element)[0];

	/// skip null elements since they are on a boundary
	if (*opposing_el == ElementNull)
	continue;

	/// skip this element if already added
	if (std::find(elem_list.begin(), elem_list.end(), *opposing_el) !=
	elem_list.end())
	continue;

	/// only regular elements have to be checked
	if (opposing_el->kind == _ek_regular)
	elements_to_check.push(*opposing_el);

	elem_list.push_back(*opposing_el);

	#ifndef AKANTU_NDEBUG
	const Array<UInt> & conn_elem =
	mesh.getConnectivity(opposing_el->type, opposing_el->ghost_type);

	AKANTU_DEBUG_ASSERT(
	hasElement(conn_elem, *opposing_el, subfacet_connectivity),
	"Subfacet doesn't belong to this element");
	#endif
	}

	/// erased checked element from the list
	elements_to_check.pop();
	}

	AKANTU_DEBUG_OUT();
	return facet_matched;
	}

	/* -------------------------------------------------------------------------- */
	void MeshUtils::buildSegmentToNodeType(const Mesh & mesh, Mesh & mesh_facets) {
	buildAllFacets(mesh, mesh_facets, 1);

	UInt spatial_dimension = mesh.getSpatialDimension();

	const ElementTypeMapArray<UInt> & element_to_rank =
	mesh.getElementSynchronizer().getPrankToElement();
	Int local_rank = StaticCommunicator::getStaticCommunicator().whoAmI();

	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(1, ghost_type);
	Mesh::type_iterator end = mesh_facets.lastType(1, ghost_type);
	for (; it != end; ++it) {
	ElementType type = *it;
	UInt nb_segments = mesh_facets.getNbElement(type, ghost_type);

	// allocate the data
	Array<Int> & segment_to_nodetype = *(mesh_facets.getDataPointer<Int>(
	"segment_to_nodetype", type, ghost_type));

	std::set<Element> connected_elements;
	const Array<std::vector<Element>> & segment_to_2Delement =
	mesh_facets.getElementToSubelement(type, ghost_type);

	// loop over segments
	for (UInt s = 0; s < nb_segments; ++s) {

	// determine the elements connected to the segment
	connected_elements.clear();
	const std::vector<Element> & twoD_elements = segment_to_2Delement(s);

	if (spatial_dimension == 2) {
	// if 2D just take the elements connected to the segments
	connected_elements.insert(twoD_elements.begin(), twoD_elements.end());
	} else if (spatial_dimension == 3) {
	// if 3D a second loop is needed to get to the 3D elements
	std::vector<Element>::const_iterator facet = twoD_elements.begin();

	for (; facet != twoD_elements.end(); ++facet) {
	const std::vector<Element> & threeD_elements =
	mesh_facets.getElementToSubelement(
	facet->type, facet->ghost_type)(facet->element);
	connected_elements.insert(threeD_elements.begin(),
	threeD_elements.end());
	}
	}

	// get the minimum processor rank associated to the connected
	// elements and verify if ghost and not ghost elements are
	// found
	Int minimum_rank = std::numeric_limits<Int>::max();

	// two booleans saying if not ghost and ghost elements are found in the
	// loop
	bool ghost_found[2];
	ghost_found[0] = false;
	ghost_found[1] = false;

	std::set<Element>::iterator connected_elements_it =
	connected_elements.begin();

	for (; connected_elements_it != connected_elements.end();
	++connected_elements_it) {
	if (*connected_elements_it == ElementNull)
	continue;
	ghost_found[connected_elements_it->ghost_type] = true;
	const Array<UInt> & el_to_rank_array = element_to_rank(
	connected_elements_it->type, connected_elements_it->ghost_type);
	minimum_rank =
	std::min(minimum_rank,
	Int(el_to_rank_array(connected_elements_it->element)));
	}

	// if no ghost elements are found the segment is local
	if (!ghost_found[1])
	segment_to_nodetype(s) = -1;
	// if no not ghost elements are found the segment is pure ghost
	else if (!ghost_found[0])
	segment_to_nodetype(s) = -3;
	// if the minimum rank is equal to the local rank, the segment is master
	else if (local_rank == minimum_rank)
	segment_to_nodetype(s) = -2;
	// if the minimum rank is less than the local rank, the segment is slave
	else if (local_rank > minimum_rank)
	segment_to_nodetype(s) = minimum_rank;
	else
	AKANTU_DEBUG_ERROR("The local rank cannot be smaller than the "
	"minimum rank if both ghost and not ghost "
	"elements are found");
	}
	}
	}
	}

	/* -------------------------------------------------------------------------- */
	UInt MeshUtils::updateLocalMasterGlobalConnectivity(Mesh & mesh,
	UInt local_nb_new_nodes) {
	StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
	Int rank = comm.whoAmI();
	Int nb_proc = comm.getNbProc();
	if (nb_proc == 1)
	return local_nb_new_nodes;

	/// resize global ids array
	Array<UInt> & nodes_global_ids = mesh.getGlobalNodesIds();
	UInt old_nb_nodes = mesh.getNbNodes() - local_nb_new_nodes;

	nodes_global_ids.resize(mesh.getNbNodes());

	/// compute the number of global nodes based on the number of old nodes
	Vector<UInt> old_local_master_nodes(nb_proc);
	for (UInt n = 0; n < old_nb_nodes; ++n)
	if (mesh.isLocalOrMasterNode(n))
	++old_local_master_nodes(rank);
	comm.allGather(old_local_master_nodes);
	UInt old_global_nodes =
	std::accumulate(old_local_master_nodes.storage(),
	old_local_master_nodes.storage() + nb_proc, 0);

	/// compute amount of local or master doubled nodes
	Vector<UInt> local_master_nodes(nb_proc);
	for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n)
	if (mesh.isLocalOrMasterNode(n))
	++local_master_nodes(rank);

	comm.allGather(local_master_nodes);

	/// update global number of nodes
	UInt total_nb_new_nodes = std::accumulate(
	local_master_nodes.storage(), local_master_nodes.storage() + nb_proc, 0);

	if (total_nb_new_nodes == 0)
	return 0;

	/// set global ids of local and master nodes
	UInt starting_index =
	std::accumulate(local_master_nodes.storage(),
	local_master_nodes.storage() + rank, old_global_nodes);

	for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) {
	if (mesh.isLocalOrMasterNode(n)) {
	nodes_global_ids(n) = starting_index;
	++starting_index;
	}
	}

	MeshAccessor mesh_accessor(mesh);
	mesh_accessor.setNbGlobalNodes(old_global_nodes + total_nb_new_nodes);
	return total_nb_new_nodes;
	}

	/* -------------------------------------------------------------------------- */
	// Deactivating -Wunused-parameter
	#if defined(__INTEL_COMPILER)
	//#pragma warning ( disable : 383 )
	#elif defined(__clang__) // test clang to be sure that when we test for gnu it
	// is only gnu
	#elif (defined(__GNUC__) \|\| defined(__GNUG__))
	#define GCC_VERSION \
	(__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
	#if GCC_VERSION > 40600
	#pragma GCC diagnostic push
	#endif
	#pragma GCC diagnostic ignored "-Wunused-parameter"
	#endif

	/* -------------------------------------------------------------------------- */
	void MeshUtils::updateElementalConnectivity(
	Mesh & mesh, UInt old_node, UInt new_node,
	const std::vector<Element> & element_list,
	__attribute__((unused)) const std::vector<Element> * facet_list) {
	AKANTU_DEBUG_IN();

	ElementType el_type = _not_defined;
	GhostType gt_type = _casper;
	Array<UInt> * conn_elem = NULL;
	#if defined(AKANTU_COHESIVE_ELEMENT)
	const Array<Element> * cohesive_facets = NULL;
	#endif
	UInt nb_nodes_per_element = 0;
	UInt * n_update = NULL;

	for (UInt el = 0; el < element_list.size(); ++el) {
	const Element & elem = element_list[el];
	if (elem.type == _not_defined)
	continue;

	if (elem.type != el_type \|\| elem.ghost_type != gt_type) {
	el_type = elem.type;
	gt_type = elem.ghost_type;
	conn_elem = &mesh.getConnectivity(el_type, gt_type);
	nb_nodes_per_element = conn_elem->getNbComponent();
	#if defined(AKANTU_COHESIVE_ELEMENT)
	if (elem.kind == _ek_cohesive)
	cohesive_facets =
	&mesh.getMeshFacets().getSubelementToElement(el_type, gt_type);
	#endif
	}

	#if defined(AKANTU_COHESIVE_ELEMENT)
	if (elem.kind == _ek_cohesive) {

	AKANTU_DEBUG_ASSERT(
	facet_list != NULL,
	"Provide a facet list in order to update cohesive elements");

	/// loop over cohesive element's facets
	for (UInt f = 0, n = 0; f < 2; ++f, n += nb_nodes_per_element / 2) {
	const Element & facet = (*cohesive_facets)(elem.element, f);

	/// skip facets if not present in the list
	if (std::find(facet_list->begin(), facet_list->end(), facet) ==
	facet_list->end())
	continue;

	n_update = std::find(
	conn_elem->storage() + elem.element * nb_nodes_per_element + n,
	conn_elem->storage() + elem.element * nb_nodes_per_element + n +
	nb_nodes_per_element / 2,
	old_node);

	AKANTU_DEBUG_ASSERT(n_update !=
	conn_elem->storage() +
	elem.element * nb_nodes_per_element + n +
	nb_nodes_per_element / 2,
	"Node not found in current element");

	/// update connectivity
	*n_update = new_node;
	}
	} else {
	#endif
	n_update =
	std::find(conn_elem->storage() + elem.element * nb_nodes_per_element,
	conn_elem->storage() + elem.element * nb_nodes_per_element +
	nb_nodes_per_element,
	old_node);

	AKANTU_DEBUG_ASSERT(n_update != conn_elem->storage() +
	elem.element * nb_nodes_per_element +
	nb_nodes_per_element,
	"Node not found in current element");

	/// update connectivity
	*n_update = new_node;
	#if defined(AKANTU_COHESIVE_ELEMENT)
	}
	#endif
	}

	AKANTU_DEBUG_OUT();
	}

	// Reactivating -Wunused-parameter
	#if defined(__INTEL_COMPILER)
	//#pragma warning ( disable : 383 )
	#elif defined(__clang__) // test clang to be sure that when we test for gnu it
	// is only gnu
	#elif defined(__GNUG__)
	#if GCC_VERSION > 40600
	#pragma GCC diagnostic pop
	#else
	#pragma GCC diagnostic warning "-Wunused-parameter"
	#endif
	#endif
	/* -------------------------------------------------------------------------- */

	} // akantu

	// LocalWords: ElementType

mesh_utils.ccNo OneTemporaryActions

File Metadata

mesh_utils.ccView Options

Event Timeline

mesh_utils.cc
No OneTemporary
Actions

mesh_utils.cc
View Options