libMesh
location_maps.C
Go to the documentation of this file.
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>
48 {
49  // This function must be run on all processors at once
50  // for non-serial meshes
51  if (!mesh.is_serial())
52  libmesh_parallel_only(mesh.comm());
53 
54  LOG_SCOPE("init()", "LocationMap");
55 
56  // Clear the old map
57  _map.clear();
58 
59  // Cache a bounding box
60  _lower_bound.clear();
61  _lower_bound.resize(LIBMESH_DIM, std::numeric_limits<Real>::max());
62  _upper_bound.clear();
63  _upper_bound.resize(LIBMESH_DIM, -std::numeric_limits<Real>::max());
64 
65  for (auto & node : mesh.node_ptr_range())
66  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
67  {
68  // Expand the bounding box if necessary
69  _lower_bound[i] = std::min(_lower_bound[i],
70  (*node)(i));
71  _upper_bound[i] = std::max(_upper_bound[i],
72  (*node)(i));
73  }
74 
75  // On a parallel mesh we might not yet have a full bounding box
76  if (!mesh.is_serial())
77  {
78  mesh.comm().min(_lower_bound);
79  mesh.comm().max(_upper_bound);
80  }
81 
82  this->fill(mesh);
83 }
84 
85 
86 
87 template <typename T>
89 {
90  this->_map.emplace(this->key(this->point_of(t)), &t);
91 }
92 
93 
94 
95 template <>
97 {
98  return node;
99 }
100 
101 
102 
103 template <>
105 {
106  return elem.vertex_average();
107 }
108 
109 
110 
111 template <typename T>
113  const Real tol)
114 {
115  LOG_SCOPE("find()", "LocationMap");
116 
117  // Look for a likely key in the multimap
118  unsigned int pointkey = this->key(p);
119 
120  // Look for the exact key first
121  for (const auto & pr : as_range(_map.equal_range(pointkey)))
122  if (p.absolute_fuzzy_equals(this->point_of(*(pr.second)), tol))
123  return pr.second;
124 
125  // Look for neighboring bins' keys next
126  for (int xoffset = -1; xoffset != 2; ++xoffset)
127  {
128  for (int yoffset = -1; yoffset != 2; ++yoffset)
129  {
130  for (int zoffset = -1; zoffset != 2; ++zoffset)
131  {
132  auto key_pos = _map.equal_range(pointkey +
133  xoffset*chunkmax*chunkmax +
134  yoffset*chunkmax +
135  zoffset);
136  for (const auto & pr : as_range(key_pos))
137  if (p.absolute_fuzzy_equals(this->point_of(*(pr.second)), tol))
138  return pr.second;
139  }
140  }
141  }
142 
143  return nullptr;
144 }
145 
146 
147 
148 template <typename T>
149 unsigned int LocationMap<T>::key(const Point & p)
150 {
151  Real xscaled = 0., yscaled = 0., zscaled = 0.;
152 
153  Real deltax = _upper_bound[0] - _lower_bound[0];
154 
155  if (std::abs(deltax) > TOLERANCE)
156  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  Real deltay = _upper_bound[1] - _lower_bound[1];
161 
162  if (std::abs(deltay) > TOLERANCE)
163  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  Real deltaz = _upper_bound[2] - _lower_bound[2];
169 
170  if (std::abs(deltaz) > TOLERANCE)
171  zscaled = (p(2) - _lower_bound[2])/deltaz;
172 #endif
173 
174  unsigned int n0 = static_cast<unsigned int> (chunkfloat * xscaled),
175  n1 = static_cast<unsigned int> (chunkfloat * yscaled),
176  n2 = static_cast<unsigned int> (chunkfloat * zscaled);
177 
178  return chunkmax*chunkmax*n0 + chunkmax*n1 + n2;
179 }
180 
181 
182 
183 template <>
185 {
186  // Populate the nodes map
187  for (auto & node : mesh.node_ptr_range())
188  this->insert(*node);
189 }
190 
191 
192 
193 template <>
195 {
196  // Populate the elem map
197  for (auto & elem : mesh.active_element_ptr_range())
198  this->insert(*elem);
199 }
200 
201 
202 
203 template class LIBMESH_EXPORT LocationMap<Elem>;
204 template class LIBMESH_EXPORT LocationMap<Node>;
205 
206 } // namespace libMesh
T * find(const Point &, const Real tol=TOLERANCE)
A Node is like a Point, but with more information.
Definition: node.h:52
static constexpr Real TOLERANCE
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
const Parallel::Communicator & comm() const
Point point_of(const T &) const
The libMesh namespace provides an interface to certain functionality in the library.
This is the MeshBase class.
Definition: mesh_base.h:75
virtual bool is_serial() const
Definition: mesh_base.h:211
void min(const T &r, T &o, Request &req) const
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
bool absolute_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:972
Data structures that enable location-based lookups The key is a hash of the Point location...
Definition: location_maps.h:53
void fill(MeshBase &)
unsigned int key(const Point &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void max(const T &r, T &o, Request &req) const
void init(MeshBase &)
Definition: location_maps.C:47
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
Point vertex_average() const
Definition: elem.C:688