Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : 
      20             : // Local Includes
      21             : #include "libmesh/elem.h"
      22             : #include "libmesh/libmesh_logging.h"
      23             : #include "libmesh/location_maps.h"
      24             : #include "libmesh/mesh_base.h"
      25             : #include "libmesh/mesh_tools.h"
      26             : #include "libmesh/node.h"
      27             : #include "libmesh/parallel.h"
      28             : 
      29             : // C++ Includes
      30             : #include <limits>
      31             : #include <utility>
      32             : 
      33             : namespace
      34             : {
      35             : using libMesh::Real;
      36             : 
      37             : // 10 bits per coordinate, to work with 32+ bit machines
      38             : const unsigned int chunkmax = 1024;
      39             : const Real chunkfloat = 1024.0;
      40             : }
      41             : 
      42             : 
      43             : 
      44             : namespace libMesh
      45             : {
      46             : 
      47             : template <typename T>
      48          71 : void LocationMap<T>::init(MeshBase & mesh)
      49             : {
      50             :   // This function must be run on all processors at once
      51             :   // for non-serial meshes
      52          71 :   if (!mesh.is_serial())
      53           0 :     libmesh_parallel_only(mesh.comm());
      54             : 
      55           6 :   LOG_SCOPE("init()", "LocationMap");
      56             : 
      57             :   // Clear the old map
      58           2 :   _map.clear();
      59             : 
      60             :   // Cache a bounding box
      61          71 :   _bbox = MeshTools::create_nodal_bounding_box(mesh);
      62             : 
      63          71 :   this->fill(mesh);
      64          71 : }
      65             : 
      66             : 
      67             : 
      68             : template <typename T>
      69       17360 : void LocationMap<T>::insert(T & t)
      70             : {
      71       17360 :   this->_map.emplace(this->key(this->point_of(t)), &t);
      72       17360 : }
      73             : 
      74             : 
      75             : 
      76             : template <>
      77       43848 : Point LocationMap<Node>::point_of(const Node & node) const
      78             : {
      79       43848 :   return node;
      80             : }
      81             : 
      82             : 
      83             : 
      84             : template <>
      85           0 : Point LocationMap<Elem>::point_of(const Elem & elem) const
      86             : {
      87           0 :   return elem.vertex_average();
      88             : }
      89             : 
      90             : 
      91             : 
      92             : template <typename T>
      93       13244 : T * LocationMap<T>::find(const Point & p,
      94             :                          const Real tol)
      95             : {
      96        1029 :   LOG_SCOPE("find()", "LocationMap");
      97             : 
      98             :   // Look for a likely key in the multimap
      99       13244 :   unsigned int pointkey = this->key(p);
     100             : 
     101             :   // Look for the exact key first
     102       13244 :   for (const auto & pr : as_range(_map.equal_range(pointkey)))
     103       13587 :     if (p.absolute_fuzzy_equals(this->point_of(*(pr.second)), tol))
     104       13244 :       return pr.second;
     105             : 
     106             :   // Look for neighboring bins' keys next
     107           0 :   for (int xoffset = -1; xoffset != 2; ++xoffset)
     108             :     {
     109           0 :       for (int yoffset = -1; yoffset != 2; ++yoffset)
     110             :         {
     111           0 :           for (int zoffset = -1; zoffset != 2; ++zoffset)
     112             :             {
     113           0 :               auto key_pos = _map.equal_range(pointkey +
     114           0 :                                               xoffset*chunkmax*chunkmax +
     115           0 :                                               yoffset*chunkmax +
     116           0 :                                               zoffset);
     117           0 :               for (const auto & pr : as_range(key_pos))
     118           0 :                 if (p.absolute_fuzzy_equals(this->point_of(*(pr.second)), tol))
     119           0 :                   return pr.second;
     120             :             }
     121             :         }
     122             :     }
     123             : 
     124           0 :   return nullptr;
     125             : }
     126             : 
     127             : 
     128             : 
     129             : template <typename T>
     130       30604 : unsigned int LocationMap<T>::key(const Point & p)
     131             : {
     132        1029 :   Real xscaled = 0., yscaled = 0., zscaled = 0.;
     133             : 
     134       30604 :   const Real deltax = _bbox.max()(0) - _bbox.min()(0);
     135             : 
     136       30604 :   if (std::abs(deltax) > TOLERANCE)
     137       30604 :     xscaled = (p(0) - _bbox.min()(0))/deltax;
     138             : 
     139             :   // Only check y-coords if libmesh is compiled with LIBMESH_DIM>1
     140             : #if LIBMESH_DIM > 1
     141       30604 :   const Real deltay = _bbox.max()(1) - _bbox.min()(1);
     142             : 
     143       30604 :   if (std::abs(deltay) > TOLERANCE)
     144       30604 :     yscaled = (p(1) - _bbox.min()(1))/deltay;
     145             : #endif
     146             : 
     147             :   // Only check z-coords if libmesh is compiled with LIBMESH_DIM>2
     148             : #if LIBMESH_DIM > 2
     149       30604 :   const Real deltaz = _bbox.max()(2) - _bbox.min()(2);
     150             : 
     151       30604 :   if (std::abs(deltaz) > TOLERANCE)
     152       30604 :     zscaled = (p(2) - _bbox.min()(2))/deltaz;
     153             : #endif
     154             : 
     155       30604 :   unsigned int n0 = static_cast<unsigned int> (chunkfloat * xscaled),
     156       30604 :     n1 = static_cast<unsigned int> (chunkfloat * yscaled),
     157       30604 :     n2 = static_cast<unsigned int> (chunkfloat * zscaled);
     158             : 
     159       30604 :   return chunkmax*chunkmax*n0 + chunkmax*n1 + n2;
     160             : }
     161             : 
     162             : 
     163             : 
     164             : template <>
     165          71 : void LocationMap<Node>::fill(MeshBase & mesh)
     166             : {
     167             :   // Populate the nodes map
     168       33488 :   for (auto & node : mesh.node_ptr_range())
     169       17427 :     this->insert(*node);
     170          71 : }
     171             : 
     172             : 
     173             : 
     174             : template <>
     175           0 : void LocationMap<Elem>::fill(MeshBase & mesh)
     176             : {
     177             :   // Populate the elem map
     178           0 :   for (auto & elem : mesh.active_element_ptr_range())
     179           0 :     this->insert(*elem);
     180           0 : }
     181             : 
     182             : 
     183             : 
     184             : template class LIBMESH_EXPORT LocationMap<Elem>;
     185             : template class LIBMESH_EXPORT LocationMap<Node>;
     186             : 
     187             : } // namespace libMesh