diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc
index 18748db70..a41cbd05d 100644
--- a/src/model/contact_mechanics/contact_detector.cc
+++ b/src/model/contact_mechanics/contact_detector.cc
@@ -1,291 +1,294 @@
/**
* @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,
const MemoryID & memory_id)
: ContactDetector(mesh, mesh.getNodes(), id, memory_id) {}
/* -------------------------------------------------------------------------- */
ContactDetector::ContactDetector(Mesh & mesh, Array<Real> positions,
const ID & id, const MemoryID & memory_id)
: Memory(id, memory_id), Parsable(ParserType::_contact_detector, id),
mesh(mesh), positions(0, mesh.getSpatialDimension()) {
AKANTU_DEBUG_IN();
this->spatial_dimension = mesh.getSpatialDimension();
this->positions.copy(positions);
this->parseSection();
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);
}
this->projection_tolerance =
section.getParameterValue<Real>("projection_tolerance");
+ this->max_iterations =
+ section.getParameterValue<Real>("max_iterations");
this->extension_tolerance =
section.getParameterValue<Real>("extension_tolerance");
}
/* -------------------------------------------------------------------------- */
void ContactDetector::search(Array<ContactElement> & elements,
Array<Real> & gaps, Array<Real> & normals,
Array<Real> & projections) {
UInt surface_dimension = spatial_dimension - 1;
this->mesh.fillNodesToElements(surface_dimension);
this->computeMaximalDetectionDistance();
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);
this->createContactElements(elements, gaps, normals, projections);
}
/* -------------------------------------------------------------------------- */
void ContactDetector::globalSearch(SpatialGrid<UInt> & slave_grid,
SpatialGrid<UInt> & master_grid) {
auto & master_list = surface_selector->getMasterList();
auto & slave_list = surface_selector->getSlaveList();
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)) {
/// check for self contact
if (slave_node == master_node)
continue;
bool is_valid = true;
Array<Element> elements;
this->mesh.getAssociatedElements(slave_node, elements);
for (auto & elem : elements) {
if (elem.kind() != _ek_regular)
continue;
Vector<UInt> connectivity =
const_cast<const Mesh &>(this->mesh).getConnectivity(elem);
auto node_iter = std::find(connectivity.begin(), connectivity.end(),
master_node);
if (node_iter != connectivity.end()) {
is_valid = false;
break;
}
}
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 and is_valid) {
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::createContactElements(Array<ContactElement> & contact_elements,
Array<Real> & gaps, Array<Real> & normals,
Array<Real> & projections) {
auto surface_dimension = spatial_dimension - 1;
Real alpha;
switch (detection_type) {
case _explicit: {
alpha = 1.0;
break;
}
case _implicit: {
alpha = -1.0;
break;
}
default:
AKANTU_EXCEPTION(detection_type << " is not a valid contact detection type");
break;
}
for (auto & pairs : contact_pairs) {
const auto & slave_node = pairs.first;
Vector<Real> slave(spatial_dimension);
for (UInt s : arange(spatial_dimension))
slave(s) = this->positions(slave_node, s);
const auto & master_node = pairs.second;
Array<Element> elements;
this->mesh.getAssociatedElements(master_node, elements);
auto & gap = gaps.begin()[slave_node];
Vector<Real> normal(normals.begin(spatial_dimension)[slave_node]);
Vector<Real> projection(projections.begin(surface_dimension)[slave_node]);
auto index = GeometryUtils::orthogonalProjection(mesh, positions, slave, elements,
gap, projection, normal, alpha,
+ this->max_iterations,
this->projection_tolerance,
this->extension_tolerance);
// if a valid projection is not found on the patch of elements
// index is -1 or if not a valid self contact, the contact element
// is not created
if (index == UInt(-1) or !isValidSelfContact(slave_node, gap, normal)) {
gap *= 0;
normal *= 0;
continue;
}
// create contact element
contact_elements.push_back(ContactElement(slave_node, elements[index]));
}
contact_pairs.clear();
}
} // namespace akantu
diff --git a/src/model/contact_mechanics/contact_detector.hh b/src/model/contact_mechanics/contact_detector.hh
index 099c8169e..5076fa00d 100644
--- a/src/model/contact_mechanics/contact_detector.hh
+++ b/src/model/contact_mechanics/contact_detector.hh
@@ -1,213 +1,216 @@
/**
* @file contact_detection.hh
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Wed Sep 12 2018
* @date last modification: Tue Oct 23 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 "aka_bbox.hh"
#include "aka_common.hh"
#include "aka_grid_dynamic.hh"
#include "aka_memory.hh"
#include "contact_element.hh"
#include "element_class.hh"
#include "element_group.hh"
#include "fe_engine.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "parsable.hh"
#include "surface_selector.hh"
#include "geometry_utils.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_CONTACT_DETECTOR_HH__
#define __AKANTU_CONTACT_DETECTOR_HH__
namespace akantu {
enum class Surface { master, slave };
/* -------------------------------------------------------------------------- */
class ContactDetector : private Memory,
public Parsable {
/* ------------------------------------------------------------------------ */
/* Constructor/Destructors */
/* ------------------------------------------------------------------------ */
public:
ContactDetector(Mesh &, const ID & id = "contact_detector",
const MemoryID & memory_id = 0);
ContactDetector(Mesh &, Array<Real> positions,
const ID & id = "contact_detector",
const MemoryID & memory_id = 0);
~ContactDetector() = default;
/* ------------------------------------------------------------------------ */
/* Members */
/* ------------------------------------------------------------------------ */
public:
/// performs all search steps
void search(Array<ContactElement> & contact_map,
Array<Real> & gaps, Array<Real> & normals,
Array<Real> & projections);
/// performs global spatial search to construct spatial grids
void globalSearch(SpatialGrid<UInt> &, SpatialGrid<UInt> &);
/// performs local search to find closet master node to a slave node
void localSearch(SpatialGrid<UInt> &, SpatialGrid<UInt> &);
/// create contact elements
void createContactElements(Array<ContactElement> & elements,
Array<Real> & gaps, Array<Real> & normals,
Array<Real> & projections);
private:
/// reads the input file to get contact detection options
void parseSection();
/* ------------------------------------------------------------------------ */
/* Inline Methods */
/* ------------------------------------------------------------------------ */
public:
/// checks whether the natural projection is valid or not
inline bool checkValidityOfProjection(Vector<Real> &);
/// extracts the coordinates of an element
inline void coordinatesOfElement(const Element &, Matrix<Real> &);
/// computes the optimal cell size for grid
inline void computeCellSpacing(Vector<Real> &);
/// constructs a grid containing nodes lying within bounding box
inline void constructGrid(SpatialGrid<UInt> &, BBox &, const Array<UInt> &);
/// constructs the bounding box based on nodes list
inline void constructBoundingBox(BBox &, const Array<UInt> &);
/// computes the maximum in radius for a given mesh
inline void computeMaximalDetectionDistance();
/// constructs the connectivity for a contact element
inline Vector<UInt> constructConnectivity(UInt &, const Element &);
/// computes normal on an element
inline void computeNormalOnElement(const Element &, Vector<Real> &);
/// extracts vectors which forms the plane of element
inline void vectorsAlongElement(const Element &, Matrix<Real> &);
/// computes the gap between slave and its projection on master
/// surface
inline Real computeGap(Vector<Real> &, Vector<Real> &);
/// filter boundary elements
inline void filterBoundaryElements(Array<Element> & elements,
Array<Element> & boundary_elements);
/// checks whether self contact condition leads to a master element
/// which is closet but not orthogonally opposite to slave surface
//inline bool checkValidityOfSelfContact(const UInt &, ContactElement &);
/// checks whether self contact condition leads to a master element
/// which is closet but not orthogonally opposite to slave surface
inline bool isValidSelfContact(const UInt &, const Real & , const Vector<Real> &);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get the mesh
AKANTU_GET_MACRO(Mesh, mesh, Mesh &)
/// returns the maximum detection distance
AKANTU_GET_MACRO(MaximumDetectionDistance, max_dd, Real);
AKANTU_SET_MACRO(MaximumDetectionDistance, max_dd, Real);
/// returns the bounding box extension
AKANTU_GET_MACRO(MaximumBoundingBox, max_bb, Real);
AKANTU_SET_MACRO(MaximumBoundingBox, max_bb, Real);
/// returns the minimum detection distance
AKANTU_GET_MACRO(MinimumDetectionDistance, min_dd, Real);
AKANTU_SET_MACRO(MinimumDetectionDistance, min_dd, Real);
AKANTU_GET_MACRO_NOT_CONST(Positions, positions, Array<Real> &);
AKANTU_SET_MACRO(Positions, positions, Array<Real>);
AKANTU_GET_MACRO_NOT_CONST(SurfaceSelector, *surface_selector,
SurfaceSelector &);
AKANTU_SET_MACRO(SurfaceSelector, surface_selector,
std::shared_ptr<SurfaceSelector>);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// maximal detection distance for grid spacing
Real max_dd;
/// minimal detection distance for grid spacing
Real min_dd;
/// maximal bounding box extension
Real max_bb;
/// tolerance for finding natural projection
Real projection_tolerance;
+ /// iterations for finding natural projection
+ UInt max_iterations;
+
/// tolerance for extending a master elements on all sides
Real extension_tolerance;
/// Mesh
Mesh & mesh;
/// dimension of the model
UInt spatial_dimension{0};
/// node selector for selecting master and slave nodes
std::shared_ptr<SurfaceSelector> surface_selector;
/// contact pair slave node to closet master node
std::vector<std::pair<UInt, UInt>> contact_pairs;
/// contains the updated positions of the nodes
Array<Real> positions;
/// type of detection explicit/implicit
DetectionType detection_type;
};
} // namespace akantu
#include "contact_detector_inline_impl.cc"
#endif /* __AKANTU_CONTACT_DETECTOR_HH__ */
diff --git a/src/model/contact_mechanics/geometry_utils.cc b/src/model/contact_mechanics/geometry_utils.cc
index 1153d8e4e..87b9ec1a5 100644
--- a/src/model/contact_mechanics/geometry_utils.cc
+++ b/src/model/contact_mechanics/geometry_utils.cc
@@ -1,438 +1,440 @@
/**
* @file geometry_utils.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Mon Mmay 20 2019
* @date last modification: Mon May 20 2019
*
* @brief Implementation of various utilities needed for contact geometry
*
* @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 "geometry_utils.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
void GeometryUtils::normal(const Mesh & mesh, const Array<Real> & positions,
const Element & element, Vector<Real> & normal, bool outward) {
UInt spatial_dimension = mesh.getSpatialDimension();
UInt surface_dimension = spatial_dimension - 1;
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type);
Matrix<Real> coords(spatial_dimension, nb_nodes_per_element);
UInt * elem_val = mesh.getConnectivity(element.type, _not_ghost).storage();
mesh.extractNodalValuesFromElement(positions, coords.storage(),
elem_val + element.element * nb_nodes_per_element,
nb_nodes_per_element, spatial_dimension);
Matrix<Real> vectors(spatial_dimension, surface_dimension);
switch (spatial_dimension) {
case 1: {
normal[0] = 1;
break;
}
case 2: {
vectors(0) = Vector<Real>(coords(1)) - Vector<Real>(coords(0));
Math::normal2(vectors.storage(), normal.storage());
break;
}
case 3: {
vectors(0) = Vector<Real>(coords(1)) - Vector<Real>(coords(0));
vectors(1) = Vector<Real>(coords(2)) - Vector<Real>(coords(0));
Math::normal3(vectors(0).storage(), vectors(1).storage(), normal.storage());
break;
}
default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); }
}
// to ensure that normal is always outwards from master surface
if (outward) {
const auto & element_to_subelement =
mesh.getElementToSubelement(element.type)(element.element);
Vector<Real> outside(spatial_dimension);
mesh.getBarycenter(element, outside);
// check if mesh facets exists for cohesive elements contact
Vector<Real> inside(spatial_dimension);
if (mesh.isMeshFacets()) {
mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside);
} else {
mesh.getBarycenter(element_to_subelement[0], inside);
}
Vector<Real> inside_to_outside = outside - inside;
auto projection = inside_to_outside.dot(normal);
if (projection < 0) {
normal *= -1.0;
}
}
}
/* -------------------------------------------------------------------------- */
void GeometryUtils::covariantBasis(const Mesh & mesh, const Array<Real> & positions,
const Element & element, const Vector<Real> & normal,
Vector<Real> & natural_coord,
Matrix<Real> & tangents) {
UInt spatial_dimension = mesh.getSpatialDimension();
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(positions, nodes_coord.storage(),
elem_val + element.element * nb_nodes_per_element,
nb_nodes_per_element, spatial_dimension);
UInt surface_dimension = spatial_dimension - 1;
Matrix<Real> dnds(surface_dimension, nb_nodes_per_element);
#define GET_SHAPE_DERIVATIVES_NATURAL(type) \
ElementClass<type>::computeDNDS(natural_coord, dnds)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL);
#undef GET_SHAPE_DERIVATIVES_NATURAL
tangents.mul<false, true>(dnds, nodes_coord);
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();
// to ensure that direction of tangents are correct, cross product
// of tangents should give the normal vector computed earlier
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 ddot = 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 ddot = normal.dot(exp_normal);
if (ddot < 0) {
tang_trans(1) *= -1.0;
}
tangents = tang_trans.transpose();
break;
}
default:
break;
}
}
/* -------------------------------------------------------------------------- */
void GeometryUtils::curvature(const Mesh & mesh, const Array<Real> & positions,
const Element & element, const Vector<Real> & natural_coord,
Matrix<Real> & curvature) {
UInt spatial_dimension = mesh.getSpatialDimension();
auto surface_dimension = spatial_dimension - 1;
const ElementType & type = element.type;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage();
Matrix<Real> dn2ds2(surface_dimension * surface_dimension,
Mesh::getNbNodesPerElement(type));
#define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \
ElementClass<type>::computeDN2DS2(natural_coord, dn2ds2)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL);
#undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL
Matrix<Real> coords(spatial_dimension, nb_nodes_per_element);
mesh.extractNodalValuesFromElement(positions, coords.storage(),
elem_val + element.element * nb_nodes_per_element,
nb_nodes_per_element, spatial_dimension);
curvature.mul<false, true>(coords, dn2ds2);
}
/* -------------------------------------------------------------------------- */
UInt GeometryUtils::orthogonalProjection(const Mesh & mesh, const Array<Real> & positions,
const Vector<Real> & slave,
const Array<Element> & elements,
Real & gap, Vector<Real> & natural_projection,
Vector<Real> & normal, Real alpha,
+ UInt max_iterations,
Real projection_tolerance,
Real extension_tolerance) {
UInt index = UInt(-1);
Real min_gap = std::numeric_limits<Real>::max();
UInt spatial_dimension = mesh.getSpatialDimension();
UInt nb_same_sides{0};
UInt nb_boundary_elements{0};
UInt counter{0};
for (auto & element : elements) {
if (!GeometryUtils::isBoundaryElement(mesh, element))
continue;
nb_boundary_elements++;
Vector<Real> normal_ele(spatial_dimension);
GeometryUtils::normal(mesh, positions, element, normal_ele);
Vector<Real> master(spatial_dimension);
GeometryUtils::realProjection(mesh, positions, slave, element, normal_ele, master);
Vector<Real> xi(natural_projection.size());
GeometryUtils::naturalProjection(mesh, positions, element, master, xi,
- projection_tolerance);
+ max_iterations, projection_tolerance);
// if gap between master projection and slave point is zero, then
// it means that slave point lies on the master element, hence the
// normal from master to slave cannot be computed in that case
auto master_to_slave = slave - master;
Real temp_gap = master_to_slave.norm();
if (temp_gap != 0)
master_to_slave /= temp_gap;
// also the slave point should lie inside the master surface in
// case of explicit or outside in case of implicit, one way to
// check that is by checking the dot product of normal at each
// master element, if the values of all dot product is same then
// the slave point lies on the same side of each master element
// A alpha parameter is introduced which is 1 in case of explicit
// and -1 in case of implicit, therefor the variation (dot product
// + alpha) should be close to zero (within tolerance) for both
// cases
Real direction_tolerance = 1e-8;
auto product = master_to_slave.dot(normal_ele);
auto variation = std::abs(product + alpha);
if (variation <= direction_tolerance and
temp_gap <= min_gap and
GeometryUtils::isValidProjection(xi, extension_tolerance)) {
gap = -temp_gap;
min_gap = temp_gap;
index = counter;
natural_projection = xi;
normal = normal_ele;
}
if(temp_gap == 0 or variation <= direction_tolerance)
nb_same_sides++;
counter++;
}
// if point is not on the same side of all the boundary elements
// than it is not consider even if the closet master element is
// found
if(nb_same_sides != nb_boundary_elements)
index = UInt(-1);
return index;
}
/* -------------------------------------------------------------------------- */
void GeometryUtils::realProjection(const Mesh & mesh, const Array<Real> & positions,
const Vector<Real> & slave, const Element & element,
const Vector<Real> & normal,
Vector<Real> & projection) {
UInt spatial_dimension = mesh.getSpatialDimension();
const ElementType & type = element.type;
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type);
UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage();
Matrix<Real> nodes_coord(spatial_dimension, nb_nodes_per_element);
mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(),
elem_val + element.element * nb_nodes_per_element,
nb_nodes_per_element, spatial_dimension);
Vector<Real> point(nodes_coord(0));
Real alpha = (slave - point).dot(normal);
projection = slave - alpha * normal;
}
/* -------------------------------------------------------------------------- */
void GeometryUtils::realProjection(const Mesh & mesh, const Array<Real> & positions,
const Element & element, const Vector<Real> & natural_coord,
Vector<Real> & projection) {
UInt spatial_dimension = mesh.getSpatialDimension();
const ElementType & type = element.type;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element.type);
Vector<Real> shapes(nb_nodes_per_element);
#define GET_SHAPE_NATURAL(type) \
ElementClass<type>::computeShapes(natural_coord, shapes)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL);
#undef GET_SHAPE_NATURAL
Matrix<Real> nodes_coord(spatial_dimension, nb_nodes_per_element);
UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage();
mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(),
elem_val + element.element * nb_nodes_per_element,
nb_nodes_per_element, spatial_dimension);
projection.mul<false>(nodes_coord, shapes);
}
/* -------------------------------------------------------------------------- */
void GeometryUtils::naturalProjection(const Mesh & mesh, const Array<Real> & positions,
const Element & element,
Vector<Real> & real_projection,
Vector<Real> & natural_projection,
+ UInt max_iterations,
Real projection_tolerance) {
UInt spatial_dimension = mesh.getSpatialDimension();
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(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, projection_tolerance)
+ natural_projection, max_iterations, projection_tolerance)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_COORDINATE);
#undef GET_NATURAL_COORDINATE
}
/* -------------------------------------------------------------------------- */
void GeometryUtils::contravariantBasis(const Matrix<Real> & covariant,
Matrix<Real> & contravariant) {
auto inv_A = GeometryUtils::contravariantMetricTensor(covariant);
contravariant.mul<false, false>(inv_A, covariant);
}
/* -------------------------------------------------------------------------- */
Matrix<Real> GeometryUtils::covariantMetricTensor(const Matrix<Real> & covariant_bases) {
Matrix<Real> A(covariant_bases.rows(), covariant_bases.rows());
A.mul<false, true>(covariant_bases, covariant_bases);
return A;
}
/* -------------------------------------------------------------------------- */
Matrix<Real> GeometryUtils::contravariantMetricTensor(const Matrix<Real> & covariant_bases) {
auto A = GeometryUtils::covariantMetricTensor(covariant_bases);
auto inv_A = A.inverse();
return inv_A;
}
/* -------------------------------------------------------------------------- */
Matrix<Real> GeometryUtils::covariantCurvatureTensor(const Mesh & mesh,
const Array<Real> & positions,
const Element & element,
const Vector<Real> & natural_coord,
const Vector<Real> & normal) {
UInt spatial_dimension = mesh.getSpatialDimension();
auto surface_dimension = spatial_dimension - 1;
const ElementType & type = element.type;
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage();
Matrix<Real> dn2ds2(surface_dimension * surface_dimension,
nb_nodes_per_element);
#define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \
ElementClass<type>::computeDN2DS2(natural_coord, dn2ds2)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL);
#undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL
Matrix<Real> coords(spatial_dimension, nb_nodes_per_element);
mesh.extractNodalValuesFromElement(positions, coords.storage(),
elem_val + element.element * nb_nodes_per_element,
nb_nodes_per_element, spatial_dimension);
Matrix<Real> curvature(spatial_dimension, surface_dimension*surface_dimension);
curvature.mul<false, true>(coords, dn2ds2);
Matrix<Real> H(surface_dimension, surface_dimension);
UInt i = 0;
for (auto alpha : arange(surface_dimension)) {
for (auto beta : arange(surface_dimension)) {
Vector<Real> temp(curvature(i));
H(alpha, beta) = temp.dot(normal);
i++;
}
}
return H;
}
}
diff --git a/src/model/contact_mechanics/geometry_utils.hh b/src/model/contact_mechanics/geometry_utils.hh
index 82a503c73..3dba0b63a 100644
--- a/src/model/contact_mechanics/geometry_utils.hh
+++ b/src/model/contact_mechanics/geometry_utils.hh
@@ -1,121 +1,125 @@
/**
* @file geometry_utils.hh
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Mon Sep 30 2019
* @date last modification: Mon Sep 30 2019
*
* @brief class to compute geometry related quantities
*
* @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 "aka_common.hh"
#include "fe_engine.hh"
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_GEOMETRY_UTILS_HH__
#define __AKANTU_GEOMETRY_UTILS_HH__
namespace akantu {
class GeometryUtils {
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// computes the normal on an element
static void normal(const Mesh & mesh, const Array<Real> & positions,
const Element & element, Vector<Real> & normal, bool outward=true);
/// computes the orthogonal projection on a set of elements and
/// returns natural projection and normal gap and index of element
static UInt orthogonalProjection(const Mesh & mesh, const Array<Real> & positions,
const Vector<Real> & slave,
const Array<Element> & elements,
Real & gap, Vector<Real> & natural_projection,
- Vector<Real> & normal, Real alpha, Real tolerance = 1e-10,
+ Vector<Real> & normal, Real alpha,
+ UInt max_iterations = 100,
+ Real tolerance = 1e-10,
Real extension_tolerance = 1e-5);
/// computes the natural projection on an element
static void naturalProjection(const Mesh & mesh, const Array<Real> & positions,
const Element & element, Vector<Real> & real_projection,
- Vector<Real> & natural_projection, Real tolerance = 1e-10);
+ Vector<Real> & natural_projection,
+ UInt max_iterations = 100,
+ Real tolerance = 1e-10);
/// computes the real projection on an element
static void realProjection(const Mesh & mesh, const Array<Real> & positions,
const Vector<Real> & slave, const Element & element,
const Vector<Real> & normal, Vector<Real> & projection);
/// computes the real projection from a natural coordinate
static void realProjection(const Mesh & mesh, const Array<Real> & positions,
const Element & element, const Vector<Real> & natural_coord,
Vector<Real> & projection);
/// computes the covariant basis/ local surface basis/ tangents on projection
/// point
static void covariantBasis(const Mesh & mesh, const Array<Real> & positions,
const Element & element, const Vector<Real> & normal,
Vector<Real> & natural_coord,
Matrix<Real> & basis);
// computes the curvature on projection
static void curvature(const Mesh & mesh, const Array<Real> & positions,
const Element & element, const Vector<Real> & natural_coord,
Matrix<Real> & curvature);
/// computes the contravariant basis on projection point
static void contravariantBasis(const Matrix<Real> & covariant,
Matrix<Real> & contravariant);
/// computes metric tesnor with covariant components
static Matrix<Real> covariantMetricTensor(const Matrix<Real> & );
/// computes metric tensor with contravariant components
static Matrix<Real> contravariantMetricTensor(const Matrix<Real> & );
// computes curvature tensor with convariant components
static Matrix<Real> covariantCurvatureTensor(const Mesh &,
const Array<Real> &,
const Element &,
const Vector<Real> &,
const Vector<Real> &);
/// checks if the element is truly a boundary element or not
inline static bool isBoundaryElement(const Mesh & mesh, const Element & element);
/// checks if the natural projection is valid for not
inline static bool isValidProjection(const Vector<Real> & projection,
Real extension_tolerance = 1e-5);
};
}
#include "geometry_utils_inline_impl.cc"
#endif /* __AKANTU_GEOMETRY_UTILS_HH__ */