libMesh
point_locator_tree.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // C++ includes
21 
22 // Local Includes
23 #include "libmesh/elem.h"
24 #include "libmesh/libmesh_logging.h"
25 #include "libmesh/mesh_base.h"
26 #include "libmesh/mesh_tools.h"
27 #include "libmesh/point_locator_tree.h"
28 #include "libmesh/tree.h"
29 
30 namespace libMesh
31 {
32 
33 
34 
35 //------------------------------------------------------------------
36 // PointLocator methods
38  const PointLocatorBase * master) :
39  PointLocatorBase (mesh,master),
40  _tree (libmesh_nullptr),
41  _element (libmesh_nullptr),
42  _out_of_mesh_mode(false),
43  _target_bin_size (200),
44  _build_type(Trees::NODES)
45 {
46  this->init(_build_type);
47 }
48 
49 
50 
52  const Trees::BuildType build_type,
53  const PointLocatorBase * master) :
54  PointLocatorBase (mesh,master),
57  _out_of_mesh_mode(false),
58  _target_bin_size (200),
59  _build_type(build_type)
60 {
61  this->init(_build_type);
62 }
63 
64 
65 
67 {
68  this->clear ();
69 }
70 
71 
72 
74 {
75  // only delete the tree when we are the master
76  if (this->_tree != libmesh_nullptr)
77  {
78  if (this->_master == libmesh_nullptr)
79  // we own the tree
80  delete this->_tree;
81  else
82  // someone else owns and therefore deletes the tree
83  this->_tree = libmesh_nullptr;
84 
85  // make sure operator () throws an assertion
86  this->_initialized = false;
87  }
88 }
89 
90 
91 
93 {
94  this->init(_build_type);
95 }
96 
97 
98 
100 {
101  libmesh_assert (!this->_tree);
102 
103  if (this->_initialized)
104  {
105  // Warn that we are already initialized
106  libMesh::err << "Warning: PointLocatorTree already initialized! Will ignore this call..." << std::endl;
107 
108  // Further warn if we try to init() again with a different build_type
109  if (_build_type != build_type)
110  {
111  libMesh::err << "Warning: PointLocatorTree is using build_type = " << _build_type << ".\n"
112  << "Your requested build_type, " << build_type << " will not be used!" << std::endl;
113  }
114  }
115 
116  else
117  {
118  // Let the requested build_type override the _build_type we were
119  // constructed with. This is no big deal since we have not been
120  // initialized before.
121  _build_type = build_type;
122 
123  if (this->_master == libmesh_nullptr)
124  {
125  LOG_SCOPE("init(no master)", "PointLocatorTree");
126 
127  if (this->_mesh.mesh_dimension() == 3)
129  else
130  {
131  // A 1D/2D mesh in 3D space needs special consideration.
132  // If the mesh is planar XY, we want to build a QuadTree
133  // to search efficiently. If the mesh is truly a manifold,
134  // then we need an octree
135 #if LIBMESH_DIM > 2
136  bool is_planar_xy = false;
137 
138  // Build the bounding box for the mesh. If the delta-z bound is
139  // negligibly small then we can use a quadtree.
140  {
142 
143  const Real
144  Dx = bbox.second(0) - bbox.first(0),
145  Dz = bbox.second(2) - bbox.first(2);
146 
147  if (std::abs(Dz/(Dx + 1.e-20)) < 1e-10)
148  is_planar_xy = true;
149  }
150 
151  if (!is_planar_xy)
153  else
154 #endif
155 #if LIBMESH_DIM > 1
157 #else
159 #endif
160  }
161  }
162 
163  else
164  {
165  // We are _not_ the master. Let our Tree point to
166  // the master's tree. But for this we first transform
167  // the master in a state for which we are friends.
168  // And make sure the master has a tree!
169  const PointLocatorTree * my_master =
170  cast_ptr<const PointLocatorTree *>(this->_master);
171 
172  if (my_master->initialized())
173  this->_tree = my_master->_tree;
174  else
175  libmesh_error_msg("ERROR: Initialize master first, then servants!");
176  }
177 
178  // Not all PointLocators may own a tree, but all of them
179  // use their own element pointer. Let the element pointer
180  // be unique for every interpolator.
181  // Suppose the interpolators are used concurrently
182  // at different locations in the mesh, then it makes quite
183  // sense to have unique start elements.
184  this->_element = libmesh_nullptr;
185  }
186 
187  // ready for take-off
188  this->_initialized = true;
189 }
190 
191 
192 
194  const std::set<subdomain_id_type> * allowed_subdomains) const
195 {
197 
198  LOG_SCOPE("operator()", "PointLocatorTree");
199 
200  // If we're provided with an allowed_subdomains list and have a cached element, make sure it complies
201  if (allowed_subdomains && this->_element && !allowed_subdomains->count(this->_element->subdomain_id())) this->_element = libmesh_nullptr;
202 
203  // First check the element from last time before asking the tree
204  if (this->_element==libmesh_nullptr || !(this->_element->contains_point(p)))
205  {
206  // ask the tree
207  this->_element = this->_tree->find_element (p,allowed_subdomains);
208 
209  if (this->_element == libmesh_nullptr)
210  {
211  // If we haven't found the element, we may want to do a linear
212  // search using a tolerance.
214  {
215  if (_verbose)
216  {
217  libMesh::out << "Performing linear search using close-to-point tolerance "
219  << std::endl;
220  }
221 
222  this->_element =
223  this->perform_linear_search(p,
224  allowed_subdomains,
225  /*use_close_to_point*/ true,
227 
228  return this->_element;
229  }
230 
231  // No element seems to contain this point. In theory, our
232  // tree now correctly handles curved elements. In
233  // out-of-mesh mode this is sometimes expected, and we can
234  // just return NULL without searching further. Out of
235  // out-of-mesh mode, something must have gone wrong.
236  libmesh_assert_equal_to (_out_of_mesh_mode, true);
237 
238  return this->_element;
239  }
240  }
241 
242  // If we found an element, it should be active
243  libmesh_assert (!this->_element || this->_element->active());
244 
245  // If we found an element and have a restriction list, they better match
246  libmesh_assert (!this->_element || !allowed_subdomains || allowed_subdomains->count(this->_element->subdomain_id()));
247 
248  // return the element
249  return this->_element;
250 }
251 
252 
254  std::set<const Elem *> & candidate_elements,
255  const std::set<subdomain_id_type> * allowed_subdomains) const
256 {
258 
259  LOG_SCOPE("operator() - Version 2", "PointLocatorTree");
260 
261  // forward call to perform_linear_search
262  candidate_elements = this->perform_fuzzy_linear_search(p, allowed_subdomains, _close_to_point_tol);
263 }
264 
265 
266 
268  const std::set<subdomain_id_type> * allowed_subdomains,
269  bool use_close_to_point,
270  Real close_to_point_tolerance) const
271 {
272  LOG_SCOPE("perform_linear_search", "PointLocatorTree");
273 
274  // The type of iterator depends on the Trees::BuildType
275  // used for this PointLocator. If it's
276  // TREE_LOCAL_ELEMENTS, we only want to double check
277  // local elements during this linear search.
282 
283  for (const auto & elem : r)
284  {
285  if (!allowed_subdomains ||
286  allowed_subdomains->count(elem->subdomain_id()))
287  {
288  if (!use_close_to_point)
289  {
290  if (elem->contains_point(p))
291  return elem;
292  }
293  else
294  {
295  if (elem->close_to_point(p, close_to_point_tolerance))
296  return elem;
297  }
298  }
299  }
300 
301  return libmesh_nullptr;
302 }
303 
304 
305 std::set<const Elem *> PointLocatorTree::perform_fuzzy_linear_search(const Point & p,
306  const std::set<subdomain_id_type> * allowed_subdomains,
307  Real close_to_point_tolerance) const
308 {
309  LOG_SCOPE("perform_fuzzy_linear_search", "PointLocatorTree");
310 
311  std::set<const Elem *> candidate_elements;
312 
313  // The type of iterator depends on the Trees::BuildType
314  // used for this PointLocator. If it's
315  // TREE_LOCAL_ELEMENTS, we only want to double check
316  // local elements during this linear search.
321 
322  for (const auto & elem : r)
323  if ((!allowed_subdomains || allowed_subdomains->count(elem->subdomain_id())) && elem->close_to_point(p, close_to_point_tolerance))
324  candidate_elements.insert(elem);
325 
326  return candidate_elements;
327 }
328 
329 
330 
332 {
333  // Out-of-mesh mode should now work properly even on meshes with
334  // non-affine elements.
335  _out_of_mesh_mode = true;
336 }
337 
338 
340 {
341  _out_of_mesh_mode = false;
342 }
343 
344 
345 void PointLocatorTree::set_target_bin_size (unsigned int target_bin_size)
346 {
347  _target_bin_size = target_bin_size;
348 }
349 
350 
352 {
353  return _target_bin_size;
354 }
355 
356 
357 } // namespace libMesh
unsigned int get_target_bin_size() const
Get the target bin size.
OStreamProxy err
The SimpleRange templated class is intended to make it easy to construct ranges from pairs of iterato...
Definition: simple_range.h:36
double abs(double a)
bool active() const
Definition: elem.h:2257
virtual void init() libmesh_override
Initializes the locator, so that the operator() methods can be used.
This is a point locator.
bool _use_close_to_point_tol
true if we will use a user-specified tolerance for locating the element.
This class defines a tree that may be used for fast point location in space.
Definition: tree.h:44
libMesh::BoundingBox create_bounding_box(const MeshBase &mesh)
The same functionality as the deprecated MeshTools::bounding_box().
Definition: mesh_tools.C:329
Real _close_to_point_tol
The tolerance to use.
Tree< 8 > OctTree
An OctTree is a tree appropriate for 3D meshes.
Definition: tree.h:134
std::set< const Elem * > perform_fuzzy_linear_search(const Point &p, const std::set< subdomain_id_type > *allowed_subdomains, Real close_to_point_tolerance=TOLERANCE) const
A method to check if "fat" point p is in multiple elements.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
bool _initialized
true when properly initialized, false otherwise.
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
virtual const Elem * operator()(const Point &p, const std::set< subdomain_id_type > *allowed_subdomains=libmesh_nullptr) const libmesh_override
Locates the element in which the point with global coordinates p is located, optionally restricted to...
The libMesh namespace provides an interface to certain functionality in the library.
virtual const Elem * find_element(const Point &p, const std::set< subdomain_id_type > *allowed_subdomains=libmesh_nullptr, Real relative_tol=TOLERANCE) const =0
PointLocatorTree(const MeshBase &mesh, const PointLocatorBase *master=libmesh_nullptr)
Constructor.
const Elem * _element
Pointer to the last element that was found by the tree.
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
Tree< 4 > QuadTree
A QuadTree is a tree appropriate for 2D meshes.
Definition: tree.h:128
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
virtual void enable_out_of_mesh_mode() libmesh_override
Enables out-of-mesh mode.
const MeshBase & _mesh
constant reference to the mesh in which the point is looked for.
This is the base class for point locators.
virtual void disable_out_of_mesh_mode() libmesh_override
Disables out-of-mesh mode (default).
Trees::BuildType _build_type
How the underlying tree is built.
unsigned int _target_bin_size
Target bin size, which gets passed to the constructor of _tree.
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
bool _out_of_mesh_mode
true if out-of-mesh mode is enabled.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool _verbose
Boolean flag to indicate whether to print out extra info.
const Elem * perform_linear_search(const Point &p, const std::set< subdomain_id_type > *allowed_subdomains, bool use_close_to_point, Real close_to_point_tolerance=TOLERANCE) const
As a fallback option, it&#39;s helpful to be able to do a linear search over the entire mesh...
OStreamProxy out
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
BuildType
enum defining how to build the tree.
Definition: tree_base.h:55
const PointLocatorBase * _master
Const pointer to our master, initialized to NULL if none given.
void set_target_bin_size(unsigned int target)
Set the target bin size.
virtual void clear() libmesh_override
Clears the locator.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
Tree< 2 > BinaryTree
A BinaryTree is a tree appropriate for 1D meshes.
Definition: tree.h:122
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const
Definition: elem.C:2448
TreeBase * _tree
Pointer to our tree.