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