Diffusion akantu (bugfix/69-segmentation-fault-caused-by-model-applybc-aka-fromtraction-trac-traction) Browse

69-segmentation…

cohesive_element_inserter_helper.cc
View Options

			/**
			* Copyright (©) 2020-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
			* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
			*
			* This file is part of Akantu
			*
			* 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_helper.hh"
			/* -------------------------------------------------------------------------- */
			#include "data_accessor.hh"
			#include "element_synchronizer.hh"
			#include "fe_engine.hh"
			#include "mesh_accessor.hh"
			/* -------------------------------------------------------------------------- */
			#include <queue>
			/* -------------------------------------------------------------------------- */

			namespace akantu {

			/* -------------------------------------------------------------------------- */
			CohesiveElementInserterHelper::CohesiveElementInserterHelper(
			Mesh & mesh, const ElementTypeMapArray<bool> & facet_insertion)
			: doubled_nodes(0, 2, "doubled_nodes"), mesh(mesh),
			mesh_facets(mesh.getMeshFacets()) {

			auto spatial_dimension = mesh_facets.getSpatialDimension();

			for (auto gt : ghost_types) {
			for (auto type : mesh_facets.elementTypes(_ghost_type = gt)) {
			nb_new_facets(type, gt) = mesh_facets.getNbElement(type, gt);
			}
			}

			std::array<Int, 2> nb_facet_to_insert{0, 0};
			// creates a vector of the facets to insert
			std::vector<Element> potential_facets_to_double;
			for (auto gt_facet : ghost_types) {
			for (auto type_facet :
			mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
			const auto & f_insertion = facet_insertion(type_facet, gt_facet);
			auto & counter = nb_facet_to_insert[gt_facet == _not_ghost ? 0 : 1];

			for (auto && data : enumerate(f_insertion)) {
			if (std::get<1>(data)) {
			auto el = std::get<0>(data);
			potential_facets_to_double.push_back({type_facet, el, gt_facet});
			++counter;
			}
			}
			}
			}

			// Defines a global order of insertion of cohesive element to ensure node
			// doubling happens in the same order, this is necessary for the global node
			// numbering
			if (mesh_facets.isDistributed()) {
			const auto & comm = mesh_facets.getCommunicator();

			ElementTypeMapArray<Int> global_orderings;

			global_orderings.initialize(mesh_facets,
			_spatial_dimension = spatial_dimension - 1,
			_with_nb_element = true, _default_value = -1);

			auto starting_index = nb_facet_to_insert[0];
			comm.exclusiveScan(starting_index);

			// define the local numbers for facet to insert
			for (auto gt_facet : ghost_types) {
			for (auto type_facet :
			mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
			for (auto data : zip(facet_insertion(type_facet, gt_facet),
			global_orderings(type_facet, gt_facet))) {
			if (std::get<0>(data)) {
			std::get<1>(data) = starting_index;
			++starting_index;
			}
			}
			}
			}

			// retrieves the order number from neighboring processors
			auto && synchronizer = mesh_facets.getElementSynchronizer();
			SimpleElementDataAccessor<Int> data_accessor(
			global_orderings, SynchronizationTag::_ce_insertion_order);
			synchronizer.synchronizeOnce(data_accessor,
			SynchronizationTag::_ce_insertion_order);

			// sort the facets to double based on this global ordering
			std::sort(potential_facets_to_double.begin(),
			potential_facets_to_double.end(),
			[&global_orderings](auto && el1, auto && el2) {
			return global_orderings(el1) < global_orderings(el2);
			});
			}

			for (auto d : arange(spatial_dimension)) {
			facets_to_double_by_dim[d] = std::make_unique<Array<Element>>(
			0, 2, "facets_to_double_" + std::to_string(d));
			}

			auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];

			MeshAccessor mesh_accessor(mesh_facets);
			auto & elements_to_subelements = mesh_accessor.getElementToSubelement();

			for (auto && facet_to_double : potential_facets_to_double) {
			auto gt_facet = facet_to_double.ghost_type;
			auto type_facet = facet_to_double.type;
			auto & elements_to_facets = elements_to_subelements(type_facet, gt_facet);
			auto & elements_to_facet = elements_to_facets(facet_to_double.element);

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

			auto new_facet = nb_new_facets(type_facet, gt_facet)++;

			facets_to_double.push_back(Vector<Element>{
			{facet_to_double, Element{type_facet, new_facet, gt_facet}}});

			/// update facet_to_element vector
			auto & element_to_update = elements_to_facet[1];
			auto el = element_to_update.element;

			auto & facets_to_elements = mesh_accessor.getSubelementToElement(
			element_to_update.type, element_to_update.ghost_type);
			auto facets_to_element =
			make_view(facets_to_elements, facets_to_elements.getNbComponent())
			.begin()[el];
			auto facet_to_update = std::find(facets_to_element.begin(),
			facets_to_element.end(), facet_to_double);

			AKANTU_DEBUG_ASSERT(facet_to_update != facets_to_element.end(),
			"Facet not found");

			facet_to_update->element = new_facet;

			/// update elements connected to facet
			const auto & first_facet_list = elements_to_facet;
			elements_to_facets.push_back(first_facet_list);

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

			elements_to_facets(facet_to_double.element)[1] = ElementNull;
			}
			}

			/* -------------------------------------------------------------------------- */
			inline auto
			CohesiveElementInserterHelper::hasElement(const Vector<Idx> & nodes_element,
			const Vector<Idx> & nodes) -> bool {
			// if one of the nodes of "nodes" is not in "nodes_element" it stops
			auto it = std::mismatch(nodes.begin(), nodes.end(), nodes_element.begin(),
			[&](auto && node, auto && /node2/) -> bool {
			auto it = std::find(nodes_element.begin(),
			nodes_element.end(), node);
			return (it != nodes_element.end());
			});

			// true if all "nodes" where found in "nodes_element"
			return (it.first == nodes.end());
			}

			/* -------------------------------------------------------------------------- */
			inline auto CohesiveElementInserterHelper::removeElementsInVector(
			const std::vector<Element> & elem_to_remove,
			std::vector<Element> & elem_list) -> bool {
			if (elem_list.size() <= elem_to_remove.size()) {
			return true;
			}

			auto el_it = elem_to_remove.begin();
			auto el_last = elem_to_remove.end();
			std::vector<Element>::iterator el_del;

			UInt deletions = 0;

			for (; el_it != el_last; ++el_it) {
			el_del = std::find(elem_list.begin(), elem_list.end(), *el_it);

			if (el_del != elem_list.end()) {
			elem_list.erase(el_del);
			++deletions;
			}
			}

			AKANTU_DEBUG_ASSERT(deletions == 0 \|\| deletions == elem_to_remove.size(),
			"Not all elements have been erased");

			return deletions == 0;
			}

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

			MeshAccessor mesh_accessor(mesh);
			for (const auto & element : element_list) {
			if (element.type == _not_defined) {
			continue;
			}

			auto && connectivity = mesh_accessor.getConnectivity(element);

			if (element.kind() == _ek_cohesive) {
			AKANTU_DEBUG_ASSERT(
			facet_list != nullptr,
			"Provide a facet list in order to update cohesive elements");

			auto && facets = mesh_facets.getSubelementToElement(element);

			auto facet_nb_nodes = connectivity.size() / 2;

			/// loop over cohesive element's facets
			for (const auto & facet : enumerate(facets)) {
			/// skip facets if not present in the list
			if (std::find(facet_list->begin(), facet_list->end(),
			std::get<1>(facet)) == facet_list->end()) {
			continue;
			}

			auto n = std::get<0>(facet);

			auto begin =
			connectivity.begin() + static_cast<UInt>(n * facet_nb_nodes);
			auto end = begin + facet_nb_nodes;

			auto it = std::find(begin, end, old_node);
			AKANTU_DEBUG_ASSERT(it != end, "Node not found in current element");

			*it = new_node;
			}
			} else {
			auto it = std::find(connectivity.begin(), connectivity.end(), old_node);
			AKANTU_DEBUG_ASSERT(it != connectivity.end(),
			"Node not found in current element");

			/// update connectivity
			*it = new_node;
			}
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void CohesiveElementInserterHelper::updateSubelementToElement(Int dim,
			bool facet_mode) {
			auto & facets_to_double = *facets_to_double_by_dim[dim];
			auto & facets_to_subfacets = elementsOfDimToElementsOfDim(
			dim + static_cast<decltype(dim)>(facet_mode), dim);

			MeshAccessor mesh_accessor(mesh_facets);

			for (auto && data :
			zip(make_view(facets_to_double, 2), facets_to_subfacets)) {
			const auto & old_subfacet = std::get<0>(data)[0];
			const auto & new_subfacet = std::get<0>(data)[1];
			auto & facet_to_subfacets = std::get<1>(data);

			for (auto & facet : facet_to_subfacets) {
			auto && subfacets = mesh_accessor.getSubelementToElement(facet);

			auto && subfacet_to_update_it =
			std::find(subfacets.begin(), subfacets.end(), old_subfacet);

			AKANTU_DEBUG_ASSERT(subfacet_to_update_it != subfacets.end(),
			"Subfacet not found");

			*subfacet_to_update_it = new_subfacet;
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			void CohesiveElementInserterHelper::updateElementToSubelement(Int dim,
			bool facet_mode) {
			auto & facets_to_double = *facets_to_double_by_dim[dim];
			auto & facets_to_subfacets = elementsOfDimToElementsOfDim(
			dim + static_cast<decltype(dim)>(facet_mode), dim);

			MeshAccessor mesh_accessor(mesh_facets);
			// resize the arrays
			mesh_accessor.getElementToSubelement().initialize(
			mesh_facets, _spatial_dimension = dim, _with_nb_element = true);

			for (auto && data :
			zip(make_view(facets_to_double, 2), facets_to_subfacets)) {
			const auto & new_facet = std::get<0>(data)[1];
			mesh_accessor.getElementToSubelement(new_facet) = std::get<1>(data);
			}
			}

			/* -------------------------------------------------------------------------- */
			void CohesiveElementInserterHelper::updateQuadraticSegments(Int dim) {
			AKANTU_DEBUG_IN();
			auto spatial_dimension = mesh.getSpatialDimension();
			auto & facets_to_double = *facets_to_double_by_dim[dim];

			MeshAccessor mesh_accessor(mesh_facets);
			auto & connectivities = mesh_accessor.getConnectivities();

			/// this ones matter only for segments in 3D
			Array<std::vector<Element>> * element_to_subfacet_double = nullptr;
			Array<std::vector<Element>> * facet_to_subfacet_double = nullptr;

			if (dim == spatial_dimension - 2) {
			element_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 2, dim);
			facet_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 1, dim);
			}

			const auto & element_to_subelement = mesh_facets.getElementToSubelement();
			std::vector<Int> middle_nodes;

			for (auto && facet_to_double : make_view(facets_to_double, 2)) {
			auto & old_facet = facet_to_double[0];
			if (old_facet.type != _segment_3) {
			continue;
			}

			auto node = connectivities(old_facet, 2);
			if (not mesh.isPureGhostNode(node)) {
			middle_nodes.push_back(node);
			}
			}

			auto n = doubled_nodes.size();
			doubleNodes(middle_nodes);

			for (auto && data : enumerate(make_view(facets_to_double, 2))) {
			auto facet = std::get<0>(data);
			auto & old_facet = std::get<1>(data)[0];
			if (old_facet.type != _segment_3) {
			continue;
			}

			auto old_node = connectivities(old_facet, 2);

			if (mesh.isPureGhostNode(old_node)) {
			continue;
			}

			auto new_node = doubled_nodes(n, 1);

			auto & new_facet = std::get<1>(data)[1];
			connectivities(new_facet, 2) = new_node;

			if (dim == spatial_dimension - 2) {
			updateElementalConnectivity(mesh_facets, old_node, new_node,
			element_to_subelement(new_facet, 0));

			updateElementalConnectivity(mesh, old_node, new_node,
			(*element_to_subfacet_double)(facet),
			&(*facet_to_subfacet_double)(facet));
			} else {
			updateElementalConnectivity(mesh, old_node, new_node,
			element_to_subelement(new_facet, 0));
			}

			++n;
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			Int CohesiveElementInserterHelper::insertCohesiveElement() {
			auto nb_new_facets = insertFacetsOnly();
			if (nb_new_facets == 0) {
			return 0;
			}

			updateCohesiveData();
			return nb_new_facets;
			}

			/* -------------------------------------------------------------------------- */
			Int CohesiveElementInserterHelper::insertFacetsOnly() {
			auto spatial_dimension = mesh.getSpatialDimension();

			if (facets_to_double_by_dim[spatial_dimension - 1]->empty()) {
			return 0;
			}

			if (spatial_dimension == 1) {
			doublePointFacet();
			} else if (spatial_dimension == 2) {
			doubleFacets<1>();
			findSubfacetToDouble<1>();

			doubleSubfacet<2>();

			} else if (spatial_dimension == 3) {
			doubleFacets<2>();
			findSubfacetToDouble<2>();

			doubleFacets<1>();
			findSubfacetToDouble<1>();

			doubleSubfacet<3>();
			}

			return facets_to_double_by_dim[spatial_dimension - 1]->size();
			}

			/* -------------------------------------------------------------------------- */
			template <Int dim> void CohesiveElementInserterHelper::doubleFacets() {
			AKANTU_DEBUG_IN();
			NewElementsEvent new_facets;
			auto spatial_dimension = mesh_facets.getSpatialDimension();

			auto & facets_to_double = *facets_to_double_by_dim[dim];
			MeshAccessor mesh_accessor(mesh_facets);
			for (auto && facet_to_double : make_view(facets_to_double, 2)) {
			auto && old_facet = facet_to_double[0];
			auto && new_facet = facet_to_double[1];

			auto & facets_connectivities =
			mesh_accessor.getConnectivity(old_facet.type, old_facet.ghost_type);
			facets_connectivities.resize(
			nb_new_facets(old_facet.type, old_facet.ghost_type));

			auto facets_connectivities_begin =
			make_view(facets_connectivities, facets_connectivities.getNbComponent())
			.begin();

			// copy the connectivities
			auto && new_conn(facets_connectivities_begin[new_facet.element]);
			auto && old_conn(facets_connectivities_begin[old_facet.element]);
			new_conn = old_conn;

			// this will fail if multiple facet types exists for a given element type
			// \TODO handle multiple sub-facet types
			auto nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(old_facet.type);

			auto & subfacets_to_facets = mesh_accessor.getSubelementToElementNC(
			old_facet.type, old_facet.ghost_type);

			subfacets_to_facets.resize(
			nb_new_facets(old_facet.type, old_facet.ghost_type), ElementNull);

			auto subfacets_to_facets_begin =
			make_view(subfacets_to_facets, nb_subfacet_per_facet).begin();

			// copy the subfacet to facets information
			auto && old_subfacets_to_facet(
			subfacets_to_facets_begin[old_facet.element]);
			auto && new_subfacet_to_facet(subfacets_to_facets_begin[new_facet.element]);
			new_subfacet_to_facet = old_subfacets_to_facet;

			for (auto & subfacet : old_subfacets_to_facet) {
			if (subfacet == ElementNull) {
			continue;
			}

			/// update facet_to_subfacet array
			mesh_accessor.getElementToSubelement(subfacet).push_back(new_facet);
			}

			new_facets.getList().push_back(new_facet);
			}

			/// update facet_to_subfacet and _segment_3 facets if any
			if (dim != spatial_dimension - 1) {
			updateSubelementToElement(dim, true);
			updateElementToSubelement(dim, true);
			updateQuadraticSegments(dim);
			} else {
			updateQuadraticSegments(dim);
			}

			mesh_facets.sendEvent(new_facets);
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template <Int dim> void CohesiveElementInserterHelper::findSubfacetToDouble() {
			AKANTU_DEBUG_IN();

			auto spatial_dimension = mesh_facets.getSpatialDimension();

			MeshAccessor mesh_accessor(mesh_facets);

			auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
			const auto & subfacets_to_facets = mesh_facets.getSubelementToElement();
			auto & elements_to_facets = mesh_accessor.getElementToSubelement();

			/// loop on every facet
			for (auto f : arange(2)) {
			for (auto && facet_to_double : make_view(facets_to_double, 2)) {
			auto old_facet = facet_to_double[f];
			auto new_facet = facet_to_double[1 - f];

			const auto & starting_element = elements_to_facets(new_facet, 0)[0];
			auto current_facet = old_facet;

			auto && subfacets_to_facet = subfacets_to_facets.get(old_facet);
			/// loop on every subfacet
			for (auto & subfacet : subfacets_to_facet) {
			if (subfacet == ElementNull) {
			continue;
			}

			if (dim == spatial_dimension - 2) {
			auto && subsubfacets_to_subfacet = subfacets_to_facets.get(subfacet);
			/// loop on every subsubfacet
			for (auto & subsubfacet : subsubfacets_to_subfacet) {
			if (subsubfacet == ElementNull) {
			continue;
			}

			auto && subsubfacet_connectivity =
			mesh_facets.getConnectivity(subsubfacet);

			std::vector<Element> element_list;
			std::vector<Element> facet_list;
			std::vector<Element> subfacet_list;

			/// check if subsubfacet is to be doubled
			if (!findElementsAroundSubfacet(
			starting_element, current_facet, subsubfacet_connectivity,
			element_list, facet_list, &subfacet_list) and
			!removeElementsInVector(subfacet_list,
			elements_to_facets(subsubfacet))) {
			Element new_subsubfacet{
			subsubfacet.type,
			nb_new_facets(subsubfacet.type, subsubfacet.ghost_type)++,
			subsubfacet.ghost_type};
			facets_to_double_by_dim[dim - 1]->push_back(
			Vector<Element>{subsubfacet, new_subsubfacet});
			elementsOfDimToElementsOfDim(dim - 1, dim - 1)
			.push_back(subfacet_list);
			elementsOfDimToElementsOfDim(dim, dim - 1).push_back(facet_list);
			elementsOfDimToElementsOfDim(dim + 1, dim - 1)
			.push_back(element_list);
			}
			}
			} else {
			std::vector<Element> element_list;
			std::vector<Element> facet_list;

			auto && subfacet_connectivity = mesh_facets.getConnectivity(subfacet);

			/// check if subfacet is to be doubled
			if (!findElementsAroundSubfacet(starting_element, current_facet,
			subfacet_connectivity, element_list,
			facet_list) and
			!removeElementsInVector(facet_list,
			elements_to_facets(subfacet))) {
			Element new_subfacet{
			subfacet.type,
			nb_new_facets(subfacet.type, subfacet.ghost_type)++,
			subfacet.ghost_type};
			facets_to_double_by_dim[dim - 1]->push_back(
			Vector<Element>{subfacet, new_subfacet});

			elementsOfDimToElementsOfDim(dim, dim - 1).push_back(facet_list);
			elementsOfDimToElementsOfDim(dim + 1, dim - 1)
			.push_back(element_list);
			}
			}
			}
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			void CohesiveElementInserterHelper::doubleNodes(
			const std::vector<Idx> & old_nodes) {
			auto & position = mesh.getNodes();
			auto spatial_dimension = mesh.getSpatialDimension();

			auto old_nb_nodes = position.size();
			position.reserve(old_nb_nodes + old_nodes.size());
			doubled_nodes.reserve(doubled_nodes.size() + old_nodes.size());

			auto position_begin = position.begin(spatial_dimension);
			for (auto && data : enumerate(old_nodes)) {
			auto n = std::get<0>(data);
			auto old_node = std::get<1>(data);
			decltype(old_node) new_node = old_nb_nodes + n;

			/// store doubled nodes
			doubled_nodes.push_back(Vector<Idx>{old_node, new_node});

			/// update position
			auto && coords = position_begin[old_node];
			position.push_back(coords);
			}
			}

			/* -------------------------------------------------------------------------- */
			bool CohesiveElementInserterHelper::findElementsAroundSubfacet(
			const Element & starting_element, const Element & end_facet,
			const Vector<Idx> & subfacet_connectivity,
			std::vector<Element> & element_list, std::vector<Element> & facet_list,
			std::vector<Element> * subfacet_list) {
			bool facet_matched = false;

			element_list.push_back(starting_element);

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

			/// keep going as long as there are elements to check
			while (not elements_to_check.empty()) {
			/// check current element
			Element & current_element = elements_to_check.front();

			auto && facets_to_element =
			mesh_facets.getSubelementToElement(current_element);

			// for every facet of the element
			for (const auto & current_facet : facets_to_element) {
			if (current_facet == ElementNull) {
			continue;
			}

			if (current_facet == end_facet) {
			facet_matched = true;
			}

			auto find_facet =
			std::find(facet_list.begin(), facet_list.end(), current_facet);

			// facet already listed or subfacet_connectivity is not in the
			// connectivity of current_facet;
			if ((find_facet != facet_list.end()) or
			not hasElement(mesh_facets.getConnectivity(current_facet),
			subfacet_connectivity)) {
			continue;
			}

			facet_list.push_back(current_facet);

			if (subfacet_list != nullptr) {
			auto && subfacets_of_facet =
			mesh_facets.getSubelementToElement(current_facet);

			/// check subfacets
			for (const auto & current_subfacet : subfacets_of_facet) {
			if (current_subfacet == ElementNull) {
			continue;
			}

			auto find_subfacet = std::find(
			subfacet_list->begin(), subfacet_list->end(), current_subfacet);
			if ((find_subfacet != subfacet_list->end()) or
			not hasElement(mesh_facets.getConnectivity(current_subfacet),
			subfacet_connectivity)) {
			continue;
			}

			subfacet_list->push_back(current_subfacet);
			}
			}

			/// consider opposing element
			const auto & elements_to_facet =
			mesh_facets.getElementToSubelement(current_facet);
			Idx opposing = 0;
			if (elements_to_facet[0] == current_element) {
			opposing = 1;
			}

			const auto & opposing_element = elements_to_facet[opposing];

			/// skip null elements since they are on a boundary
			if (opposing_element == ElementNull) {
			continue;
			}

			/// skip this element if already added
			if (std::find(element_list.begin(), element_list.end(),
			opposing_element) != element_list.end()) {
			continue;
			}

			/// only regular elements have to be checked
			if (opposing_element.kind() == _ek_regular) {
			elements_to_check.push(opposing_element);
			}

			element_list.push_back(opposing_element);

			AKANTU_DEBUG_ASSERT(hasElement(mesh.getConnectivity(opposing_element),
			subfacet_connectivity),
			"Subfacet doesn't belong to this element");
			}

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

			return facet_matched;
			}

			/* -------------------------------------------------------------------------- */
			void CohesiveElementInserterHelper::updateCohesiveData() {
			auto spatial_dimension = mesh.getSpatialDimension();
			auto third_dimension = spatial_dimension == 3;

			MeshAccessor mesh_accessor(mesh);
			MeshAccessor mesh_facets_accessor(mesh_facets);

			for (auto ghost_type : ghost_types) {
			for (auto cohesive_type : mesh.elementTypes(_element_kind = _ek_cohesive)) {
			auto nb_cohesive_elements = mesh.getNbElement(cohesive_type, ghost_type);
			nb_new_facets(cohesive_type, ghost_type) = nb_cohesive_elements;

			mesh_facets_accessor.getSubelementToElement(cohesive_type, ghost_type);
			}
			}

			auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];

			new_elements.reserve(new_elements.size() + facets_to_double.size());

			std::array<Element, 2> facets;

			auto & element_to_facet = mesh_facets_accessor.getElementToSubelement();
			auto & facets_to_cohesive_elements =
			mesh_facets_accessor.getSubelementToElement();

			for (auto && facet_to_double : make_view(facets_to_double, 2)) {
			auto & old_facet = facet_to_double[0];

			/// (in 3D cohesive elements connectivity is inverted)
			facets[third_dimension ? 1 : 0] = old_facet;
			facets[third_dimension ? 0 : 1] = facet_to_double[1];

			auto type_cohesive = FEEngine::getCohesiveElementType(old_facet.type);

			auto & facet_connectivities =
			mesh_facets.getConnectivity(old_facet.type, old_facet.ghost_type);
			auto facet_connectivity_it =
			facet_connectivities.begin(facet_connectivities.getNbComponent());

			auto cohesive_element = Element{
			type_cohesive, nb_new_facets(type_cohesive, old_facet.ghost_type)++,
			old_facet.ghost_type};

			auto & cohesives_connectivities =
			mesh_accessor.getConnectivity(type_cohesive, old_facet.ghost_type);
			Matrix<Idx> connectivity(facet_connectivities.getNbComponent(), 2);
			Vector<Element, 2> facets_to_cohesive_element;

			for (auto s : arange(2)) {
			/// store doubled facets
			facets_to_cohesive_element[s] = facets[s];

			// update connectivities
			connectivity(s) = facet_connectivity_it[facets[s].element];

			/// update element_to_facet vectors
			element_to_facet(facets[s], 0)[1] = cohesive_element;
			}

			cohesives_connectivities.push_back(connectivity);

			facets_to_cohesive_elements(type_cohesive, old_facet.ghost_type)
			.push_back(facets_to_cohesive_element);

			/// add cohesive element to the element event list
			new_elements.push_back(cohesive_element);
			}
			}

			/* -------------------------------------------------------------------------- */
			void CohesiveElementInserterHelper::doublePointFacet() {
			auto spatial_dimension = mesh.getSpatialDimension();
			if (spatial_dimension != 1) {
			return;
			}

			NewElementsEvent new_facets_event;

			auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
			const auto & element_to_facet = mesh_facets.getElementToSubelement();
			auto & position = mesh.getNodes();
			MeshAccessor mesh_facets_accessor(mesh_facets);
			MeshAccessor mesh_accessor(mesh);

			for (auto ghost_type : ghost_types) {
			for (auto facet_type : nb_new_facets.elementTypes(
			spatial_dimension - 1, ghost_type, _ek_regular)) {
			auto nb_new_element = nb_new_facets(facet_type, ghost_type);
			auto & connectivities =
			mesh_facets_accessor.getConnectivity(facet_type, ghost_type);
			connectivities.resize(nb_new_element);
			}
			}

			position.reserve(position.size() + facets_to_double.size());
			for (auto facet_to_double : make_view(facets_to_double, 2)) {
			auto & old_facet = facet_to_double[0];
			auto & new_facet = facet_to_double[1];

			auto element = element_to_facet(new_facet)[0];
			auto & facet_connectivities = mesh_facets_accessor.getConnectivity(
			old_facet.type, old_facet.ghost_type);
			auto old_node = facet_connectivities(old_facet.element);
			auto new_node = position.size();

			/// update position
			position.push_back(position(old_node));
			facet_connectivities(new_facet.element) = new_node;

			auto && segment_conectivity = mesh_accessor.getConnectivity(element);

			/// update facet connectivity
			auto it = std::find(segment_conectivity.begin(), segment_conectivity.end(),
			old_node);
			*it = new_node;

			doubled_nodes.push_back(Vector<Idx>{old_node, new_node});
			new_facets_event.getList().push_back(new_facet);
			}

			mesh_facets.sendEvent(new_facets_event);
			}

			/* -------------------------------------------------------------------------- */
			template <Int spatial_dimension>
			void CohesiveElementInserterHelper::doubleSubfacet() {
			if (spatial_dimension == 1) {
			return;
			}

			NewElementsEvent new_facets_event;

			std::vector<Idx> nodes_to_double;
			MeshAccessor mesh_accessor(mesh_facets);
			auto & connectivities = mesh_accessor.getConnectivities();

			auto & facets_to_double = *facets_to_double_by_dim[0];
			ElementTypeMap<Int> nb_new_elements;
			for (auto && facet_to_double : make_view(facets_to_double, 2)) {
			auto && old_element = facet_to_double[0];
			nb_new_elements(old_element.type, old_element.ghost_type) = 0;
			}

			for (auto && facet_to_double : make_view(facets_to_double, 2)) {
			auto && old_element = facet_to_double[0];
			++nb_new_elements(old_element.type, old_element.ghost_type);
			}

			for (auto ghost_type : ghost_types) {
			for (auto facet_type :
			nb_new_elements.elementTypes(0, ghost_type, _ek_regular)) {
			auto & connectivities =
			mesh_accessor.getConnectivity(facet_type, ghost_type);
			connectivities.resize(connectivities.size() +
			nb_new_elements(facet_type, ghost_type));
			}
			}

			for (auto && facet_to_double : make_view(facets_to_double, 2)) {
			auto & old_facet = facet_to_double(0);
			// auto & new_facet = facet_to_double(1);

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

			/// list nodes double
			nodes_to_double.push_back(connectivities(old_facet));
			}

			auto old_nb_doubled_nodes = doubled_nodes.size();
			doubleNodes(nodes_to_double);

			auto double_nodes_view = make_view(doubled_nodes, 2);

			for (auto && data :
			zip(make_view(facets_to_double, 2),
			range(double_nodes_view.begin() + old_nb_doubled_nodes,
			double_nodes_view.end()),
			arange(facets_to_double.size()))) {
			// auto & old_facet = std::get<0>(data)[0];
			auto & new_facet = std::get<0>(data)[1];

			new_facets_event.getList().push_back(new_facet);

			auto & nodes = std::get<1>(data);
			auto old_node = nodes(0);
			auto new_node = nodes(1);

			auto f = std::get<2>(data);

			/// add new nodes in connectivity
			connectivities(new_facet) = new_node;

			updateElementalConnectivity(mesh, old_node, new_node,
			elementsOfDimToElementsOfDim(2, 0)(f),
			&elementsOfDimToElementsOfDim(1, 0)(f));

			updateElementalConnectivity(mesh_facets, old_node, new_node,
			elementsOfDimToElementsOfDim(1, 0)(f));

			if (spatial_dimension == 3) {
			updateElementalConnectivity(mesh_facets, old_node, new_node,
			elementsOfDimToElementsOfDim(0, 0)(f));
			}
			}

			updateSubelementToElement(0, spatial_dimension == 2);
			updateElementToSubelement(0, spatial_dimension == 2);

			mesh_facets.sendEvent(new_facets_event);
			}

			} // namespace akantu

rAKA/src/mesh_utils/cohesive_element_inserter_helper.ccbf637b29f909bugfix/69-segmentation…

cohesive_element_inserter_helper.ccView Options

rAKA/src/mesh_utils/cohesive_element_inserter_helper.cc
bf637b29f909bugfix/69-segmentation…

cohesive_element_inserter_helper.cc
View Options