diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index bdcf4561e..134954cb8 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,286 +1,290 @@ /** * @file contact_detector.cc * * @author Mohit Pundir * * @date creation: Wed Dec 05 2018 * @date last modification: Thu Jun 24 2021 * * @brief Mother class for all detection algorithms * * * @section LICENSE * * Copyright (©) 2018-2021 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) + : ContactDetector(mesh, mesh.getNodes(), id) {} + +/* -------------------------------------------------------------------------- */ +ContactDetector::ContactDetector(Mesh & mesh, Array positions, const ID & id) : Parsable(ParserType::_contact_detector, id), - AbstractContactDetector(mesh) { + AbstractContactDetector(mesh, positions) { AKANTU_DEBUG_IN(); const Parser & parser = getStaticParser(); const ParserSection & section = *(parser.getSubSections(ParserType::_contact_detector).first); this->parseSection(section); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactDetector::parseSection(const ParserSection & section) { 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) { auto 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.emplace_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(); } } // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector.hh b/src/model/contact_mechanics/contact_detector.hh index a4fd629aa..d6195644e 100644 --- a/src/model/contact_mechanics/contact_detector.hh +++ b/src/model/contact_mechanics/contact_detector.hh @@ -1,197 +1,200 @@ /** * @file contact_detector.hh * * @author Mohit Pundir * * @date creation: Mon Dec 13 2010 * @date last modification: Thu Jun 24 2021 * * @brief Mother class for all detection algorithms * * * @section LICENSE * * Copyright (©) 2010-2021 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 "aka_bbox.hh" #include "aka_common.hh" #include "aka_grid_dynamic.hh" #include "contact_detector_shared.hh" #include "contact_element.hh" #include "element_class.hh" #include "element_group.hh" #include "fe_engine.hh" #include "geometry_utils.hh" #include "mesh.hh" #include "mesh_io.hh" #include "parsable.hh" #include "surface_selector.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_DETECTOR_HH__ #define __AKANTU_CONTACT_DETECTOR_HH__ namespace akantu { enum class Surface { master, slave }; /* -------------------------------------------------------------------------- */ class ContactDetector : public Parsable, public AbstractContactDetector { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ public: ContactDetector(Mesh & /*mesh*/, const ID & id = "contact_detector"); + ContactDetector(Mesh & /*mesh*/, Array positions, + const ID & id = "contact_detector"); + ~ContactDetector() override = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ public: /// performs all search steps void search(Array & elements, Array & gaps, Array & normals, Array & tangents, Array & projections); /// performs global spatial search to construct spatial grids void globalSearch(SpatialGrid & /*slave_grid*/, SpatialGrid & /*master_grid*/); /// performs local search to find closet master node to a slave node void localSearch(SpatialGrid & /*slave_grid*/, SpatialGrid & /*master_grid*/); /// create contact elements void createContactElements(Array & elements, Array & gaps, Array & normals, Array & tangents, Array & projections); private: /// reads the input file to get contact detection options void parseSection(const ParserSection & section) override; /* ------------------------------------------------------------------------ */ /* Inline Methods */ /* ------------------------------------------------------------------------ */ public: /// checks whether the natural projection is valid or not inline bool checkValidityOfProjection(Vector & projection) const; /// computes the optimal cell size for grid inline void computeCellSpacing(Vector & spacing) const; /// constructs a grid containing nodes lying within bounding box inline void constructGrid(SpatialGrid & grid, BBox & bbox, const Array & nodes_list) const; /// constructs the bounding box based on nodes list inline void constructBoundingBox(BBox & bbox, const Array & nodes_list) const; /// computes the maximum in radius for a given mesh inline void computeMaximalDetectionDistance(); /// constructs the connectivity for a contact element inline Vector constructConnectivity(UInt & slave, const Element & master) const; /// computes normal on an element inline void computeNormalOnElement(const Element & element, Vector & normal) const; /// extracts vectors which forms the plane of element inline void vectorsAlongElement(const Element & el, Matrix & vectors) const; /// computes the gap between slave and its projection on master /// surface inline Real computeGap(const Vector & slave, const Vector & master) const; /// filter boundary elements inline void filterBoundaryElements(const Array & elements, Array & boundary_elements) const; /// 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 & slave_node, const Real & gap, const Vector & normal) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// 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(SurfaceSelector, *surface_selector, SurfaceSelector &); AKANTU_SET_MACRO(SurfaceSelector, surface_selector, std::shared_ptr); /* ------------------------------------------------------------------------ */ /* 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; /// node selector for selecting master and slave nodes std::shared_ptr surface_selector; /// contact pair slave node to closet master node std::vector> contact_pairs; /// 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/contact_detector_shared.cc b/src/model/contact_mechanics/contact_detector_shared.cc index 40c0aca0e..ffd7054fb 100644 --- a/src/model/contact_mechanics/contact_detector_shared.cc +++ b/src/model/contact_mechanics/contact_detector_shared.cc @@ -1,13 +1,11 @@ #include "contact_detector_shared.hh" namespace akantu { -AbstractContactDetector::AbstractContactDetector(akantu::Mesh & mesh) +AbstractContactDetector::AbstractContactDetector(akantu::Mesh & mesh, Array initial_positions) : mesh(mesh), spatial_dimension(mesh.getSpatialDimension()), - positions(0, spatial_dimension) { - - this->positions.copy(mesh.getNodes()); + positions(initial_positions) { } } // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector_shared.hh b/src/model/contact_mechanics/contact_detector_shared.hh index b62829d66..51d181463 100644 --- a/src/model/contact_mechanics/contact_detector_shared.hh +++ b/src/model/contact_mechanics/contact_detector_shared.hh @@ -1,103 +1,106 @@ #ifndef __AKANTU_CONTACT_DETECTOR_SHARED_HH__ #define __AKANTU_CONTACT_DETECTOR_SHARED_HH__ #include "aka_common.hh" #include "aka_grid_dynamic.hh" #include "fe_engine.hh" namespace akantu { /// Base class for contact detectors class AbstractContactDetector { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ public: - AbstractContactDetector(Mesh & mesh); + AbstractContactDetector(Mesh & mesh, Array initial_positions); + virtual ~AbstractContactDetector() = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ protected: + // TODO: move body to tmpl header file + /// Return a new vector with the positions of a node. inline Vector getNodePosition(UInt node) const { Vector position(spatial_dimension); for (UInt s : arange(spatial_dimension)) { position(s) = this->positions(node, s); } return position; } /// Fill the matrix with the coordinates of an element. inline void coordinatesOfElement(const Element & el, Matrix & coords) const { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); const Vector connect = mesh.getConnectivity(el.type, _not_ghost) .begin(nb_nodes_per_element)[el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; for (UInt s : arange(spatial_dimension)) { coords(s, n) = this->positions(node, s); } } } /// Compute the minimum and maximum element sizes. template inline auto computeElementSizes(vector_type && nodes) const { struct { Real min_size; Real max_size; } out_values; out_values = { std::numeric_limits::max(), std::numeric_limits::min() }; for (auto node : nodes) { Array elements; this->mesh.getAssociatedElements(node, elements); for (auto element : elements) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element.type); Matrix elem_coords(spatial_dimension, nb_nodes_per_element); this->coordinatesOfElement(element, elem_coords); Real elem_size = FEEngine::getElementInradius(elem_coords, element.type); out_values.max_size = std::max(out_values.max_size, elem_size); out_values.min_size = std::min(out_values.min_size, elem_size); } } return out_values; } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the mesh AKANTU_GET_MACRO(Mesh, mesh, Mesh &) AKANTU_GET_MACRO_NOT_CONST(Positions, positions, Array &); AKANTU_SET_MACRO(Positions, positions, Array); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Mesh Mesh & mesh; /// dimension of the model const UInt spatial_dimension; /// contains the updated positions of the nodes - Array positions; + Array positions; // TODO: get on the fly, storing is wasteful }; } // namespace akantu #endif // __AKANTU_CONTACT_DETECTOR_SHARED_HH__ diff --git a/src/model/contact_mechanics_internodes/contact_detector_internodes.cc b/src/model/contact_mechanics_internodes/contact_detector_internodes.cc index 3400c6b8e..74273c7a1 100644 --- a/src/model/contact_mechanics_internodes/contact_detector_internodes.cc +++ b/src/model/contact_mechanics_internodes/contact_detector_internodes.cc @@ -1,335 +1,335 @@ /** * @file contact_detector_internodes.cc * * @author Moritz Waldleben * * @date creation: Thu Jul 09 2022 * @date last modification: Thu Jul 17 2022 * * @brief Algorithm to detetect contact nodes * * * @section LICENSE * * Copyright (©) 2010-2021 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_internodes.hh" #include "contact_detector_shared.hh" #include "element_group.hh" #include "node_group.hh" #include "parsable.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactDetectorInternodes::ContactDetectorInternodes(Mesh & mesh, const ID & id) : Parsable(ParserType::_contact_detector, id), - AbstractContactDetector(mesh) { + AbstractContactDetector(mesh, mesh.getNodes()) { const Parser & parser = getStaticParser(); const ParserSection & section = *(parser.getSubSections(ParserType::_contact_detector).first); this->parseSection(section); auto & initial_master_node_group = mesh.createNodeGroup("initial_contact_master_nodes"); auto & initial_slave_node_group = mesh.createNodeGroup("initial_contact_slave_nodes"); auto & master_node_group = mesh.createNodeGroup("contact_master_nodes"); auto & slave_node_group = mesh.createNodeGroup("contact_slave_nodes"); initial_master_node_group.append( mesh.getElementGroup(id_master_nodes).getNodeGroup()); initial_slave_node_group.append(mesh.getElementGroup(id_slave_nodes).getNodeGroup()); master_node_group.append( mesh.getElementGroup(id_master_nodes).getNodeGroup()); slave_node_group.append(mesh.getElementGroup(id_slave_nodes).getNodeGroup()); } /* -------------------------------------------------------------------------- */ void ContactDetectorInternodes::parseSection(const ParserSection & section) { this->id_master_nodes = section.getParameterValue("master"); this->id_slave_nodes = section.getParameterValue("slave"); } // TODO: why are these initial node groups a thing? /* -------------------------------------------------------------------------- */ NodeGroup & ContactDetectorInternodes::getInitialMasterNodeGroup() { return mesh.getNodeGroup("initial_contact_master_nodes"); } /* -------------------------------------------------------------------------- */ NodeGroup & ContactDetectorInternodes::getInitialSlaveNodeGroup() { return mesh.getNodeGroup("initial_contact_slave_nodes"); } /* -------------------------------------------------------------------------- */ NodeGroup & ContactDetectorInternodes::getMasterNodeGroup() { return mesh.getNodeGroup("contact_master_nodes"); } /* -------------------------------------------------------------------------- */ NodeGroup & ContactDetectorInternodes::getSlaveNodeGroup() { return mesh.getNodeGroup("contact_slave_nodes"); } /* -------------------------------------------------------------------------- */ void ContactDetectorInternodes::findContactNodes() { auto & master_node_group = getMasterNodeGroup(); auto & slave_node_group = getSlaveNodeGroup(); auto element_sizes = computeElementSizes(arange(mesh.getNbNodes())); Real security_factor = 3; // TODO: allow user to pick this via material file Real grid_spacing = element_sizes.max_size * security_factor; bool still_isolated_nodes = true; while (still_isolated_nodes) { auto master_grid = constructGrid(master_node_group, grid_spacing); auto slave_grid = constructGrid(slave_node_group, grid_spacing); auto && nb_slave_nodes_inside_radius = computeRadiuses(master_radiuses, master_node_group, master_grid, slave_node_group, slave_grid); auto && nb_master_nodes_inside_radius = computeRadiuses(slave_radiuses, slave_node_group, slave_grid, master_node_group, master_grid); still_isolated_nodes = false; if (master_node_group.applyNodeFilter([&](auto && node) { return nb_master_nodes_inside_radius[node] > 0; }) > 0) { still_isolated_nodes = true; } if (slave_node_group.applyNodeFilter([&](auto && node) { return nb_slave_nodes_inside_radius[node] > 0; }) > 0) { still_isolated_nodes = true; } // TODO: These calls are useless when we're only removing nodes... master_node_group.optimize(); slave_node_group.optimize(); } } /* -------------------------------------------------------------------------- */ SpatialGrid ContactDetectorInternodes::constructGrid(const akantu::NodeGroup & node_group, akantu::Real spacing) const { SpatialGrid grid(spatial_dimension); auto spacing_vector = Vector(spatial_dimension, spacing); grid.setSpacing(spacing_vector); for (UInt node : node_group) { grid.insert(node, getNodePosition(node)); } return grid; } /* -------------------------------------------------------------------------- */ Matrix ContactDetectorInternodes::constructInterpolationMatrix( const NodeGroup & ref_node_group, const NodeGroup & eval_node_group, Array eval_radiuses) { auto && phi_eval_eval = constructPhiMatrix(eval_node_group, eval_node_group, eval_radiuses); auto && phi_ref_eval = constructPhiMatrix(ref_node_group, eval_node_group, eval_radiuses); Matrix phi_eval_eval_inv(eval_node_group.size(), eval_node_group.size(), 1.); phi_eval_eval_inv.inverse(phi_eval_eval); auto && interpol_ref_eval = phi_ref_eval * phi_eval_eval_inv; Vector ones(eval_node_group.size(), 1); ones.set(1.0); for (UInt i : arange(ref_node_group.size())) { Real row_sum = 0.; for (UInt k : arange(eval_node_group.size())) { row_sum = row_sum + interpol_ref_eval(i, k); } for (UInt j : arange(eval_node_group.size())) { interpol_ref_eval(i, j) = interpol_ref_eval(i, j) / row_sum; } } // extended to 2D Matrix interpol_ref_eval_ext(spatial_dimension*interpol_ref_eval.rows(), spatial_dimension*interpol_ref_eval.cols(), 0.); for (UInt i : arange(interpol_ref_eval.rows())) { for (UInt j : arange(interpol_ref_eval.cols())) { for (int dim = 0; dim < spatial_dimension; dim++) { interpol_ref_eval_ext(spatial_dimension*i+dim, spatial_dimension*j+dim) = interpol_ref_eval(i, j); } } } return interpol_ref_eval_ext; } /* -------------------------------------------------------------------------- */ Matrix ContactDetectorInternodes::constructPhiMatrix(const NodeGroup & ref_node_group, const NodeGroup & eval_node_group, Array & eval_radiuses) { auto && positions = mesh.getNodes(); auto && positions_view = make_view(positions, spatial_dimension).begin(); Array distances(eval_node_group.size()); Matrix phi(ref_node_group.size(), eval_node_group.size()); phi.set(0.); for (auto && ref_node_data : enumerate(ref_node_group.getNodes())) { auto i = std::get<0>(ref_node_data); auto ref_node = std::get<1>(ref_node_data); computeDistancesToRefNode(ref_node, eval_node_group, distances); for (UInt j : arange(eval_node_group.size())) { if (distances(j) <= eval_radiuses(j)) { phi(i, j) = computeRadialBasisInterpolation(distances(j), eval_radiuses(j)); } } } return phi; } /* -------------------------------------------------------------------------- */ std::map ContactDetectorInternodes::computeRadiuses( Array & attack_radiuses, const NodeGroup & ref_node_group, const SpatialGrid & ref_grid, const NodeGroup & eval_node_group, const SpatialGrid & eval_grid) { Real c = 0.5; Real C = 0.95; Real d = 0.05; // maximum number of support nodes UInt f = std::floor(1 / (std::pow(1 - c, 4) * (1 + 4 * c))); std::vector temp_nodes; attack_radiuses.resize(ref_node_group.size()); std::map nb_neighboor_nodes_inside_radiuses; std::map nb_opposite_nodes_inside_radiuses; UInt nb_iter = 0; UInt max_nb_supports = std::numeric_limits::max(); while (max_nb_supports > f && nb_iter < MAX_RADIUS_ITERATIONS) { for (auto && ref_node_data : enumerate(ref_node_group.getNodes())) { auto j = std::get<0>(ref_node_data); auto ref_node = std::get<1>(ref_node_data); // compute radius of attack, i.e. distance to closest neighbor node Real attack_radius = std::numeric_limits::max(); for (auto neighboor_node : ref_grid.setToNeighboring(getNodePosition(ref_node), temp_nodes)) { if (neighboor_node != ref_node) { Real distance = computeDistance(ref_node, neighboor_node); attack_radius = std::min(distance, attack_radius); } } Real correction_radius = std::sqrt(d * d + 0.25 * std::pow(attack_radius, 2)); correction_radius = std::max(attack_radius, correction_radius); if (correction_radius > attack_radius / c) { // TODO: shouldn't this be correction_radius * c ? attack_radius = correction_radius / c; } else { attack_radius = correction_radius; } attack_radiuses(j) = attack_radius; // compute number of neighboor nodes inside attack radius for (auto neighbor_node : eval_grid.setToNeighboring(getNodePosition(ref_node), temp_nodes)) { Real distance = computeDistance(ref_node, neighbor_node); if (distance <= attack_radius) { nb_neighboor_nodes_inside_radiuses[neighbor_node]++; } } // compute number of opposite nodes inside C * attack_radius for (auto opposite_node : eval_grid.setToNeighboring(getNodePosition(ref_node), temp_nodes)) { Real distance = computeDistance(ref_node, opposite_node); if (distance < C * attack_radius) { nb_opposite_nodes_inside_radiuses[opposite_node]++; } } } // maximum number of neighboors inside radius // aka maximum number of supports max_nb_supports = 0; for (const auto & entry : nb_neighboor_nodes_inside_radiuses) { max_nb_supports = std::max(max_nb_supports, entry.second); } // correct maximum number of support nodes if (max_nb_supports > f) { c = 0.5 * (1 + c); f = floor(1 / (pow(1 - c, 4) * (1 + 4 * c))); nb_neighboor_nodes_inside_radiuses.clear(); nb_opposite_nodes_inside_radiuses.clear(); nb_iter++; } } if (nb_iter == MAX_RADIUS_ITERATIONS) { AKANTU_EXCEPTION("Could not find suitable radii, maximum number of iterations (" << nb_iter << ") was exceeded"); } return nb_opposite_nodes_inside_radiuses; } /* -------------------------------------------------------------------------- */ void ContactDetectorInternodes::computeDistancesToRefNode( UInt & ref_node, const NodeGroup & eval_node_group, Array & out_array) { auto && positions = mesh.getNodes(); auto && positions_view = make_view(positions, spatial_dimension).begin(); Vector ref_pos = positions_view[ref_node]; for (auto && eval_node_data : enumerate(eval_node_group.getNodes())) { auto i = std::get<0>(eval_node_data); auto eval_node = std::get<1>(eval_node_data); auto && pos_eval = positions_view[eval_node]; out_array(i) = ref_pos.distance(pos_eval); } } /* -------------------------------------------------------------------------- */ Real ContactDetectorInternodes::computeRadialBasisInterpolation( const Real distance, const Real radius) { /// rescaled radial basis function: Wendland Real ratio = distance / radius; Real phi_of_x = std::pow(1 - ratio, 4) * (1 + 4 * ratio); return phi_of_x; } } // namespace akantu