contact_detector.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Nov 8, 13:21

contact_detector.cc
View Options

	/**
	* @file contact_detector.cc
	*
	* @author Mohit Pundir <mohit.pundir@epfl.ch>
	*
	* @date creation: Wed Sep 12 2018
	* @date last modification: Fri Sep 21 2018
	*
	* @brief Mother class for all detection algorithms
	*
	* @section LICENSE
	*
	* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
	* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
	*
	* Akantu is free software: you can redistribute it and/or modify it under the
	* terms of the GNU Lesser General Public License as published by the Free
	* Software Foundation, either version 3 of the License, or (at your option) any
	* later version.
	*
	* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
	* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
	* details.
	*
	* You should have received a copy of the GNU Lesser General Public License
	* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
	*
	*/

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

	namespace akantu {

	/* -------------------------------------------------------------------------- */
	ContactDetector::ContactDetector(Mesh & mesh, const ID & id, UInt memory_id)
	: ContactDetector(mesh, mesh.getNodes(), id, memory_id) {

	}

	/* -------------------------------------------------------------------------- */
	ContactDetector::ContactDetector(Mesh & mesh, Array<Real> positions, const ID & id, UInt memory_id)
	: Memory(id, memory_id),
	Parsable(ParserType::_contact_detector, id),
	mesh(mesh) {

	AKANTU_DEBUG_IN();

	this->spatial_dimension = mesh.getSpatialDimension();

	this->positions = positions;

	this->mesh.fillNodesToElements(this->spatial_dimension - 1);

	this->parseSection();

	this->computeMaximalDetectionDistance();

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void ContactDetector::parseSection() {

	const Parser & parser = getStaticParser();

	const ParserSection & section =
	*(parser.getSubSections(ParserType::_contact_detector).first);

	auto type = section.getParameterValue<std::string>("type");

	if (type == "implicit") {
	this->detection_type = _implicit;
	}
	else if (type == "explicit"){
	this->detection_type = _explicit;
	}
	else {
	AKANTU_ERROR("Unknown detection type : " << type);
	}

	surfaces[Surface::master] = section.getParameterValue<std::string>("master");
	surfaces[Surface::slave ] = section.getParameterValue<std::string>("slave");

	two_pass_algorithm = section.getParameterValue<bool>("two_pass_algorithm");
	}

	/* -------------------------------------------------------------------------- */
	void ContactDetector::search(std::map<UInt, ContactElement> & contact_map) {

	contact_pairs.clear();

	SpatialGrid<UInt> master_grid(spatial_dimension);
	SpatialGrid<UInt> slave_grid(spatial_dimension);

	this->globalSearch(slave_grid, master_grid);

	this->localSearch(slave_grid, master_grid);

	//if (two_pass_algorithm) {
	// this->localSearch(master_grid, slave_grid);
	//}

	this->constructContactMap(contact_map);
	}

	/* -------------------------------------------------------------------------- */
	void ContactDetector::globalSearch(SpatialGrid<UInt> & slave_grid,
	SpatialGrid<UInt> & master_grid) {

	auto & master_list =
	mesh.getElementGroup(surfaces[Surface::master]).getNodeGroup().getNodes();

	auto & slave_list =
	mesh.getElementGroup(surfaces[Surface::slave]).getNodeGroup().getNodes();

	BBox bbox_master(spatial_dimension);
	this->constructBoundingBox(bbox_master, master_list);

	BBox bbox_slave(spatial_dimension);
	this->constructBoundingBox(bbox_slave, slave_list);

	auto && bbox_intersection =
	bbox_master.intersection(bbox_slave);

	AKANTU_DEBUG_INFO( "Intersection BBox "
	<< bbox_intersection );

	Vector<Real> center(spatial_dimension);
	bbox_intersection.getCenter(center);

	Vector<Real> spacing(spatial_dimension);
	this->computeCellSpacing(spacing);

	master_grid.setCenter(center);
	master_grid.setSpacing(spacing);
	this->constructGrid(master_grid, bbox_intersection, master_list);

	slave_grid.setCenter(center);
	slave_grid.setSpacing(spacing);
	this->constructGrid(slave_grid, bbox_intersection, slave_list);

	// search slave grid nodes in contactelement array and if they exits
	// and still have orthogonal projection on its associated master
	// facetremove it from the spatial grid or do not consider it for
	// local search, maybe better option will be to have spatial grid of
	// type node info and one of the variable of node info should be
	// facet already exits
	// so contact elements will be updated based on the above
	// consideration , this means only those contact elements will be
	// keep whose slave node is still in intersection bbox and still has
	// projection in its master facet
	// also if slave node is already exists in contact element and
	// orthogonal projection does not exits then search the associated
	// master facets with the current master facets within a given
	// radius , this is subjected to computational cost as searching
	// neighbbor cells can be more effective.
	}

	/* -------------------------------------------------------------------------- */
	void ContactDetector::localSearch(SpatialGrid<UInt> & slave_grid,
	SpatialGrid<UInt> & master_grid) {

	// local search
	// out of these array check each cell for closet node in that cell
	// and neighbouring cells find the actual orthogonally closet
	// check the projection of slave node on master facets connected to
	// the closet master node, if yes update the contact element with
	// slave node and master node and master surfaces connected to the
	// master node
	// these master surfaces will be needed later to update contact
	// elements

	/// find the closet master node for each slave node
	for (auto && cell_id : slave_grid) {
	/// loop over all the slave nodes of the current cell
	for (auto && slave_node: slave_grid.getCell(cell_id)) {

	bool pair_exists = false;

	Vector<Real> pos(spatial_dimension);
	for (UInt s: arange(spatial_dimension))
	pos(s) = this->positions(slave_node, s);

	Real closet_distance = std::numeric_limits<Real>::max();
	UInt closet_master_node;

	/// loop over all the neighboring cells of the current cell
	for (auto && neighbor_cell : cell_id.neighbors()) {
	/// loop over the data of neighboring cells from master grid
	for (auto && master_node : master_grid.getCell(neighbor_cell)) {

	Vector<Real> pos2(spatial_dimension);
	for (UInt s: arange(spatial_dimension))
	pos2(s) = this->positions(master_node, s);

	Real distance = pos.distance(pos2);

	if (distance <= closet_distance) {
	closet_master_node = master_node;
	closet_distance = distance;
	pair_exists = true;
	}
	}
	}

	if (pair_exists)
	contact_pairs.push_back( std::make_pair(slave_node,
	closet_master_node));
	}
	}
	}

	/* -------------------------------------------------------------------------- */
	void ContactDetector::constructContactMap(std::map<UInt, ContactElement> & contact_map) {

	auto surface_dimension = spatial_dimension - 1;

	std::map<UInt, ContactElement> previous_contact_map = contact_map;
	contact_map.clear();

	auto get_index = [&](auto & gaps, auto & projections) {

	UInt index;
	Real gap_min = std::numeric_limits<Real>::max();

	UInt counter = 0;
	for (auto && values : zip(gaps,
	make_view(projections, surface_dimension))) {
	auto & gap = std::get<0>(values);
	auto & projection = std::get<1>(values);

	/// todo adhoc fix to assign a master element in case the
	/// projection does not lie in the extended element. As it is
	/// tolerance based
	bool is_valid = this->checkValidityOfProjection(projection);

	if (is_valid and gap <= gap_min) {
	gap_min = gap;
	index = counter;
	}
	counter++;
	}

	if (index >= gaps.size()) {
	auto gap_min_it = std::min_element(gaps.begin(), gaps.end());
	auto index_it = std::find(gaps.begin(), gaps.end(), *gap_min_it);
	index = *index_it;
	}

	return index;
	};


	auto get_connectivity = [&](auto & slave, auto & master) {
	Vector<UInt> master_conn = this->mesh.getConnectivity(master);
	Vector<UInt> elem_conn(master_conn.size() + 1);

	elem_conn[0] = slave;
	for (UInt i = 1; i < elem_conn.size(); ++i) {
	elem_conn[i] = master_conn[i-1];
	}

	return elem_conn;
	};

	for (auto & pairs : contact_pairs) {

	const auto & slave_node = pairs.first;
	const auto & master_node = pairs.second;

	Array<Element> elements;
	this->mesh.getAssociatedElements(master_node, elements);

	Array<Real> gaps(elements.size(), 1, "gaps");
	Array<Real> normals(elements.size(), spatial_dimension, "normals");
	Array<Real> projections(elements.size(), surface_dimension, "projections");

	this->computeOrthogonalProjection(slave_node, elements,
	normals, gaps, projections);

	auto index = get_index(gaps, projections);

	auto connectivity = get_connectivity(slave_node, elements[index]);

	contact_map[slave_node].setMaster(elements[index]);
	contact_map[slave_node].setGap(gaps[index]);
	contact_map[slave_node].setNormal(Vector<Real>(normals.begin(spatial_dimension)[index], true));
	contact_map[slave_node].setProjection(Vector<Real>(projections.begin(surface_dimension)[index], true));
	contact_map[slave_node].setConnectivity(connectivity);

	Matrix<Real> tangents(surface_dimension, spatial_dimension);
	this->computeTangentsOnElement(contact_map[slave_node].master,
	contact_map[slave_node].projection,
	tangents);

	switch (spatial_dimension) {
	case 2: {
	Vector<Real> e_z(3);
	e_z[0] = 0.;
	e_z[1] = 0.;
	e_z[2] = 1.;

	Vector<Real> tangent(3);
	tangent[0] = tangents(0, 0);
	tangent[1] = tangents(0, 1);
	tangent[2] = 0.;

	auto exp_normal = e_z.crossProduct(tangent);

	auto & cal_normal = contact_map[slave_node].normal;

	auto ddot = cal_normal.dot(exp_normal);
	if (ddot < 0) {
	tangents *= -1.0;
	}

	break;
	}
	case 3: {
	auto tang_trans = tangents.transpose();
	auto tang1 = Vector<Real>(tang_trans(0));
	auto tang2 = Vector<Real>(tang_trans(1));

	auto tang1_cross_tang2 = tang1.crossProduct(tang2);
	auto exp_normal = tang1_cross_tang2 / tang1_cross_tang2.norm();

	auto & cal_normal = contact_map[slave_node].normal;

	auto ddot = cal_normal.dot(exp_normal);
	if (ddot < 0) {
	tang_trans(1) *= -1.0;
	}

	tangents = tang_trans.transpose();
	break;
	}
	default:
	break;
	}

	contact_map[slave_node].setTangent(tangents);

	auto search = previous_contact_map.find(slave_node);
	if (search != previous_contact_map.end()) {
	auto previous_projection = previous_contact_map[slave_node].getPreviousProjection();
	contact_map[slave_node].setPreviousProjection(previous_projection);

	auto previous_traction = previous_contact_map[slave_node].getTraction();
	contact_map[slave_node].setTraction(previous_traction);
	}
	else {
	Vector<Real> previous_projection(surface_dimension, 0.);
	contact_map[slave_node].setPreviousProjection(previous_projection);

	Vector<Real> previous_traction(surface_dimension, 0.);
	contact_map[slave_node].setTraction(previous_traction);
	}
	}

	contact_pairs.clear();
	}

	/* -------------------------------------------------------------------------- */
	void ContactDetector::computeOrthogonalProjection(const UInt & node,
	const Array<Element> & elements,
	Array<Real> & normals, Array<Real> & gaps,
	Array<Real> & projections) {

	Vector<Real> query(spatial_dimension);
	for (UInt s: arange(spatial_dimension))
	query(s) = this->positions(node, s);

	for (auto && values :
	zip( elements,
	gaps,
	make_view(normals , spatial_dimension),
	make_view(projections, spatial_dimension - 1))) {

	const auto & element = std::get<0>(values);
	auto & gap = std::get<1>(values);
	auto & normal = std::get<2>(values);
	auto & projection = std::get<3>(values);

	this->computeNormalOnElement(element, normal);

	Vector<Real> real_projection(spatial_dimension);
	this->computeProjectionOnElement(element, normal, query,
	projection, real_projection);

	gap = this->computeGap(query, real_projection, normal);
	}

	}


	/* -------------------------------------------------------------------------- */
	void ContactDetector::computeProjectionOnElement(const Element & element,
	const Vector<Real> & normal,
	const Vector<Real> & query,
	Vector<Real> & natural_projection,
	Vector<Real> & real_projection) {

	UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type);

	Matrix<Real> coords(spatial_dimension, nb_nodes_per_element);
	this->coordinatesOfElement(element, coords);

	Vector<Real> point(coords(0));
	Real alpha = (query - point).dot(normal);
	real_projection = query - alpha * normal;

	this->computeNaturalProjection(element, real_projection, natural_projection);
	}

	/* -------------------------------------------------------------------------- */
	void ContactDetector::computeNaturalProjection(const Element & element, Vector<Real> & real_projection,
	Vector<Real> & natural_projection) {

	const ElementType & type = element.type;
	UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
	UInt * elem_val = mesh.getConnectivity(type,
	_not_ghost).storage();
	Matrix<Real> nodes_coord(spatial_dimension, nb_nodes_per_element);

	mesh.extractNodalValuesFromElement(this->positions, nodes_coord.storage(),
	elem_val + element.element * nb_nodes_per_element,
	nb_nodes_per_element, spatial_dimension);

	#define GET_NATURAL_COORDINATE(type) \
	ElementClass<type>::inverseMap(real_projection, nodes_coord, natural_projection)
	AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_COORDINATE);
	#undef GET_NATURAL_COORDINATE
	}



	/* -------------------------------------------------------------------------- */
	void ContactDetector::computeTangentsOnElement(const Element & el, Vector<Real> & projection, Matrix<Real> & tangents) {

	const ElementType & type = el.type;

	UInt nb_nodes_master = Mesh::getNbNodesPerElement(type);

	Vector<Real> shapes(nb_nodes_master);
	Matrix<Real> shapes_derivatives(spatial_dimension - 1, nb_nodes_master);

	#define GET_SHAPES_NATURAL(type) \
	ElementClass<type>::computeShapes(projection, shapes)
	AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL);
	#undef GET_SHAPES_NATURAL

	#define GET_SHAPE_DERIVATIVES_NATURAL(type) \
	ElementClass<type>::computeDNDS(projection, shapes_derivatives)
	AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL);
	#undef GET_SHAPE_DERIVATIVES_NATURAL


	Matrix<Real> coords(spatial_dimension, nb_nodes_master);
	coordinatesOfElement(el, coords);

	tangents.mul<false, true>(shapes_derivatives, coords);

	auto temp_tangents = tangents.transpose();
	for (UInt i = 0; i < spatial_dimension -1; ++i) {
	auto temp = Vector<Real>(temp_tangents(i));
	temp_tangents(i) = temp.normalize();
	}

	tangents = temp_tangents.transpose();
	}


	} // akantu

contact_detector.ccNo OneTemporaryActions

File Metadata

contact_detector.ccView Options

Event Timeline

contact_detector.cc
No OneTemporary
Actions

contact_detector.cc
View Options