ntrf_contact.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, Jul 14, 16:35

ntrf_contact.cc
View Options

	/**
	* @file ntrf_contact.cc
	* @author David Kammer <david.kammer@epfl.ch>
	* @date Thu Mar 14 14:16:07 2013
	*
	* @brief
	*
	* @section LICENSE
	*
	* Copyright (©) 2010-2011 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 "ntrf_contact.hh"

	__BEGIN_SIMTOOLS__

	/* -------------------------------------------------------------------------- */
	NTRFContact::NTRFContact(SolidMechanicsModel & model,
	const ContactID & id,
	const MemoryID & memory_id) :
	Memory(memory_id), id(id), model(model),
	slaves(0,1,0,id+":slaves",std::numeric_limits<UInt>::quiet_NaN(),"slaves"),
	contact_pressure(0,model.getSpatialDimension(),0,id+":contact_pressure",
	std::numeric_limits<Real>::quiet_NaN(),"contact_pressure"),
	is_in_contact(0,1,false,id+":is_in_contact",false,"is_in_contact"),
	lumped_boundary(0,1,0,id+":lumped_boundary",
	std::numeric_limits<Real>::quiet_NaN(),"lumped_boundary"),
	reference_point(0,model.getSpatialDimension()),
	normal(0,model.getSpatialDimension()),
	contact_surfaces()
	{
	AKANTU_DEBUG_IN();

	FEM & boundary_fem = this->model.getFEMBoundary();
	boundary_fem.initShapeFunctions();

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::setReferencePoint(Real x, Real y, Real z) {
	AKANTU_DEBUG_IN();

	Real coord[3];
	coord[0] = x;
	coord[1] = y;
	coord[2] = z;

	this->reference_point.resize(1);

	UInt dim = this->model.getSpatialDimension();
	for (UInt d=0; d<dim; ++d)
	this->reference_point(0,d) = coord[d];

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::setNormal(Real x, Real y, Real z) {
	AKANTU_DEBUG_IN();

	UInt dim = this->model.getSpatialDimension();

	Real coord[3];
	coord[0] = x;
	coord[1] = y;
	coord[2] = z;

	if (dim == 2)
	Math::normalize2(coord);
	else if (dim == 3)
	Math::normalize3(coord);
	else
	AKANTU_DEBUG_ERROR("setNormal in NTRFContact not implemented for dimension " << dim);

	this->normal.resize(1);

	for (UInt d=0; d<dim; ++d)
	this->normal(0,d) = coord[d];


	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::addSurface(Surface surf) {
	AKANTU_DEBUG_IN();

	UInt dim = this->model.getSpatialDimension();
	this->contact_surfaces.insert(surf);

	// get all nodes for all surfaces
	CSR<UInt> all_surface_nodes;
	MeshUtils::buildNodesPerSurface(this->model.getFEM().getMesh(), all_surface_nodes);

	// find slave nodes
	for(CSR<UInt>::iterator snode = all_surface_nodes.begin(surf);
	snode != all_surface_nodes.end(surf);
	++snode) {
	this->addNode(*snode);
	}

	// synchronize with depending nodes
	updateInternalData();
	syncArrays(_added);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::addNode(UInt node) {
	AKANTU_DEBUG_IN();

	UInt dim = this->model.getSpatialDimension();

	// add to node arrays
	this->slaves.push_back(node);

	// set contact as false
	this->is_in_contact.push_back(false);

	// before initializing
	// set contact pressure, normal, lumped_boundary to Nan
	this->contact_pressure.push_back(std::numeric_limits<Real>::quiet_NaN());
	this->lumped_boundary.push_back(std::numeric_limits<Real>::quiet_NaN());

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::addNodes(Array<UInt> & nodes) {
	AKANTU_DEBUG_IN();

	UInt nb_nodes = nodes.getSize();
	UInt nb_compo = nodes.getNbComponent();
	for (UInt n=0; n<nb_nodes; ++n) {
	for (UInt c=0; c<nb_compo; ++c) {
	this->addNode(nodes(n,c));
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::registerSyncronizedArray(SyncronizedArrayBase & array) {
	AKANTU_DEBUG_IN();

	this->slaves.registerDependingArray(array);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::dumpRestart(const std::string & file_name) const {
	AKANTU_DEBUG_IN();

	this->slaves.dumpRestartFile(file_name);
	this->is_in_contact.dumpRestartFile(file_name);
	this->contact_pressure.dumpRestartFile(file_name);
	this->lumped_boundary.dumpRestartFile(file_name);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::readRestart(const std::string & file_name) {
	AKANTU_DEBUG_IN();

	this->slaves.readRestartFile(file_name);
	this->is_in_contact.readRestartFile(file_name);
	this->contact_pressure.readRestartFile(file_name);
	this->lumped_boundary.readRestartFile(file_name);

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::updateInternalData() {
	AKANTU_DEBUG_IN();

	updateLumpedBoundary();

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::updateLumpedBoundary() {
	AKANTU_DEBUG_IN();

	// set all values in lumped_boundary to zero
	this->lumped_boundary.clear();

	UInt dim = this->model.getSpatialDimension();
	UInt nb_contact_nodes = getNbContactNodes();

	const FEM & boundary_fem = this->model.getFEMBoundary();

	const Mesh & mesh = this->model.getFEM().getMesh();
	Mesh::type_iterator it = mesh.firstType(dim-1);
	Mesh::type_iterator last = mesh.lastType(dim-1);
	for (; it != last; ++it) {
	// get elements connected to each node
	CSR<UInt> node_to_element;
	MeshUtils::buildNode2ElementsByElementType(mesh, *it, node_to_element);

	UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it);
	const Array<UInt> & connectivity = mesh.getConnectivity(*it);

	// get shapes and compute integral
	const Array<Real> & shapes = boundary_fem.getShapes(*it);
	Array<Real> area(mesh.getNbElement(*it),nb_nodes_per_element);
	boundary_fem.integrate(shapes,area,nb_nodes_per_element,*it);

	// get surface id information
	const Array<UInt> & surface_id = mesh.getSurfaceID(*it);
	std::set<UInt>::iterator pos;
	std::set<UInt>::iterator end = this->contact_surfaces.end();

	if (this->contact_surfaces.size() == 0)
	std::cerr << "No surfaces in ntrf contact. You have to define the lumped boundary by yourself." << std::endl;

	// loop over contact nodes
	for (UInt i=0; i<nb_contact_nodes; ++i) {
	UInt node = this->slaves(i);

	CSR<UInt>::iterator elem = node_to_element.begin(node);
	// loop over all elements connected to this node
	for (; elem != node_to_element.end(node); ++elem) {
	UInt e = *elem;

	// if element is not at interface continue
	pos = this->contact_surfaces.find(surface_id(e));
	if (pos == end)
	continue;

	// loop over all points of this element
	for (UInt q=0; q<nb_nodes_per_element; ++q) {
	if (connectivity(e,q) == node) {
	this->lumped_boundary(i) += area(e,q);
	}
	}
	}
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::computeContactPressure() {
	AKANTU_DEBUG_IN();

	UInt dim = this->model.getSpatialDimension();
	Real delta_t = this->model.getTimeStep();
	UInt nb_contact_nodes = getNbContactNodes();

	// pre-compute the acceleration
	// (not increment acceleration, because residual is still Kf)
	Array<Real> acceleration(this->model.getFEM().getMesh().getNbNodes(),dim);
	this->model.solveLumped(acceleration,
	this->model.getMass(),
	this->model.getResidual(),
	this->model.getBoundary(),
	this->model.getF_M2A());

	const Array<Real> & residual = this->model.getResidual();
	const Array<Real> & mass = this->model.getMass();

	// compute relative normal fields of displacement, velocity and acceleration
	Array<Real> r_disp(0,1);
	Array<Real> r_velo(0,1);
	Array<Real> r_acce(0,1);
	Array<Real> r_old_acce(0,1);
	computeNormalGap(r_disp);
	computeRelativeNormalField(this->model.getVelocity(), r_velo);
	computeRelativeNormalField(acceleration, r_acce);
	computeRelativeNormalField(this->model.getAcceleration(), r_old_acce);

	AKANTU_DEBUG_ASSERT(r_disp.getSize() == nb_contact_nodes,
	"computeNormalGap does not give back arrays "
	<< "size == nb_contact_nodes. nb_contact_nodes = "
	<< nb_contact_nodes << " \| array size = "
	<< r_disp.getSize());
	AKANTU_DEBUG_ASSERT(r_velo.getSize() == nb_contact_nodes,
	"computeRelativeNormalField does not give back arrays "
	<< "size == nb_contact_nodes. nb_contact_nodes = "
	<< nb_contact_nodes << " \| array size = "
	<< r_velo.getSize());

	// compute gap array for all nodes
	Array<Real> gap(nb_contact_nodes, 1);
	Real * gap_p = gap.storage();
	Real * r_disp_p = r_disp.storage();
	Real * r_velo_p = r_velo.storage();
	Real * r_acce_p = r_acce.storage();
	Real * r_old_acce_p = r_old_acce.storage();
	for (UInt i=0; i<nb_contact_nodes; ++i) {
	gap_p = r_disp_p + delta_t * r_velo_p + delta_t delta_t * r_acce_p - 0.5 delta_t * delta_t * *r_old_acce_p;
	// increment pointers
	gap_p++;
	r_disp_p++;
	r_velo_p++;
	r_acce_p++;
	r_old_acce_p++;
	}

	// check if gap is negative -> is in contact
	for (UInt n=0; n<nb_contact_nodes; ++n) {
	if (gap(n) < 0.) {
	UInt node = this->slaves(n);
	for (UInt d=0; d<dim; ++d) {
	this->contact_pressure(n,d) = mass(node) / delta_t / delta_t
	/ this->lumped_boundary(n) * gap(n) * this->normal(0,d);
	/*
	// from the ntn_contact:
	this->contact_pressure(n,d) = this->impedance(n) * gap(n) / (2 * delta_t)
	* this->normal(0,d);
	*/
	}
	this->is_in_contact(n) = true;
	}
	else {
	for (UInt d=0; d<dim; ++d)
	this->contact_pressure(n,d) = 0.;
	this->is_in_contact(n) = false;
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::applyContactPressure() {
	AKANTU_DEBUG_IN();

	UInt nb_contact_nodes = getNbContactNodes();
	UInt dim = this->model.getSpatialDimension();

	Array<Real> & residual = this->model.getResidual();

	for (UInt n=0; n<nb_contact_nodes; ++n) {
	UInt node = this->slaves(n);
	for (UInt d=0; d<dim; ++d) {
	residual(node, d) -= this->lumped_boundary(n) * this->contact_pressure(n,d);
	}
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::computeRelativeTangentialField(const Array<Real> & field,
	Array<Real> & rel_tang_field) const {
	AKANTU_DEBUG_IN();

	// resize arrays to zero
	rel_tang_field.resize(0);

	UInt dim = this->model.getSpatialDimension();
	Real * field_p = field.storage();
	Real * normals_p = this->normal.storage();

	UInt nb_contact_nodes = this->slaves.getSize();
	for (UInt n=0; n<nb_contact_nodes; ++n) {
	// nodes
	UInt node = this->slaves(n);

	// compute relative field to master
	Real rel_array[dim];
	for (UInt d=0; d<dim; ++d) {
	rel_array[d] = field_p[node*dim + d];
	}

	// compute dot product with normal of master
	Real dot_prod = Math::vectorDot(normals_p, rel_array, dim);

	// compute the tangential projection of the relative field to the master
	for (UInt d=0; d<dim; ++d)
	rel_array[d] -= dot_prod * normals_p[d];

	rel_tang_field.push_back(rel_array);
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::computeNormalGap(Array<Real> & gap) const {
	AKANTU_DEBUG_IN();

	gap.resize(0);

	UInt dim = this->model.getSpatialDimension();
	Real * cur_pos = this->model.getCurrentPosition().storage();
	Real * normals_p = this->normal.storage();

	UInt nb_contact_nodes = this->getNbContactNodes();
	for (UInt n=0; n<nb_contact_nodes; ++n) {
	// nodes
	UInt node = this->slaves(n);

	// compute relative field to master
	Real rel_array[dim];
	for (UInt d=0; d<dim; ++d) {
	rel_array[d] = cur_pos[node*dim + d] - this->reference_point(0,d);
	}

	// compute dot product with normal of master
	Real dot_prod = Math::vectorDot(normals_p, rel_array, dim);

	gap.push_back(dot_prod);
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::computeRelativeNormalField(const Array<Real> & field,
	Array<Real> & rel_normal_field) const {
	AKANTU_DEBUG_IN();

	// resize arrays to zero
	rel_normal_field.resize(0);

	UInt dim = this->model.getSpatialDimension();
	Real * field_p = field.storage();
	Real * normals_p = this->normal.storage();

	UInt nb_contact_nodes = this->getNbContactNodes();
	for (UInt n=0; n<nb_contact_nodes; ++n) {
	// nodes
	UInt node = this->slaves(n);

	// compute dot product with normal of master
	Real dot_prod = Math::vectorDot(normals_p, &(field_p[node*dim]), dim);

	rel_normal_field.push_back(dot_prod);
	}

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	Int NTRFContact::getNodeIndex(UInt node) const {
	AKANTU_DEBUG_IN();

	AKANTU_DEBUG_OUT();
	return this->slaves.find(node);
	}

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

	stream << space << "NTRFContact [" << std::endl;
	stream << space << " + id : " << id << std::endl;
	stream << space << " + slaves : " << std::endl;
	this->slaves.printself(stream, indent + 2);
	stream << space << " + contact_pressure : " << std::endl;
	this->contact_pressure.printself(stream, indent + 2);

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

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void NTRFContact::syncArrays(SyncChoice sync_choice) {
	AKANTU_DEBUG_IN();

	this->slaves.syncElements(sync_choice);
	this->is_in_contact.syncElements(sync_choice);
	this->contact_pressure.syncElements(sync_choice);
	this->lumped_boundary.syncElements(sync_choice);

	AKANTU_DEBUG_OUT();
	}

	__END_SIMTOOLS__

ntrf_contact.ccNo OneTemporaryActions

File Metadata

ntrf_contact.ccView Options

Event Timeline

ntrf_contact.cc
No OneTemporary
Actions

ntrf_contact.cc
View Options