diff --git a/src/model/contact_mechanics/geometry_utils.cc b/src/model/contact_mechanics/geometry_utils.cc index 9f3e0ddbc..d875cda48 100644 --- a/src/model/contact_mechanics/geometry_utils.cc +++ b/src/model/contact_mechanics/geometry_utils.cc @@ -1,431 +1,434 @@ /** * @file geometry_utils.cc * * @author Mohit Pundir * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "geometry_utils.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void GeometryUtils::normal(const Mesh & mesh, const Array & positions, const Element & element, Vector & normal, bool outward) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); Matrix 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 vectors(spatial_dimension, surface_dimension); switch (spatial_dimension) { case 1: { normal[0] = 1; break; } case 2: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); Math::normal2(vectors.storage(), normal.storage()); break; } case 3: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); vectors(1) = Vector(coords(2)) - Vector(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 outside(spatial_dimension); mesh.getBarycenter(element, outside); // check if mesh facets exists for cohesive elements contact Vector inside(spatial_dimension); if (mesh.isMeshFacets()) { mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); } else { mesh.getBarycenter(element_to_subelement[0], inside); } Vector 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 & positions, const Element & element, const Vector & normal, Vector & natural_coord, Matrix & 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 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 dnds(surface_dimension, nb_nodes_per_element); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(natural_coord, dnds) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL tangents.mul(dnds, nodes_coord); auto temp_tangents = tangents.transpose(); for (UInt i = 0; i < spatial_dimension - 1; ++i) { auto temp = Vector(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 e_z(3); e_z[0] = 0.; e_z[1] = 0.; e_z[2] = 1.; Vector 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(tang_trans(0)); auto tang2 = Vector(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 & positions, const Element & element, const Vector & natural_coord, Matrix & 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 dn2ds2(surface_dimension * surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ ElementClass::computeDN2DS2(natural_coord, dn2ds2) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); #undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL Matrix 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(coords, dn2ds2); } /* -------------------------------------------------------------------------- */ UInt GeometryUtils::orthogonalProjection(const Mesh & mesh, const Array & positions, const Vector & slave, const Array & elements, Real & gap, Vector & natural_projection, - Vector & normal, Real alpha, Real tolerance, + Vector & normal, Real alpha, + Real projection_tolerance, Real extension_tolerance) { UInt index = UInt(-1); Real min_gap = std::numeric_limits::max(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_same_sides{0}; UInt counter{0}; for (auto & element : elements) { if (!GeometryUtils::isBoundaryElement(mesh, element)) continue; Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, positions, element, normal_ele); Vector master(spatial_dimension); GeometryUtils::realProjection(mesh, positions, slave, element, normal_ele, master); Vector xi(natural_projection.size()); - GeometryUtils::naturalProjection(mesh, positions, element, master, xi, tolerance); + GeometryUtils::naturalProjection(mesh, positions, element, master, xi, + 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 (withon tolerance) for both + // + alpha) should be close to zero (within tolerance) for both // cases - Real tolerance = 1e-8; + Real direction_tolerance = 1e-8; auto product = master_to_slave.dot(normal_ele); auto variation = std::abs(product + alpha); - if (variation <= tolerance and + 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 <= tolerance) + if(temp_gap == 0 or variation <= direction_tolerance) nb_same_sides++; counter++; } - // if point is not onthe same side of all elements than it si not + // if point is not on the same side of all elements than it si not // consider even if the closet master element is found if(nb_same_sides != elements.size()) index = UInt(-1); return index; } /* -------------------------------------------------------------------------- */ void GeometryUtils::realProjection(const Mesh & mesh, const Array & positions, const Vector & slave, const Element & element, const Vector & normal, Vector & 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 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 point(nodes_coord(0)); Real alpha = (slave - point).dot(normal); projection = slave - alpha * normal; } /* -------------------------------------------------------------------------- */ void GeometryUtils::realProjection(const Mesh & mesh, const Array & positions, const Element & element, const Vector & natural_coord, Vector & projection) { UInt spatial_dimension = mesh.getSpatialDimension(); const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element.type); Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(natural_coord, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL Matrix 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(nodes_coord, shapes); } /* -------------------------------------------------------------------------- */ void GeometryUtils::naturalProjection(const Mesh & mesh, const Array & positions, const Element & element, Vector & real_projection, - Vector & natural_projection, Real tolerance) { + Vector & natural_projection, + 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 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::inverseMap(real_projection, nodes_coord, \ - natural_projection, tolerance) + natural_projection, projection_tolerance) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_COORDINATE); #undef GET_NATURAL_COORDINATE } /* -------------------------------------------------------------------------- */ void GeometryUtils::contravariantBasis(const Matrix & covariant, Matrix & contravariant) { auto inv_A = GeometryUtils::contravariantMetricTensor(covariant); contravariant.mul(inv_A, covariant); } /* -------------------------------------------------------------------------- */ Matrix GeometryUtils::covariantMetricTensor(const Matrix & covariant_bases) { Matrix A(covariant_bases.rows(), covariant_bases.rows()); A.mul(covariant_bases, covariant_bases); return A; } /* -------------------------------------------------------------------------- */ Matrix GeometryUtils::contravariantMetricTensor(const Matrix & covariant_bases) { auto A = GeometryUtils::covariantMetricTensor(covariant_bases); auto inv_A = A.inverse(); return inv_A; } /* -------------------------------------------------------------------------- */ Matrix GeometryUtils::covariantCurvatureTensor(const Mesh & mesh, const Array & positions, const Element & element, const Vector & natural_coord, const Vector & 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 dn2ds2(surface_dimension * surface_dimension, nb_nodes_per_element); #define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ ElementClass::computeDN2DS2(natural_coord, dn2ds2) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); #undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL Matrix 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 curvature(spatial_dimension, surface_dimension*surface_dimension); curvature.mul(coords, dn2ds2); Matrix H(surface_dimension, surface_dimension); UInt i = 0; for (auto alpha : arange(surface_dimension)) { for (auto beta : arange(surface_dimension)) { Vector temp(curvature(i)); H(alpha, beta) = temp.dot(normal); i++; } } return H; } }