point_locator_tree.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
25 #include "libmesh/mesh_base.h"
26 #include "libmesh/mesh_tools.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 (NULL),
41  _element (NULL),
42  _out_of_mesh_mode(false)
43 {
44  this->init(Trees::NODES);
45 }
46 
47 
48 
49 
51  const Trees::BuildType build_type,
52  const PointLocatorBase* master) :
53  PointLocatorBase (mesh,master),
54  _tree (NULL),
55  _element (NULL),
56  _out_of_mesh_mode(false)
57 {
58  this->init(build_type);
59 }
60 
61 
62 
63 
65 {
66  this->clear ();
67 }
68 
69 
70 
71 
73 {
74  // only delete the tree when we are the master
75  if (this->_tree != NULL)
76  {
77  if (this->_master == NULL)
78  // we own the tree
79  delete this->_tree;
80  else
81  // someone else owns and therefore deletes the tree
82  this->_tree = NULL;
83  }
84 }
85 
86 
87 
88 
89 
91 {
92  libmesh_assert (!this->_tree);
93 
94  if (this->_initialized)
95  {
96  libMesh::err << "ERROR: Already initialized! Will ignore this call..."
97  << std::endl;
98  }
99 
100  else
101 
102  {
103 
104  if (this->_master == NULL)
105  {
106  START_LOG("init(no master)", "PointLocatorTree");
107 
108  if (this->_mesh.mesh_dimension() == 3)
109  _tree = new Trees::OctTree (this->_mesh, 200, build_type);
110  else
111  {
112  // A 1D/2D mesh in 3D space needs special consideration.
113  // If the mesh is planar XY, we want to build a QuadTree
114  // to search efficiently. If the mesh is truly a manifold,
115  // then we need an octree
116 #if LIBMESH_DIM > 2
117  bool is_planar_xy = false;
118 
119  // Build the bounding box for the mesh. If the delta-z bound is
120  // negligibly small then we can use a quadtree.
121  {
123 
124  const Real
125  Dx = bbox.second(0) - bbox.first(0),
126  Dz = bbox.second(2) - bbox.first(2);
127 
128  if (std::abs(Dz/(Dx + 1.e-20)) < 1e-10)
129  is_planar_xy = true;
130  }
131 
132  if (!is_planar_xy)
133  _tree = new Trees::OctTree (this->_mesh, 200, build_type);
134  else
135 #endif
136 #if LIBMESH_DIM > 1
137  _tree = new Trees::QuadTree (this->_mesh, 200, build_type);
138 #else
139  _tree = new Trees::BinaryTree (this->_mesh, 200, build_type);
140 #endif
141  }
142 
143  STOP_LOG("init(no master)", "PointLocatorTree");
144  }
145 
146  else
147 
148  {
149  // We are _not_ the master. Let our Tree point to
150  // the master's tree. But for this we first transform
151  // the master in a state for which we are friends.
152  // And make sure the master @e has a tree!
153  const PointLocatorTree* my_master =
154  libmesh_cast_ptr<const PointLocatorTree*>(this->_master);
155 
156  if (my_master->initialized())
157  this->_tree = my_master->_tree;
158  else
159  {
160  libMesh::err << "ERROR: Initialize master first, then servants!"
161  << std::endl;
162  libmesh_error();
163  }
164  }
165 
166 
167  // Not all PointLocators may own a tree, but all of them
168  // use their own element pointer. Let the element pointer
169  // be unique for every interpolator.
170  // Suppose the interpolators are used concurrently
171  // at different locations in the mesh, then it makes quite
172  // sense to have unique start elements.
173  this->_element = NULL;
174  }
175 
176 
177  // ready for take-off
178  this->_initialized = true;
179 }
180 
181 
182 
183 
184 
185 const Elem* PointLocatorTree::operator() (const Point& p) const
186 {
187  libmesh_assert (this->_initialized);
188 
189  START_LOG("operator()", "PointLocatorTree");
190 
191  // First check the element from last time before asking the tree
192  if (this->_element==NULL || !(this->_element->contains_point(p)))
193  {
194  // ask the tree
195  this->_element = this->_tree->find_element (p);
196 
197  if (this->_element == NULL)
198  {
199  /* No element seems to contain this point. If out-of-mesh
200  mode is enabled, just return NULL. If not, however, we
201  have to perform a linear search before we call \p
202  libmesh_error() since in the case of curved elements, the
203  bounding box computed in \p TreeNode::insert(const
204  Elem*) might be slightly inaccurate. */
205  if(!_out_of_mesh_mode)
206  {
207  START_LOG("linear search", "PointLocatorTree");
208  MeshBase::const_element_iterator pos = this->_mesh.active_elements_begin();
209  const MeshBase::const_element_iterator end_pos = this->_mesh.active_elements_end();
210 
211  for ( ; pos != end_pos; ++pos)
212  if ((*pos)->contains_point(p))
213  {
214  STOP_LOG("linear search", "PointLocatorTree");
215  STOP_LOG("operator()", "PointLocatorTree");
216  return this->_element = (*pos);
217  }
218 
219 /*
220  if (this->_element == NULL)
221  {
222  libMesh::err << std::endl
223  << " ******** Serious Problem. Could not find an Element "
224  << "in the Mesh"
225  << std:: endl
226  << " ******** that contains the Point "
227  << p;
228  libmesh_error();
229  }
230 */
231  STOP_LOG("linear search", "PointLocatorTree");
232  }
233  }
234  }
235 
236  // If we found an element, it should be active
237  libmesh_assert (!this->_element || this->_element->active());
238 
239  STOP_LOG("operator()", "PointLocatorTree");
240 
241  // return the element
242  return this->_element;
243 }
244 
245 void PointLocatorTree::enable_out_of_mesh_mode (void)
246 {
247  /* Out-of-mesh mode is currently only supported if all of the
248  elements have affine mappings. The reason is that for quadratic
249  mappings, it is not easy to construct a relyable bounding box of
250  the element, and thus, the fallback linear search in \p
251  operator() is required. Hence, out-of-mesh mode would be
252  extremely slow. */
253  if(!_out_of_mesh_mode)
254  {
255 #ifdef DEBUG
256  MeshBase::const_element_iterator pos = this->_mesh.active_elements_begin();
257  const MeshBase::const_element_iterator end_pos = this->_mesh.active_elements_end();
258  for ( ; pos != end_pos; ++pos)
259  if (!(*pos)->has_affine_map())
260  {
261  libMesh::err << "ERROR: Out-of-mesh mode is currently only supported if all elements have affine mappings." << std::endl;
262  libmesh_error();
263  }
264 #endif
265 
266  _out_of_mesh_mode = true;
267  }
268 }
269 
270 void PointLocatorTree::disable_out_of_mesh_mode (void)
271 {
272  _out_of_mesh_mode = false;
273 }
274 
275 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo