diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index eefb3ff3a..25a89c55d 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,360 +1,362 @@ /** * @file contact_detector.cc * * @author Mohit Pundir * * @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 . * */ /* -------------------------------------------------------------------------- */ #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 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("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("projection_tolerance"); this->max_iterations = section.getParameterValue("max_iterations"); this->extension_tolerance = section.getParameterValue("extension_tolerance"); } /* -------------------------------------------------------------------------- */ void ContactDetector::search(Array & elements, Array & gaps, Array & normals, Array & tangents, Array & projections) { UInt surface_dimension = spatial_dimension - 1; this->mesh.fillNodesToElements(surface_dimension); this->computeMaximalDetectionDistance(); contact_pairs.clear(); SpatialGrid master_grid(spatial_dimension); SpatialGrid slave_grid(spatial_dimension); this->globalSearch(slave_grid, master_grid); this->localSearch(slave_grid, master_grid); this->createContactElements(elements, gaps, normals, tangents, projections); } /* -------------------------------------------------------------------------- */ void ContactDetector::globalSearch(SpatialGrid & slave_grid, SpatialGrid & 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 center(spatial_dimension); bbox_intersection.getCenter(center); Vector 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 & slave_grid, SpatialGrid & 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 pos(spatial_dimension); for (UInt s : arange(spatial_dimension)) pos(s) = this->positions(slave_node, s); Real closet_distance = std::numeric_limits::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 elements; this->mesh.getAssociatedElements(slave_node, elements); for (auto & elem : elements) { if (elem.kind() != _ek_regular) continue; Vector connectivity = const_cast(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 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 & contact_elements, Array & gaps, Array & normals, Array & tangents, Array & 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 slave(spatial_dimension); for (UInt s : arange(spatial_dimension)) slave(s) = this->positions(slave_node, s); const auto & master_node = pairs.second; Array elements; this->mesh.getAssociatedElements(master_node, elements); auto & gap = gaps.begin()[slave_node]; Vector normal(normals.begin(spatial_dimension)[slave_node]); Vector projection(projections.begin(surface_dimension)[slave_node]); Matrix tangent(tangents.begin(surface_dimension, spatial_dimension)[slave_node]); auto index = GeometryUtils::orthogonalProjection(mesh, positions, slave, elements, gap, projection, normal, tangent, 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; + projection *= 0; + tangent *= 0; continue; } // create contact element contact_elements.push_back(ContactElement(slave_node, elements[index])); } contact_pairs.clear(); } /* -------------------------------------------------------------------------- */ /*void ContactDetector::createContactElements(Array & contact_elements, Array & gaps, Array & normals, Array & 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 slave(spatial_dimension); for (UInt s : arange(spatial_dimension)) slave(s) = this->positions(slave_node, s); const auto & master_node = pairs.second; Array elements; this->mesh.getAssociatedElements(master_node, elements); auto & gap = gaps.begin()[slave_node]; Vector normal(normals.begin(spatial_dimension)[slave_node]); Vector 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/geometry_utils.cc b/src/model/contact_mechanics/geometry_utils.cc index 4751dfeca..2ed946fb7 100644 --- a/src/model/contact_mechanics/geometry_utils.cc +++ b/src/model/contact_mechanics/geometry_utils.cc @@ -1,847 +1,847 @@ /** * @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::normal(const Mesh & mesh, const Element & element, Matrix & tangents, Vector & normal, bool outward) { UInt spatial_dimension = mesh.getSpatialDimension(); // to ensure that normal is always outwards from master surface we // compute a direction vector form inside of element attached to the // suurface element Vector inside_to_outside(spatial_dimension); 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); } inside_to_outside = outside - inside; //auto projection = inside_to_outside.dot(normal); //if (projection < 0) { // normal *= -1.0; //} } // to ensure that direction of tangents are correct, cross product // of tangents should give be in the same direction as of inside to outside switch (spatial_dimension) { case 2: { normal[0] = -tangents(0, 1); normal[1] = tangents(0, 0); auto ddot = inside_to_outside.dot(normal); if (ddot < 0) { tangents *= -1.0; normal *= -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); normal = tang1_cross_tang2 / tang1_cross_tang2.norm(); auto ddot = inside_to_outside.dot(normal); if (ddot < 0) { tang_trans(1) *= -1.0; normal *= -1.0; } tangents = tang_trans.transpose(); break; } default: break; } } /* -------------------------------------------------------------------------- */ 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::covariantBasis(const Mesh & mesh, const Array & positions, const Element & element, 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(); } /* -------------------------------------------------------------------------- */ 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, UInt max_iterations, Real projection_tolerance, Real extension_tolerance) { UInt index = UInt(-1); Real min_gap = std::numeric_limits::max(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; UInt nb_same_sides{0}; UInt nb_boundary_elements{0}; UInt counter{0}; auto & contact_group = mesh.getElementGroup("contact_surface"); for (auto & element : elements) { // filter out elements which are not there in the element group // contact surface created by the surface selector and is stored // in the mesh or mesh_facet, if a element is not there it // returnas UInt(-1) auto & elements_of_type = contact_group.getElements(element.type); if(elements_of_type.find(element.element) == UInt(-1)) continue; //if (!GeometryUtils::isBoundaryElement(mesh, element)) // continue; nb_boundary_elements++; // find the natural coordinate corresponding to the minimum gap // between slave node and master element /* 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, max_iterations, projection_tolerance);*/ Vector master(spatial_dimension); Vector xi(natural_projection.size()); GeometryUtils::naturalProjection(mesh, positions, element, slave, master, xi, max_iterations, projection_tolerance); Matrix tangent_ele(surface_dimension, spatial_dimension); GeometryUtils::covariantBasis(mesh, positions, element, xi, tangent_ele); Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, element, tangent_ele, normal_ele); // 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; } /* -------------------------------------------------------------------------- */ UInt GeometryUtils::orthogonalProjection(const Mesh & mesh, const Array & positions, const Vector & slave, const Array & elements, Real & gap, Vector & natural_projection, Vector & normal, Matrix & tangent, Real alpha, UInt max_iterations, Real projection_tolerance, Real extension_tolerance) { UInt index = UInt(-1); Real min_gap = std::numeric_limits::max(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; UInt nb_same_sides{0}; UInt nb_boundary_elements{0}; auto & contact_group = mesh.getElementGroup("contact_surface"); UInt counter{0}; for (auto & element : elements) { // filter out elements which are not there in the element group // contact surface created by the surface selector and is stored // in the mesh or mesh_facet, if a element is not there it // returnas UInt(-1) auto & elements_of_type = contact_group.getElements(element.type); if(elements_of_type.find(element.element) == UInt(-1)) continue; //if (!GeometryUtils::isBoundaryElement(mesh, element)) // continue; nb_boundary_elements++; // find the natural coordinate corresponding to the minimum gap // between slave node and master element /* 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, max_iterations, projection_tolerance);*/ Vector master(spatial_dimension); Vector xi_ele(natural_projection.size()); GeometryUtils::naturalProjection(mesh, positions, element, slave, master, xi_ele, max_iterations, projection_tolerance); Matrix tangent_ele(surface_dimension, spatial_dimension); GeometryUtils::covariantBasis(mesh, positions, element, xi_ele, tangent_ele); Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, element, tangent_ele, normal_ele); // 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_ele, extension_tolerance)) { gap = -temp_gap; min_gap = temp_gap; index = counter; natural_projection = xi_ele; normal = normal_ele; tangent = tangent_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); + //if(nb_same_sides != nb_boundary_elements) + // 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, const Vector & slave_coords, Vector & master_coords, Vector & natural_projection, UInt max_iterations, Real projection_tolerance) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt 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 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); // initial guess natural_projection.clear(); // obhjective function computed on the natural_guess Matrix f(surface_dimension, 1); // jacobian matrix computed on the natural_guess Matrix J(surface_dimension, surface_dimension); // Jinv = J^{-1} Matrix Jinv(surface_dimension, surface_dimension); // dxi = \xi_{k+1} - \xi_{k} in the iterative process Matrix dxi(surface_dimension, 1); // gradient at natural projection Matrix gradient(surface_dimension, spatial_dimension); // second derivative at natural peojection Matrix double_gradient(surface_dimension, surface_dimension); // second derivative of shape function at natural projection Matrix dn2ds2(surface_dimension * surface_dimension, nb_nodes_per_element); auto compute_double_gradient = [&double_gradient, &dn2ds2, &nodes_coord, surface_dimension, spatial_dimension](UInt & alpha, UInt & beta) { auto index = alpha * surface_dimension + beta; Vector d_alpha_beta(spatial_dimension); auto dn2ds2_transpose = dn2ds2.transpose(); Vector dn2ds2_alpha_beta(dn2ds2_transpose(index)); d_alpha_beta.mul(nodes_coord, dn2ds2_alpha_beta); return d_alpha_beta; }; /* --------------------------- */ /* init before iteration loop */ /* --------------------------- */ // do interpolation auto update_f = [&f, &master_coords, &natural_projection, &nodes_coord, &slave_coords, &gradient, surface_dimension, spatial_dimension, nb_nodes_per_element, type]() { // compute real coords on natural projection Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(natural_projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL master_coords.mul(nodes_coord, shapes); auto distance = slave_coords - master_coords; // first derivative of shape function at natural projection Matrix dnds(surface_dimension, nb_nodes_per_element); // computing gradient at projection point #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(natural_projection, dnds) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL gradient.mul(dnds, nodes_coord); // gradient transpose at natural projection Matrix gradient_transpose(surface_dimension, spatial_dimension); gradient_transpose = gradient.transpose(); // loop over surface dimensions for(auto alpha : arange(surface_dimension)) { Vector gradient_alpha(gradient_transpose(alpha)); f(alpha, 0) = -2.*gradient_alpha.dot(distance); } // compute initial error auto error = f.norm(); return error; }; auto projection_error = update_f(); /* --------------------------- */ /* iteration loop */ /* --------------------------- */ UInt iterations{0}; while(projection_tolerance < projection_error and iterations < max_iterations) { // compute covariant components of metric tensor auto a = GeometryUtils::covariantMetricTensor(gradient); // computing second derivative at natural projection #define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ ElementClass::computeDN2DS2(natural_projection, dn2ds2) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); #undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL // real coord - physical guess auto distance = slave_coords - master_coords; // computing Jacobian J for(auto alpha : arange(surface_dimension)) { for(auto beta : arange(surface_dimension)) { auto dgrad_alpha_beta = compute_double_gradient(alpha, beta); J(alpha, beta) = 2.*(a(alpha, beta) - dgrad_alpha_beta.dot(distance)); } } Jinv.inverse(J); // compute increment dxi.mul(Jinv, f, -1.0); // update our guess natural_projection += Vector(dxi(0)); projection_error = update_f(); iterations++; } /* #define GET_NATURAL_COORDINATE(type) \ ElementClass::inverseMap(slave_coords, nodes_coord, \ natural_projection, max_iterations, 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; } }