Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 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/node.h"
      26             : #include "libmesh/parallel.h"
      27             : 
      28             : // C++ Includes
      29             : #include <limits>
      30             : #include <utility>
      31             : 
      32             : namespace
      33             : {
      34             : using libMesh::Real;
      35             : 
      36             : // 10 bits per coordinate, to work with 32+ bit machines
      37             : const unsigned int chunkmax = 1024;
      38             : const Real chunkfloat = 1024.0;
      39             : }
      40             : 
      41             : 
      42             : 
      43             : namespace libMesh
      44             : {
      45             : 
      46             : template <typename T>
      47          71 : void LocationMap<T>::init(MeshBase & mesh)
      48             : {
      49             :   // This function must be run on all processors at once
      50             :   // for non-serial meshes
      51          71 :   if (!mesh.is_serial())
      52           0 :     libmesh_parallel_only(mesh.comm());
      53             : 
      54           4 :   LOG_SCOPE("init()", "LocationMap");
      55             : 
      56             :   // Clear the old map
      57           2 :   _map.clear();
      58             : 
      59             :   // Cache a bounding box
      60          71 :   _lower_bound.clear();
      61          71 :   _lower_bound.resize(LIBMESH_DIM, std::numeric_limits<Real>::max());
      62          71 :   _upper_bound.clear();
      63          71 :   _upper_bound.resize(LIBMESH_DIM, -std::numeric_limits<Real>::max());
      64             : 
      65       33555 :   for (auto & node : mesh.node_ptr_range())
      66       69440 :     for (unsigned int i=0; i != LIBMESH_DIM; ++i)
      67             :       {
      68             :         // Expand the bounding box if necessary
      69       52368 :         _lower_bound[i] = std::min(_lower_bound[i],
      70       54138 :                                    (*node)(i));
      71       52080 :         _upper_bound[i] = std::max(_upper_bound[i],
      72        4116 :                                    (*node)(i));
      73             :       }
      74             : 
      75             :   // On a parallel mesh we might not yet have a full bounding box
      76          71 :   if (!mesh.is_serial())
      77             :     {
      78          64 :       mesh.comm().min(_lower_bound);
      79          64 :       mesh.comm().max(_upper_bound);
      80             :     }
      81             : 
      82          71 :   this->fill(mesh);
      83          71 : }
      84             : 
      85             : 
      86             : 
      87             : template <typename T>
      88       17360 : void LocationMap<T>::insert(T & t)
      89             : {
      90       17360 :   this->_map.emplace(this->key(this->point_of(t)), &t);
      91       17360 : }
      92             : 
      93             : 
      94             : 
      95             : template <>
      96       43848 : Point LocationMap<Node>::point_of(const Node & node) const
      97             : {
      98       43848 :   return node;
      99             : }
     100             : 
     101             : 
     102             : 
     103             : template <>
     104           0 : Point LocationMap<Elem>::point_of(const Elem & elem) const
     105             : {
     106           0 :   return elem.vertex_average();
     107             : }
     108             : 
     109             : 
     110             : 
     111             : template <typename T>
     112       13244 : T * LocationMap<T>::find(const Point & p,
     113             :                          const Real tol)
     114             : {
     115         686 :   LOG_SCOPE("find()", "LocationMap");
     116             : 
     117             :   // Look for a likely key in the multimap
     118       13244 :   unsigned int pointkey = this->key(p);
     119             : 
     120             :   // Look for the exact key first
     121       13244 :   for (const auto & pr : as_range(_map.equal_range(pointkey)))
     122       13587 :     if (p.absolute_fuzzy_equals(this->point_of(*(pr.second)), tol))
     123       13244 :       return pr.second;
     124             : 
     125             :   // Look for neighboring bins' keys next
     126           0 :   for (int xoffset = -1; xoffset != 2; ++xoffset)
     127             :     {
     128           0 :       for (int yoffset = -1; yoffset != 2; ++yoffset)
     129             :         {
     130           0 :           for (int zoffset = -1; zoffset != 2; ++zoffset)
     131             :             {
     132           0 :               auto key_pos = _map.equal_range(pointkey +
     133           0 :                                               xoffset*chunkmax*chunkmax +
     134           0 :                                               yoffset*chunkmax +
     135           0 :                                               zoffset);
     136           0 :               for (const auto & pr : as_range(key_pos))
     137           0 :                 if (p.absolute_fuzzy_equals(this->point_of(*(pr.second)), tol))
     138           0 :                   return pr.second;
     139             :             }
     140             :         }
     141             :     }
     142             : 
     143           0 :   return nullptr;
     144             : }
     145             : 
     146             : 
     147             : 
     148             : template <typename T>
     149       30604 : unsigned int LocationMap<T>::key(const Point & p)
     150             : {
     151        1029 :   Real xscaled = 0., yscaled = 0., zscaled = 0.;
     152             : 
     153       30604 :   Real deltax = _upper_bound[0] - _lower_bound[0];
     154             : 
     155       30604 :   if (std::abs(deltax) > TOLERANCE)
     156       30604 :     xscaled = (p(0) - _lower_bound[0])/deltax;
     157             : 
     158             :   // Only check y-coords if libmesh is compiled with LIBMESH_DIM>1
     159             : #if LIBMESH_DIM > 1
     160       30604 :   Real deltay = _upper_bound[1] - _lower_bound[1];
     161             : 
     162       30604 :   if (std::abs(deltay) > TOLERANCE)
     163       30604 :     yscaled = (p(1) - _lower_bound[1])/deltay;
     164             : #endif
     165             : 
     166             :   // Only check z-coords if libmesh is compiled with LIBMESH_DIM>2
     167             : #if LIBMESH_DIM > 2
     168       30604 :   Real deltaz = _upper_bound[2] - _lower_bound[2];
     169             : 
     170       30604 :   if (std::abs(deltaz) > TOLERANCE)
     171       30604 :     zscaled = (p(2) - _lower_bound[2])/deltaz;
     172             : #endif
     173             : 
     174       30604 :   unsigned int n0 = static_cast<unsigned int> (chunkfloat * xscaled),
     175       30604 :     n1 = static_cast<unsigned int> (chunkfloat * yscaled),
     176       30604 :     n2 = static_cast<unsigned int> (chunkfloat * zscaled);
     177             : 
     178       30604 :   return chunkmax*chunkmax*n0 + chunkmax*n1 + n2;
     179             : }
     180             : 
     181             : 
     182             : 
     183             : template <>
     184          71 : void LocationMap<Node>::fill(MeshBase & mesh)
     185             : {
     186             :   // Populate the nodes map
     187       33488 :   for (auto & node : mesh.node_ptr_range())
     188       17427 :     this->insert(*node);
     189          71 : }
     190             : 
     191             : 
     192             : 
     193             : template <>
     194           0 : void LocationMap<Elem>::fill(MeshBase & mesh)
     195             : {
     196             :   // Populate the elem map
     197           0 :   for (auto & elem : mesh.active_element_ptr_range())
     198           0 :     this->insert(*elem);
     199           0 : }
     200             : 
     201             : 
     202             : 
     203             : template class LIBMESH_EXPORT LocationMap<Elem>;
     204             : template class LIBMESH_EXPORT LocationMap<Node>;
     205             : 
     206             : } // namespace libMesh