SCIRun  5.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QuadSurfMesh.h
Go to the documentation of this file.
1 /*
2  For more information, please see: http://software.sci.utah.edu
3 
4  The MIT License
5 
6  Copyright (c) 2009 Scientific Computing and Imaging Institute,
7  University of Utah.
8 
9 
10  Permission is hereby granted, free of charge, to any person obtaining a
11  copy of this software and associated documentation files (the "Software"),
12  to deal in the Software without restriction, including without limitation
13  the rights to use, copy, modify, merge, publish, distribute, sublicense,
14  and/or sell copies of the Software, and to permit persons to whom the
15  Software is furnished to do so, subject to the following conditions:
16 
17  The above copyright notice and this permission notice shall be included
18  in all copies or substantial portions of the Software.
19 
20  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26  DEALINGS IN THE SOFTWARE.
27 */
28 
29 #ifndef CORE_DATATYPES_QUADSURFMESH_H
30 #define CORE_DATATYPES_QUADSURFMESH_H 1
31 
32 /// Include what kind of support we want to have
33 /// Need to fix this and couple it sci-defs
34 #include <Core/Datatypes/Legacy/Field/MeshSupport.h>
35 
36 #include <Core/Containers/StackVector.h>
37 
38 #include <Core/GeometryPrimitives/SearchGridT.h>
39 #include <Core/GeometryPrimitives/BBox.h>
40 #include <Core/GeometryPrimitives/CompGeom.h>
41 #include <Core/GeometryPrimitives/Point.h>
42 #include <Core/GeometryPrimitives/Transform.h>
43 #include <Core/GeometryPrimitives/Vector.h>
44 
45 #include <Core/Basis/Locate.h>
46 #include <Core/Basis/QuadBilinearLgn.h>
47 #include <Core/Basis/QuadBiquadraticLgn.h>
48 #include <Core/Basis/QuadBicubicHmt.h>
49 
50 #include <Core/Datatypes/Legacy/Field/Mesh.h>
51 #include <Core/Datatypes/Legacy/Field/FieldIterator.h>
52 #include <Core/Datatypes/Legacy/Field/FieldRNG.h>
53 #include <Core/Datatypes/Legacy/Field/VMesh.h>
54 #include <Core/Datatypes/Mesh/VirtualMeshFacade.h>
55 
56 #include <Core/Thread/Mutex.h>
57 #include <Core/Thread/ConditionVariable.h>
58 #include <boost/thread.hpp>
59 #include <Core/Persistent/PersistentSTL.h>
60 
61 /// Needed for some specialized functions
62 #include <set>
63 
64 #include <Core/Datatypes/Legacy/Field/share.h>
65 
66 namespace SCIRun {
67 
68 /////////////////////////////////////////////////////
69 // Declarations for virtual interface
70 
71 template <class Basis> class QuadSurfMesh;
72 
73 /// make sure any other mesh other than the preinstantiate ones
74 /// returns no virtual interface. Altering this behavior will allow
75 /// for dynamically compiling the interface if needed.
76 template<class MESH>
77 VMesh* CreateVQuadSurfMesh(MESH* mesh) { return (0); }
78 
79 #if (SCIRUN_QUADSURF_SUPPORT > 0)
80 /// These declarations are needed for a combined dynamic compilation as
81 /// as well as virtual functions solution.
82 /// Declare that these can be found in a library that is already
83 /// precompiled. So dynamic compilation will not instantiate them again.
84 SCISHARE VMesh* CreateVQuadSurfMesh(QuadSurfMesh<Core::Basis::QuadBilinearLgn<Core::Geometry::Point> >* mesh);
85 #if (SCIRUN_QUADRATIC_SUPPORT > 0)
86 SCISHARE VMesh* CreateVQuadSurfMesh(QuadSurfMesh<Core::Basis::QuadBiquadraticLgn<Core::Geometry::Point> >* mesh);
87 #endif
88 #if (SCIRUN_CUBIC_SUPPORT > 0)
89 SCISHARE VMesh* CreateVQuadSurfMesh(QuadSurfMesh<Core::Basis::QuadBicubicHmt<Core::Geometry::Point> >* mesh);
90 #endif
91 
92 #endif
93 /////////////////////////////////////////////////////
94 
95 
96 /////////////////////////////////////////////////////
97 // Declarations for QuadSurfMesh class
98 
99 
100 template <class Basis>
101 class QuadSurfMesh : public Mesh
102 {
103 
104 /// Make sure the virtual interface has access
105 template<class MESH> friend class VQuadSurfMesh;
106 template<class MESH> friend class VMeshShared;
107 template<class MESH> friend class VUnstructuredMesh;
108 
109 public:
110  // Types that change depending on 32 or 64bits
111  typedef SCIRun::index_type under_type;
112  typedef SCIRun::index_type index_type;
113  typedef SCIRun::size_type size_type;
114  typedef SCIRun::mask_type mask_type;
115 
116  typedef boost::shared_ptr<QuadSurfMesh<Basis> > handle_type;
117  typedef Basis basis_type;
118 
119  /// Index and Iterator types required for Mesh Concept.
120  struct Node {
121  typedef NodeIndex<under_type> index_type;
122  typedef NodeIterator<under_type> iterator;
123  typedef NodeIndex<under_type> size_type;
124  typedef StackVector<index_type, 4> array_type;
125  };
126 
127  struct Edge {
128  typedef EdgeIndex<under_type> index_type;
129  typedef EdgeIterator<under_type> iterator;
130  typedef EdgeIndex<under_type> size_type;
131  typedef std::vector<index_type> array_type;
132  };
133 
134  struct Face {
135  typedef FaceIndex<under_type> index_type;
136  typedef FaceIterator<under_type> iterator;
137  typedef FaceIndex<under_type> size_type;
138  typedef std::vector<index_type> array_type;
139  };
140 
141  struct Cell {
142  typedef CellIndex<under_type> index_type;
143  typedef CellIterator<under_type> iterator;
144  typedef CellIndex<under_type> size_type;
145  typedef std::vector<index_type> array_type;
146  };
147 
148  /// Elem refers to the most complex topological object
149  /// DElem refers to object just below Elem in the topological hierarchy
150  typedef Face Elem;
151  typedef Edge DElem;
152 
153  /// Somehow the information of how to interpolate inside an element
154  /// ended up in a separate class, as they need to share information
155  /// this construction was created to transfer data.
156  /// Hopefully in the future this class will disappear again.
157  friend class ElemData;
158 
159  class ElemData
160  {
161  public:
162  typedef typename QuadSurfMesh<Basis>::index_type index_type;
163 
164  ElemData(const QuadSurfMesh<Basis>& msh, const index_type ind) :
165  mesh_(msh),
166  index_(ind)
167  {
168  /// Linear and Constant Basis never use edges_
169  if (basis_type::polynomial_order() > 1)
170  {
171  mesh_.get_edges_from_face(edges_, index_);
172  }
173  }
174 
175  // the following designed to coordinate with ::get_nodes
176  inline
177  index_type node0_index() const
178  {
179  return mesh_.faces_[index_ * 4];
180  }
181  inline
182  index_type node1_index() const
183  {
184  return mesh_.faces_[index_ * 4 + 1];
185  }
186  inline
187  index_type node2_index() const
188  {
189  return mesh_.faces_[index_ * 4 + 2];
190  }
191  inline
192  index_type node3_index() const
193  {
194  return mesh_.faces_[index_ * 4 + 3];
195  }
196 
197  /// changed the indexing as we now use unique indices
198  inline
199  index_type edge0_index() const
200  {
201  return edges_[0];
202  }
203  inline
204  index_type edge1_index() const
205  {
206  return edges_[1];
207  }
208  inline
209  index_type edge2_index() const
210  {
211  return edges_[2];
212  }
213  inline
214  index_type edge3_index() const
215  {
216  return edges_[3];
217  }
218 
219  inline
220  const Core::Geometry::Point &node0() const
221  {
222  return mesh_.points_[node0_index()];
223  }
224  inline
225  const Core::Geometry::Point &node1() const
226  {
227  return mesh_.points_[node1_index()];
228  }
229  inline
230  const Core::Geometry::Point &node2() const
231  {
232  return mesh_.points_[node2_index()];
233  }
234  inline
235  const Core::Geometry::Point &node3() const
236  {
237  return mesh_.points_[node3_index()];
238  }
239 
240  private:
241  /// reference to the mesh
242  const QuadSurfMesh<Basis> &mesh_;
243  /// copy of element index
244  const index_type index_;
245  /// need edges for quadratic meshes
246  typename Edge::array_type edges_;
247  };
248 
249 
250  friend class Synchronize;
251 
252  class Synchronize //: public Runnable
253  {
254  public:
255  Synchronize(QuadSurfMesh<Basis>* mesh, mask_type sync) :
256  mesh_(mesh), sync_(sync)
257  {
258  }
259 
260  void operator()()
261  {
262  run();
263  }
264 
265  void run()
266  {
267  mesh_->synchronize_lock_.lock();
268  // block out all the tables that are already synchronized
269  sync_ &= ~(mesh_->synchronized_);
270  // block out all the tables that are already being computed
271  sync_ &= ~(mesh_->synchronizing_);
272  // Now sync_ contains what this thread will synchronize
273  // Denote what this thread will synchronize
274  mesh_->synchronizing_ |= sync_;
275  // Other threads now know what this tread will be doing
276  mesh_->synchronize_lock_.unlock();
277 
278  // Sync node neighbors
279  if (sync_ & Mesh::NODE_NEIGHBORS_E) mesh_->compute_node_neighbors();
280 
281  if (sync_ & Mesh::EDGES_E) mesh_->compute_edges();
282  if (sync_ & Mesh::NORMALS_E) mesh_->compute_normals();
283 
284  if (sync_ & Mesh::BOUNDING_BOX_E) mesh_->compute_bounding_box();
285 
286  // These depend on the bounding box being synchronized
287  if (sync_ & (Mesh::NODE_LOCATE_E|Mesh::ELEM_LOCATE_E))
288  {
289  {
290  Core::Thread::UniqueLock lock(mesh_->synchronize_lock_.get());
291  while(!(mesh_->synchronized_ & Mesh::BOUNDING_BOX_E))
292  {
293  mesh_->synchronize_cond_.wait(lock);
294  }
295  }
296  if (sync_ & Mesh::NODE_LOCATE_E)
297  {
298  mesh_->compute_node_grid();
299  }
300  if (sync_ & Mesh::ELEM_LOCATE_E)
301  {
302  mesh_->compute_elem_grid();
303  }
304  }
305 
306  mesh_->synchronize_lock_.lock();
307  // Mark the ones that were just synchronized
308  mesh_->synchronized_ |= sync_;
309  // Unmark the the ones that were done
310  mesh_->synchronizing_ &= ~(sync_);
311  mesh_->synchronize_lock_.unlock();
312  /// Tell other threads we are done
313  mesh_->synchronize_cond_.conditionBroadcast();
314  }
315 
316  private:
317  QuadSurfMesh<Basis>* mesh_;
318  mask_type sync_;
319  };
320 
321 
322  //////////////////////////////////////////////////////////////////
323 
324  /// Construct a new mesh
325  QuadSurfMesh();
326 
327  /// Copy a mesh, needed for detaching the mesh from a field
328  QuadSurfMesh(const QuadSurfMesh& copy);
329 
330  /// Clone function for detaching the mesh and automatically generating
331  /// a new version if needed.
332  virtual QuadSurfMesh *clone() const { return new QuadSurfMesh(*this); }
333 
334  /// Destructor
335  virtual ~QuadSurfMesh();
336 
337  MeshFacadeHandle getFacade() const
338  {
339  return boost::make_shared<Core::Datatypes::VirtualMeshFacade<VMesh>>(vmesh_);
340  }
341 
342  /// Access point to virtual interface
343  virtual VMesh* vmesh() { return vmesh_.get(); }
344 
345  /// This one should go at some point, should be reroute through the
346  /// virtual interface
347  virtual int basis_order() { return (basis_.polynomial_order()); }
348 
349  /// Topological dimension
350  virtual int dimensionality() const { return 2; }
351 
352  /// What kind of mesh is this
353  /// structured = no connectivity data
354  /// regular = no node location data
355  virtual int topology_geometry() const
356  { return (Mesh::UNSTRUCTURED | Mesh::IRREGULAR); }
357 
358  /// Get the bounding box of the field
359  virtual Core::Geometry::BBox get_bounding_box() const;
360 
361  /// Return the transformation that takes a 0-1 space bounding box
362  /// to the current bounding box of this mesh.
363  virtual void get_canonical_transform(Core::Geometry::Transform &t) const;
364 
365  /// Core::Geometry::Transform a field (transform all nodes using this transformation matrix)
366  virtual void transform(const Core::Geometry::Transform &t);
367 
368  /// Check whether mesh can be altered by adding nodes or elements
369  virtual bool is_editable() const { return (true); }
370 
371  /// Has this mesh normals.
372  // Note: normals point inward.
373  /// @todo: this is inconsistent with TriSurfMesh - should both surfaces
374  // have consistent normal direction?
375  virtual bool has_normals() const { return (true); }
376 
377  /// Compute tables for doing topology, these need to be synchronized
378  ///before doing a lot of operations.
379  virtual bool synchronize(mask_type mask);
380  virtual bool unsynchronize(mask_type mask);
381  bool clear_synchronization();
382 
383  /// Get the basis class.
384  Basis& get_basis() { return basis_; }
385 
386  /// begin/end iterators
387  void begin(typename Node::iterator &) const;
388  void begin(typename Edge::iterator &) const;
389  void begin(typename Face::iterator &) const;
390  void begin(typename Cell::iterator &) const;
391 
392  void end(typename Node::iterator &) const;
393  void end(typename Edge::iterator &) const;
394  void end(typename Face::iterator &) const;
395  void end(typename Cell::iterator &) const;
396 
397  void size(typename Node::size_type &) const;
398  void size(typename Edge::size_type &) const;
399  void size(typename Face::size_type &) const;
400  void size(typename Cell::size_type &) const;
401 
402  /// These are here to convert indices to unsigned int
403  /// counters. Some how the decision was made to use multi
404  /// dimensional indices in some fields, these functions
405  /// should deal with different pointer types.
406  /// Use the virtual interface to avoid all this non sense.
407  void to_index(typename Node::index_type &index, index_type i) const
408  { index = i; }
409  void to_index(typename Edge::index_type &index, index_type i) const
410  { index = i; }
411  void to_index(typename Face::index_type &index, index_type i) const
412  { index = i; }
413  void to_index(typename Cell::index_type&, index_type) const
414  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
415 
416  /// Get the child elements of the given index.
417  void get_nodes(typename Node::array_type &array,
418  typename Node::index_type idx) const
419  { array.resize(1); array[0]= idx; }
420  void get_nodes(typename Node::array_type &array,
421  typename Edge::index_type idx) const
422  { get_nodes_from_edge(array,idx); }
423  void get_nodes(typename Node::array_type &array,
424  typename Face::index_type idx) const
425  { get_nodes_from_face(array,idx); }
426  void get_nodes(typename Node::array_type&,
427  typename Cell::index_type) const
428  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
429 
430  void get_edges(typename Edge::array_type &array,
431  typename Node::index_type idx)
432  { get_edges_from_node(array,idx); }
433  void get_edges(typename Edge::array_type &array,
434  typename Edge::index_type idx) const
435  { array.resize(1); array[0]= idx; }
436  void get_edges(typename Edge::array_type &array,
437  typename Face::index_type idx) const
438  { get_edges_from_face(array,idx); }
439  void get_edges(typename Edge::array_type&,
440  typename Cell::index_type) const
441  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
442 
443  void get_faces(typename Face::array_type &array,
444  typename Node::index_type idx) const
445  { get_faces_from_node(array,idx); }
446  void get_faces(typename Face::array_type &array,
447  typename Edge::index_type idx) const
448  { get_faces_from_edge(array,idx); }
449  void get_faces(typename Face::array_type &array,
450  typename Face::index_type idx) const
451  { array.resize(1); array[0]= idx; }
452  void get_faces(typename Face::array_type&,
453  typename Cell::index_type) const
454  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
455 
456  void get_cells(typename Cell::array_type&,
457  typename Node::index_type) const
458  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
459  void get_cells(typename Cell::array_type&,
460  typename Edge::index_type) const
461  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
462  void get_cells(typename Cell::array_type&,
463  typename Face::index_type) const
464  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
465  void get_cells(typename Cell::array_type&,
466  typename Cell::index_type) const
467  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
468 
469  void get_elems(typename Elem::array_type &array,
470  typename Node::index_type idx) const
471  { get_faces_from_node(array,idx); }
472  void get_elems(typename Elem::array_type &array,
473  typename Edge::index_type idx) const
474  { get_faces_from_edge(array,idx); }
475  void get_elems(typename Elem::array_type &array,
476  typename Face::index_type idx) const
477  { array.resize(1); array[0]= idx; }
478  void get_elems(typename Face::array_type&,
479  typename Cell::index_type) const
480  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
481 
482  void get_delems(typename DElem::array_type &array,
483  typename Node::index_type idx) const
484  { get_edges_from_node(array,idx); }
485  void get_delems(typename DElem::array_type &array,
486  typename Edge::index_type idx) const
487  { array.resize(1); array[0]= idx; }
488  void get_delems(typename DElem::array_type &array,
489  typename Face::index_type idx) const
490  { get_edges_from_face(array,idx); }
491  void get_delems(typename DElem::array_type&,
492  typename Cell::index_type) const
493  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
494 
495  /// Generate the list of points that make up a sufficiently accurate
496  /// piecewise linear approximation of an edge.
497  template<class VECTOR, class INDEX>
498  void pwl_approx_edge(std::vector<VECTOR > &coords,
499  INDEX ci,
500  unsigned int which_edge,
501  unsigned int div_per_unit) const
502  {
503  basis_.approx_edge(which_edge, div_per_unit, coords);
504  }
505 
506  /// Generate the list of points that make up a sufficiently accurate
507  /// piecewise linear approximation of an face.
508  template<class VECTOR, class INDEX>
509  void pwl_approx_face(std::vector<std::vector<VECTOR > > &coords,
510  INDEX ci,
511  unsigned int which_face,
512  unsigned int div_per_unit) const
513  {
514  basis_.approx_face(which_face, div_per_unit, coords);
515  }
516 
517  /// get the center point (in object space) of an element
518  void get_center(Core::Geometry::Point &result, typename Node::index_type idx) const
519  { get_node_center(result, idx); }
520  void get_center(Core::Geometry::Point &result, typename Edge::index_type idx) const
521  { get_edge_center(result, idx); }
522  void get_center(Core::Geometry::Point &result, typename Face::index_type idx) const
523  { get_face_center(result, idx); }
524  void get_center(Core::Geometry::Point&, typename Cell::index_type) const
525  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
526 
527  /// Get the size of an element (length, area, volume)
528  double get_size(typename Node::index_type /*idx*/) const
529  { return 0.0; }
530  double get_size(typename Edge::index_type idx) const
531  {
532  typename Node::array_type arr;
533  get_nodes_from_edge(arr, idx);
534  return (point(arr[0]) - point(arr[1])).length();
535  }
536  double get_size(typename Face::index_type idx) const
537  {
538  typename Node::array_type ra;
539  get_nodes_from_face(ra,idx);
540  const Core::Geometry::Point &p0 = point(ra[0]);
541  const Core::Geometry::Point &p1 = point(ra[1]);
542  const Core::Geometry::Point &p2 = point(ra[2]);
543  const Core::Geometry::Point &p3 = point(ra[3]);
544  return ((Cross(p0-p1,p2-p1)).length()+(Cross(p2-p3,p0-p3)).length()+
545  (Cross(p3-p0,p1-p0)).length()+(Cross(p1-p2,p3-p2)).length())*0.25;
546  }
547  double get_size(typename Cell::index_type) const
548  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
549 
550  /// More specific names for get_size
551  double get_length(typename Edge::index_type idx) const
552  { return get_size(idx); }
553  double get_area(typename Face::index_type idx) const
554  { return get_size(idx); }
555  double get_volume(typename Cell::index_type) const
556  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
557 
558 
559 
560  /// Get neighbors of an element or a node
561 
562  /// THIS ONE IS FLAWED AS IN 3D SPACE MULTIPLE EDGES CAN CONNECTED THROUGH
563  /// ONE EDGE
564  bool get_neighbor(typename Elem::index_type &neighbor,
565  typename Elem::index_type elem,
566  typename DElem::index_type delem) const
567  { return(get_elem_neighbor(neighbor,elem,delem)); }
568 
569  /// These are more general implementations
570  void get_neighbors(std::vector<typename Node::index_type> &array,
571  typename Node::index_type node)
572  { get_node_neighbors(array,node); }
573  bool get_neighbors(std::vector<typename Elem::index_type> &array,
574  typename Elem::index_type elem,
575  typename DElem::index_type delem) const
576  { return(get_elem_neighbors(array,elem,delem)); }
577  void get_neighbors(typename Elem::array_type &array,
578  typename Elem::index_type elem) const
579  { get_elem_neighbors(array,elem); }
580 
581 
582  /// Locate a point in a mesh, find which is the closest node
583  bool locate(typename Node::index_type &node, const Core::Geometry::Point &p) const
584  { return (locate_node(node,p)); }
585  bool locate(typename Edge::index_type &edge, const Core::Geometry::Point &p) const
586  { return (locate_edge(edge,p)); }
587  bool locate(typename Face::index_type &face, const Core::Geometry::Point &p) const
588  { return (locate_elem(face,p)); }
589  bool locate(typename Cell::index_type&, const Core::Geometry::Point &) const
590  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
591 
592  bool locate(typename Elem::index_type &elem,
593  std::vector<double>& coords,
594  const Core::Geometry::Point &p) const
595  { return (locate_elem(elem,coords,p)); }
596 
597 
598  /// These should become obsolete soon, they do not follow the concept
599  /// of the basis functions....
600  int get_weights(const Core::Geometry::Point &p, typename Node::array_type &l, double *w);
601  int get_weights(const Core::Geometry::Point & , typename Edge::array_type & , double * )
602  { ASSERTFAIL("QuadSurfMesh: get_weights(edges) is not supported."); }
603  int get_weights(const Core::Geometry::Point &p, typename Face::array_type &l, double *w);
604  int get_weights(const Core::Geometry::Point & , typename Cell::array_type & , double * )
605  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
606 
607  /// Access the nodes of the mesh
608  void get_point(Core::Geometry::Point &p, typename Node::index_type i) const
609  { p = points_[i]; }
610  void set_point(const Core::Geometry::Point &p, typename Node::index_type i)
611  { points_[i] = p; }
612 
613  void get_random_point(Core::Geometry::Point &, typename Elem::index_type, FieldRNG &rng) const;
614 
615  /// Normals for visualizations
616  void get_normal(Core::Geometry::Vector &n, typename Node::index_type i) const
617  {
618  ASSERTMSG(synchronized_ & Mesh::NORMALS_E,
619  "Must call synchronize NORMALS_E on QuadSurfMesh first");
620  n = normals_[i];
621  }
622 
623  /// Get the normals at the outside of the element
624  template<class VECTOR, class INDEX1, class INDEX2>
625  void get_normal(Core::Geometry::Vector &result, VECTOR &coords,
626  INDEX1 eidx,INDEX2 /*fidx*/)
627  {
628 /*
629  * THE ISSUE IS WHAT THE NORMAL IS AT THE NODE, THIS FUNCTION TAKES AN ELEMENT
630  * INDEX AS WEL, HENCE IT SHOULD RETURN THE NORMAL OF THAT FACE, NOT AN AVERAGED
631  * ONE OVER NEIGHBORING FACES
632 
633  if (basis_.polynomial_order() < 2)
634  {
635  ASSERTMSG(synchronized_ & Mesh::NORMALS_E, "Must call synchronize NORMALS_E on TriSurfMesh first");
636 
637  typename Node::array_type arr(3);
638  get_nodes_from_face(arr, eidx);
639 
640  const double c0_0 = fabs(coords[0]);
641  const double c1_0 = fabs(coords[1]);
642  const double c0_1 = fabs(coords[0] - 1.0L);
643  const double c1_1 = fabs(coords[1] - 1.0L);
644 
645  if (c0_0 < 1e-7 && c1_0 < 1e-7) {
646  // arr[0]
647  result = normals_[arr[0]];
648  return;
649  } else if (c0_1 < 1e-7 && c1_0 < 1e-7) {
650  // arr[1]
651  result = normals_[arr[1]];
652  return;
653  } else if (c0_1 < 1e-7 && c1_1 < 1e-7) {
654  // arr[2]
655  result = normals_[arr[2]];
656  return;
657  } else if (c0_0 < 1e-7 && c1_1 < 1e-7) {
658  // arr[3]
659  result = normals_[arr[3]];
660  return;
661  }
662  }
663 */
664 
665  ElemData ed(*this, eidx);
666  std::vector<Core::Geometry::Point> Jv;
667  basis_.derivate(coords, ed, Jv);
668  result = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1]));
669  result.normalize();
670  }
671 
672  /// Add a new node to the mesh
673  typename Node::index_type add_point(const Core::Geometry::Point &p);
674  typename Node::index_type add_node(const Core::Geometry::Point &p)
675  { return(add_point(p)); }
676 
677  /// Add a new element to the mesh
678  template <class ARRAY>
679  typename Elem::index_type add_elem(ARRAY a)
680  {
681  ASSERTMSG(a.size() == 4, "Tried to add non-quad element.");
682  ASSERTMSG(order_face_nodes(a[0],a[1],a[2],a[3]), "add_elem: element that is being added is invalid");
683 
684  faces_.push_back(static_cast<typename Node::index_type>(a[0]));
685  faces_.push_back(static_cast<typename Node::index_type>(a[1]));
686  faces_.push_back(static_cast<typename Node::index_type>(a[2]));
687  faces_.push_back(static_cast<typename Node::index_type>(a[3]));
688  return static_cast<typename Elem::index_type>((static_cast<index_type>(faces_.size()) - 1) >> 2);
689  }
690 
691  void node_reserve(size_type s) { points_.reserve(static_cast<size_t>(s)); }
692  void elem_reserve(size_type s) { faces_.reserve(static_cast<size_t>(s*4)); }
693  void resize_nodes(size_type s) { points_.resize(static_cast<size_t>(s)); }
694  void resize_elems(size_type s) { faces_.resize(static_cast<size_t>(s*4)); }
695 
696 
697  /// Get the local coordinates for a certain point within an element
698  /// This function uses a couple of newton iterations to find the local
699  /// coordinate of a point
700  template<class VECTOR, class INDEX>
701  bool get_coords(VECTOR &coords, const Core::Geometry::Point &p, INDEX idx) const
702  {
703  ElemData ed(*this, idx);
704  return basis_.get_coords(coords, p, ed);
705  }
706 
707  /// Find the location in the global coordinate system for a local coordinate
708  /// This function is the opposite of get_coords.
709  template<class VECTOR, class INDEX>
710  void interpolate(Core::Geometry::Point &pt, const VECTOR &coords, INDEX idx) const
711  {
712  ElemData ed(*this, idx);
713  pt = basis_.interpolate(coords, ed);
714  }
715 
716  /// Interpolate the derivate of the function, This infact will return the
717  /// jacobian of the local to global coordinate transformation. This function
718  /// is mainly intended for the non linear elements
719  template<class VECTOR1, class INDEX, class VECTOR2>
720  void derivate(const VECTOR1 &coords, INDEX idx, VECTOR2 &J) const
721  {
722  ElemData ed(*this, idx);
723  basis_.derivate(coords, ed, J);
724  }
725 
726  /// Get the determinant of the jacobian, which is the local volume of an element
727  /// and is intended to help with the integration of functions over an element.
728  template<class VECTOR, class INDEX>
729  double det_jacobian(const VECTOR& coords, INDEX idx) const
730  {
731  double J[9];
732  jacobian(coords,idx,J);
733  return (DetMatrix3x3(J));
734  }
735 
736  /// Get the jacobian of the transformation. In case one wants the non inverted
737  /// version of this matrix. This is currently here for completeness of the
738  /// interface
739  template<class VECTOR, class INDEX>
740  void jacobian(const VECTOR& coords, INDEX idx, double* J) const
741  {
742  StackVector<Core::Geometry::Point,2> Jv;
743  ElemData ed(*this,idx);
744  basis_.derivate(coords,ed,Jv);
745  Core::Geometry::Vector Jv2 = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1]));
746  Jv2.normalize();
747  J[0] = Jv[0].x();
748  J[1] = Jv[0].y();
749  J[2] = Jv[0].z();
750  J[3] = Jv[1].x();
751  J[4] = Jv[1].y();
752  J[5] = Jv[1].z();
753  J[6] = Jv2.x();
754  J[7] = Jv2.y();
755  J[8] = Jv2.z();
756  }
757 
758  /// Get the inverse jacobian of the transformation. This one is needed to
759  /// translate local gradients into global gradients. Hence it is crucial for
760  /// calculating gradients of fields, or constructing finite elements.
761  template<class VECTOR, class INDEX>
762  double inverse_jacobian(const VECTOR& coords, INDEX idx, double* Ji) const
763  {
764  StackVector<Core::Geometry::Point,2> Jv;
765  ElemData ed(*this,idx);
766  basis_.derivate(coords,ed,Jv);
767  double J[9];
768  Core::Geometry::Vector Jv2 = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1]));
769  Jv2.normalize();
770  J[0] = Jv[0].x();
771  J[1] = Jv[0].y();
772  J[2] = Jv[0].z();
773  J[3] = Jv[1].x();
774  J[4] = Jv[1].y();
775  J[5] = Jv[1].z();
776  J[6] = Jv2.x();
777  J[7] = Jv2.y();
778  J[8] = Jv2.z();
779 
780  return (InverseMatrix3x3(J,Ji));
781  }
782 
783 
784  template<class INDEX>
785  double scaled_jacobian_metric(INDEX idx) const
786  {
787  StackVector<Core::Geometry::Point,3> Jv;
788  ElemData ed(*this,idx);
789 
790  double temp;
791 
792  basis_.derivate(basis_.unit_center,ed,Jv);
793  Jv.resize(3);
794  Core::Geometry::Vector v = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1]));
795  v.normalize();
796  Jv[2] = v.asPoint();
797  double min_jacobian = ScaledDetMatrix3P(Jv);
798 
799  size_t num_vertices = basis_.number_of_vertices();
800  for (size_t j=0;j < num_vertices;j++)
801  {
802  basis_.derivate(basis_.unit_vertices[j],ed,Jv);
803  Jv.resize(3);
804  v = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1])); v.normalize();
805  Jv[2] = v.asPoint();
806  temp = ScaledDetMatrix3P(Jv);
807  if(temp < min_jacobian) min_jacobian = temp;
808  }
809 
810  return (min_jacobian);
811  }
812 
813 
814  template<class INDEX>
815  double jacobian_metric(INDEX idx) const
816  {
817  StackVector<Core::Geometry::Point,3> Jv;
818  ElemData ed(*this,idx);
819 
820  double temp;
821 
822  basis_.derivate(basis_.unit_center,ed,Jv);
823  Jv.resize(3);
824  Core::Geometry::Vector v = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1])); v.normalize();
825  Jv[2] = v.asPoint();
826  double min_jacobian = DetMatrix3P(Jv);
827 
828  size_t num_vertices = basis_.number_of_vertices();
829  for (size_t j=0;j < num_vertices;j++)
830  {
831  basis_.derivate(basis_.unit_vertices[j],ed,Jv);
832  Jv.resize(3);
833  v = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1])); v.normalize();
834  Jv[2] = v.asPoint();
835  temp = DetMatrix3P(Jv);
836  if(temp < min_jacobian) min_jacobian = temp;
837  }
838 
839  return (min_jacobian);
840  }
841 
842  template<class INDEX>
843  double inscribed_circumscribed_radius_metric(INDEX idx) const
844  {
845  return (0.0);
846  }
847 
848  template <class INDEX>
849  bool find_closest_node(double& pdist, Core::Geometry::Point &result,
850  INDEX &node, const Core::Geometry::Point &p) const
851  {
852  return(find_closest_node(pdist,result,node,p,-1.0));
853  }
854 
855  template <class INDEX>
856  bool find_closest_node(double& pdist, Core::Geometry::Point &result,
857  INDEX &node, const Core::Geometry::Point &p,double maxdist) const
858  {
859  if (maxdist < 0.0) maxdist = DBL_MAX; else maxdist = maxdist*maxdist;
860  typename Node::size_type sz; size(sz);
861 
862  /// If there are no nodes we cannot find the closest one
863  if (sz == 0) return (false);
864 
865  if (node >= 0 && node < sz)
866  {
867  Core::Geometry::Point point = points_[node];
868  double dist = (point-p).length2();
869 
870  if ( dist < epsilon2_ )
871  {
872  result = point;
873  pdist = sqrt(dist);
874  return (true);
875  }
876  }
877 
878  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
879  "QuadSurfMesh::find_closest_node requires synchronize(NODE_LOCATE_E).")
880 
881  // get grid sizes
882  const size_type ni = node_grid_->get_ni()-1;
883  const size_type nj = node_grid_->get_nj()-1;
884  const size_type nk = node_grid_->get_nk()-1;
885 
886  // Convert to grid coordinates.
887  index_type bi, bj, bk, ei, ej, ek;
888  node_grid_->unsafe_locate(bi, bj, bk, p);
889 
890  // Clamp to closest point on the grid.
891  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
892  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
893  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
894 
895  ei = bi; ej = bj; ek = bk;
896 
897  double dmin = maxdist;
898  bool found = true;
899  bool found_one = false;
900 
901  do
902  {
903  found = true;
904  /// This looks incorrect - but it is correct
905  /// We need to do a full shell without any elements that are closer
906  /// to make sure there no closer elements in neighboring searchgrid cells
907 
908  for (index_type i = bi; i <= ei; i++)
909  {
910  if (i < 0 || i > ni) continue;
911  for (index_type j = bj; j <= ej; j++)
912  {
913  if (j < 0 || j > nj) continue;
914  for (index_type k = bk; k <= ek; k++)
915  {
916  if (k < 0 || k > nk) continue;
917  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
918  {
919  if (node_grid_->min_distance_squared(p, i, j, k) < dmin)
920  {
921  typename SearchGridT<index_type>::iterator it, eit;
922  found = false;
923  node_grid_->lookup_ijk(it,eit, i, j, k);
924 
925  while (it != eit)
926  {
927  const Core::Geometry::Point point = points_[*it];
928  const double dist = (p-point).length2();
929 
930  if (dist < dmin)
931  {
932  found_one = true;
933  result = point;
934  node = INDEX(*it);
935  dmin = dist;
936 
937  /// If we are closer than eps^2 we found a node close enough
938  if (dmin < epsilon2_)
939  {
940  pdist = sqrt(dmin);
941  return (true);
942  }
943  }
944  ++it;
945  }
946  }
947  }
948  }
949  }
950  }
951  bi--;ei++;
952  bj--;ej++;
953  bk--;ek++;
954  }
955  while (!found) ;
956 
957  if (!found_one) return (false);
958 
959  pdist = sqrt(dmin);
960  return (true);
961  }
962 
963  template <class ARRAY>
964  bool find_closest_nodes(ARRAY &nodes, const Core::Geometry::Point &p, double maxdist) const
965  {
966  nodes.clear();
967 
968  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
969  "QuadSurfMesh::find_closest_node requires synchronize(NODE_LOCATE_E).")
970 
971  // get grid sizes
972  const size_type ni = node_grid_->get_ni()-1;
973  const size_type nj = node_grid_->get_nj()-1;
974  const size_type nk = node_grid_->get_nk()-1;
975 
976  // Convert to grid coordinates.
977  index_type bi, bj, bk, ei, ej, ek;
978 
979  Core::Geometry::Point max = p+Core::Geometry::Vector(maxdist,maxdist,maxdist);
980  Core::Geometry::Point min = p+Core::Geometry::Vector(-maxdist,-maxdist,-maxdist);
981 
982  node_grid_->unsafe_locate(bi, bj, bk, min);
983  node_grid_->unsafe_locate(ei, ej, ek, max);
984 
985  // Clamp to closest point on the grid.
986  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
987  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
988  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
989 
990  if (ei > ni) ei = ni; if (ei < 0) ei = 0;
991  if (ej > nj) ej = nj; if (ej < 0) ej = 0;
992  if (ek > nk) ek = nk; if (ek < 0) ek = 0;
993 
994  double maxdist2 = maxdist*maxdist;
995 
996  for (index_type i = bi; i <= ei; i++)
997  {
998  for (index_type j = bj; j <= ej; j++)
999  {
1000  for (index_type k = bk; k <= ek; k++)
1001  {
1002  if (node_grid_->min_distance_squared(p, i, j, k) < maxdist2)
1003  {
1004  typename SearchGridT<index_type>::iterator it, eit;
1005  node_grid_->lookup_ijk(it, eit, i, j, k);
1006 
1007  while (it != eit)
1008  {
1009  const Core::Geometry::Point point = points_[*it];
1010  const double dist = (p-point).length2();
1011 
1012  if (dist < maxdist2)
1013  {
1014  nodes.push_back(*it);
1015  }
1016  ++it;
1017  }
1018  }
1019  }
1020  }
1021  }
1022 
1023  return(nodes.size() > 0);
1024  }
1025 
1026  template <class ARRAY1, class ARRAY2>
1027  bool find_closest_nodes(ARRAY1 &distances,ARRAY2 &nodes, const Core::Geometry::Point &p, double maxdist) const
1028  {
1029  nodes.clear();
1030  distances.clear();
1031 
1032  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
1033  "QuadSurfMesh::find_closest_node requires synchronize(NODE_LOCATE_E).")
1034 
1035  // get grid sizes
1036  const size_type ni = node_grid_->get_ni()-1;
1037  const size_type nj = node_grid_->get_nj()-1;
1038  const size_type nk = node_grid_->get_nk()-1;
1039 
1040  // Convert to grid coordinates.
1041  index_type bi, bj, bk, ei, ej, ek;
1042 
1043  Core::Geometry::Point max = p+Core::Geometry::Vector(maxdist,maxdist,maxdist);
1044  Core::Geometry::Point min = p+Core::Geometry::Vector(-maxdist,-maxdist,-maxdist);
1045 
1046  node_grid_->unsafe_locate(bi, bj, bk, min);
1047  node_grid_->unsafe_locate(ei, ej, ek, max);
1048 
1049  // Clamp to closest point on the grid.
1050  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
1051  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
1052  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
1053 
1054  if (ei > ni) ei = ni; if (ei < 0) ei = 0;
1055  if (ej > nj) ej = nj; if (ej < 0) ej = 0;
1056  if (ek > nk) ek = nk; if (ek < 0) ek = 0;
1057 
1058  double maxdist2 = maxdist*maxdist;
1059 
1060  for (index_type i = bi; i <= ei; i++)
1061  {
1062  for (index_type j = bj; j <= ej; j++)
1063  {
1064  for (index_type k = bk; k <= ek; k++)
1065  {
1066  if (node_grid_->min_distance_squared(p, i, j, k) < maxdist2)
1067  {
1068  typename SearchGridT<index_type>::iterator it, eit;
1069  node_grid_->lookup_ijk(it, eit, i, j, k);
1070 
1071  while (it != eit)
1072  {
1073  const Core::Geometry::Point point = points_[*it];
1074  const double dist = (p-point).length2();
1075 
1076  if (dist < maxdist2)
1077  {
1078  nodes.push_back(*it);
1079  distances.push_back(dist);
1080  }
1081  ++it;
1082  }
1083  }
1084  }
1085  }
1086  }
1087 
1088  return(nodes.size() > 0);
1089  }
1090 
1091 
1092 
1093 
1094  /// This function will find the closest element and the location on that
1095  /// element that is the closest
1096 
1097  template <class INDEX, class ARRAY>
1098  bool find_closest_elem(double& pdist,
1099  Core::Geometry::Point &result,
1100  ARRAY& coords,
1101  INDEX &elem,
1102  const Core::Geometry::Point &p) const
1103  {
1104  return(find_closest_elem(pdist,result,coords,elem,p,-1.0));
1105  }
1106 
1107  template <class INDEX, class ARRAY>
1108  bool find_closest_elem(double& pdist,
1109  Core::Geometry::Point &result,
1110  ARRAY& coords,
1111  INDEX &elem,
1112  const Core::Geometry::Point &p,
1113  double maxdist) const
1114  {
1115  if (maxdist < 0.0) maxdist = DBL_MAX; else maxdist = maxdist*maxdist;
1116 
1117  typename Elem::size_type sz; size(sz);
1118 
1119  /// If there are no nodes we cannot find the closest one
1120  if (sz == 0) return (false);
1121 
1122  /// Test the one in face that is an initial guess
1123  if (elem >= 0 && elem < sz)
1124  {
1125  const index_type idx = elem * 4;
1126  est_closest_point_on_quad(result, p,
1127  points_[faces_[idx ]],
1128  points_[faces_[idx+1]],
1129  points_[faces_[idx+2]],
1130  points_[faces_[idx+3]]);
1131 
1132  double dist = (p-result).length2();
1133  if ( dist < epsilon2_ )
1134  {
1135  /// As we computed an estimate, we use the Newton's method in the basis functions
1136  /// compute a more refined solution. This function may slow down computation.
1137  /// This piece of code will calculate the coordinates in the local element framework
1138  /// (the newton's method of finding a minimum), then it will project this back
1139  /// THIS CODE SHOULD BE FURTHER OPTIMIZED
1140 
1141  ElemData ed(*this,elem);
1142  basis_.get_coords(coords,result,ed);
1143 
1144  result = basis_.interpolate(coords,ed);
1145  double dmin = (result-p).length2();
1146  pdist = sqrt(dmin);
1147  return (true);
1148  }
1149  }
1150 
1151  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
1152  "QuadSurfMesh::find_closest_elem requires synchronize(ELEM_LOCATE_E).")
1153  // get grid sizes
1154  const size_type ni = elem_grid_->get_ni()-1;
1155  const size_type nj = elem_grid_->get_nj()-1;
1156  const size_type nk = elem_grid_->get_nk()-1;
1157 
1158  // Convert to grid coordinates.
1159  index_type bi, bj, bk, ei, ej, ek;
1160  elem_grid_->unsafe_locate(bi, bj, bk, p);
1161 
1162  // Clamp to closest point on the grid.
1163  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
1164  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
1165  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
1166 
1167  ei = bi; ej = bj; ek = bk;
1168 
1169  double dmin = maxdist;
1170  bool found = true;
1171  bool found_one = false;
1172 
1173  do
1174  {
1175  found = true;
1176  /// This looks incorrect - but it is correct
1177  /// We need to do a full shell without any elements that are closer
1178  /// to make sure there no closer elements in neighboring searchgrid cells
1179 
1180  for (index_type i = bi; i <= ei; i++)
1181  {
1182  if (i < 0 || i > ni) continue;
1183  for (index_type j = bj; j <= ej; j++)
1184  {
1185  if (j < 0 || j > nj) continue;
1186  for (index_type k = bk; k <= ek; k++)
1187  {
1188  if (k < 0 || k > nk) continue;
1189  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
1190  {
1191  if (elem_grid_->min_distance_squared(p, i, j, k) < dmin)
1192  {
1193  found = false;
1194  typename SearchGridT<index_type>::iterator it, eit;
1195  elem_grid_->lookup_ijk(it, eit, i, j, k);
1196 
1197  while (it != eit)
1198  {
1199  Core::Geometry::Point r;
1200  const index_type idx = (*it) * 4;
1201  est_closest_point_on_quad(r, p,
1202  points_[faces_[idx ]],
1203  points_[faces_[idx+1]],
1204  points_[faces_[idx+2]],
1205  points_[faces_[idx+3]]);
1206  const double dist = (p - r).length2();
1207  if (dist < dmin)
1208  {
1209  found_one = true;
1210  result = r;
1211  elem = INDEX(*it);
1212  dmin = dist;
1213  if (dmin < epsilon2_)
1214  {
1215  /// As we computed an estimate, we use the Newton's method in the basis functions
1216  /// compute a more refined solution. This function may slow down computation.
1217  /// This piece of code will calculate the coordinates in the local element framework
1218  /// (the newton's method of finding a minimum), then it will project this back
1219  /// THIS CODE SHOULD BE FURTHER OPTIMIZED
1220 
1221  ElemData ed(*this,elem);
1222  basis_.get_coords(coords,result,ed);
1223 
1224  result = basis_.interpolate(coords,ed);
1225  dmin = (result-p).length2();
1226 
1227  pdist = sqrt(dmin);
1228  return (true);
1229  }
1230  }
1231  ++it;
1232  }
1233  }
1234  }
1235  }
1236  }
1237  }
1238  bi--;ei++;
1239  bj--;ej++;
1240  bk--;ek++;
1241  }
1242  while (!found);
1243 
1244  if (!found_one) return (false);
1245 
1246  /// As we computed an estimate, we use the Newton's method in the basis functions
1247  /// compute a more refined solution. This function may slow down computation.
1248  /// This piece of code will calculate the coordinates in the local element framework
1249  /// (the newton's method of finding a minimum), then it will project this back
1250  /// THIS CODE SHOULD BE FURTHER OPTIMIZED
1251 
1252  ElemData ed(*this,elem);
1253  basis_.get_coords(coords,result,ed);
1254 
1255  result = basis_.interpolate(coords,ed);
1256  dmin = (result-p).length2();
1257 
1258  pdist = sqrt(dmin);
1259 
1260  return (true);
1261  }
1262 
1263  /// This function will find the closest element and the location on that
1264  /// element that is the closest; it takes an additional repelpos argument
1265  /// which defines an exclusion zone that's not considered. A repelnormal
1266  /// can also be specified in order to restrict the faces to those that are
1267  /// pointing in a different direction from the repel point.
1268  template <class INDEX, class ARRAY>
1269  bool find_closest_elem_far_from(double& pdist,
1270  Core::Geometry::Point &result,
1271  ARRAY &coords,
1272  INDEX &face,
1273  const Core::Geometry::Point &p,
1274  const Core::Geometry::Point &repelpos,
1275  double mindistfromrepel,
1276  const Core::Geometry::Vector &repelnormal,
1277  double mindotfromrepel,
1278  double maxdist) const
1279  {
1280  if (maxdist < 0.0) maxdist = DBL_MAX; else maxdist = maxdist*maxdist;
1281  typename Elem::size_type sz; size(sz);
1282 
1283  /// If there are no nodes we cannot find the closest one
1284  if (sz == 0) return (false);
1285 
1286  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
1287  "QuadSurfMesh::find_closest_elem requires synchronize(ELEM_LOCATE_E).")
1288 
1289  // get grid sizes
1290  const size_type ni = elem_grid_->get_ni()-1;
1291  const size_type nj = elem_grid_->get_nj()-1;
1292  const size_type nk = elem_grid_->get_nk()-1;
1293 
1294  // Convert to grid coordinates.
1295  index_type bi, ei, bj, ej, bk, ek;
1296  elem_grid_->unsafe_locate(bi, bj, bk, p);
1297 
1298  // Clamp to closest point on the grid.
1299  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
1300  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
1301  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
1302 
1303  ei = bi; ej = bj; ek = bk;
1304 
1305  double dmin = maxdist;
1306  bool found = true;
1307  bool found_one = false;
1308  double mindistfromrepelsquared = mindistfromrepel * mindistfromrepel;
1309  do
1310  {
1311  found = true;
1312  /// We need to do a full shell without any elements that are closer
1313  /// to make sure there no closer elements in neighboring searchgrid cells
1314  for (index_type i = bi; i <= ei; i++)
1315  {
1316  if (i < 0 || i > ni) continue;
1317  for (index_type j = bj; j <= ej; j++)
1318  {
1319  if (j < 0 || j > nj) continue;
1320  for (index_type k = bk; k <= ek; k++)
1321  {
1322  if (k < 0 || k > nk) continue;
1323  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
1324  {
1325  if (elem_grid_->min_distance_squared(p, i, j, k) < dmin)
1326  {
1327  found = false;
1328  typename SearchGridT<index_type>::iterator it, eit;
1329  elem_grid_->lookup_ijk(it,eit, i, j, k);
1330 
1331  while (it != eit)
1332  {
1333  Core::Geometry::Point p1(points_[faces_[(*it)*4]]);
1334  Core::Geometry::Point p2(points_[faces_[(*it)*4+1]]);
1335  Core::Geometry::Point p3(points_[faces_[(*it)*4+2]]);
1336  Core::Geometry::Point p4(points_[faces_[(*it)*4+3]]);
1337  Core::Geometry::Vector v0(p2-p1);
1338  v0.normalize();
1339  Core::Geometry::Vector v1(p3-p2);
1340  v1.normalize();
1341  Core::Geometry::Vector n=Cross(v0,v1);
1342  double dot=Dot(n,repelnormal);
1343  if (dot > mindotfromrepel &&
1344  (((p1-repelpos).length2() < mindistfromrepelsquared) ||
1345  ((p2-repelpos).length2() < mindistfromrepelsquared) ||
1346  ((p3-repelpos).length2() < mindistfromrepelsquared) ||
1347  ((p4-repelpos).length2() < mindistfromrepelsquared)))
1348  { ++it; continue; }
1349  Core::Geometry::Point r;
1350  est_closest_point_on_quad(r, p, p1, p2, p3, p4);
1351  const double dtmp = (p - r).length2();
1352  if (dtmp < dmin)
1353  {
1354  found_one = true;
1355  result = r;
1356  face = INDEX(*it);
1357  dmin = dtmp;
1358 
1359  if (dmin < epsilon2_)
1360  {
1361  /// As we computed an estimate, we use the Newton's method in the basis functions
1362  /// compute a more refined solution. This function may slow down computation.
1363  /// This piece of code will calculate the coordinates in the local element framework
1364  /// (the newton's method of finding a minimum), then it will project this back
1365  /// THIS CODE SHOULD BE FURTHER OPTIMIZED
1366 
1367  ElemData ed(*this,face);
1368  basis_.get_coords(coords,result,ed);
1369 
1370  result = basis_.interpolate(coords,ed);
1371  dmin = (result-p).length2();
1372 
1373  pdist = sqrt(dmin);
1374  return (true);
1375  }
1376  }
1377 
1378  ++it;
1379  }
1380  }
1381  }
1382  }
1383  }
1384  }
1385  bi--;ei++;
1386  bj--;ej++;
1387  bk--;ek++;
1388  }
1389  while (!found);
1390 
1391  if (!found_one) return (false);
1392 
1393  /// As we computed an estimate, we use the Newton's method in the basis functions
1394  /// compute a more refined solution. This function may slow down computation.
1395  /// This piece of code will calculate the coordinates in the local element framework
1396  /// (the newton's method of finding a minimum), then it will project this back
1397  /// THIS CODE SHOULD BE FURTHER OPTIMIZED
1398 
1399  ElemData ed(*this,face);
1400  basis_.get_coords(coords,result,ed);
1401 
1402  result = basis_.interpolate(coords,ed);
1403  dmin = (result-p).length2();
1404 
1405  pdist = sqrt(dmin);
1406 
1407  return (true);
1408  }
1409 
1410  template <class INDEX>
1411  bool find_closest_elem(double& pdist,
1412  Core::Geometry::Point &result,
1413  INDEX &elem,
1414  const Core::Geometry::Point &p) const
1415  {
1416  StackVector<double,2> coords;
1417  return(find_closest_elem(pdist,result,coords,elem,p,-1.0));
1418  }
1419 
1420  template <class ARRAY>
1421  bool find_closest_elems(double& pdist, Core::Geometry::Point &result,
1422  ARRAY &elems, const Core::Geometry::Point &p) const
1423  {
1424  elems.clear();
1425 
1426  typename Elem::size_type sz; size(sz);
1427 
1428  /// If there are no nodes we cannot find the closest one
1429  if (sz == 0) return (false);
1430 
1431  // Walking the grid like this works really well if we're near the
1432  // surface. It's degenerately bad if for example the point is
1433  // placed in the center of a sphere (because then we still have to
1434  // test all the faces, but with the grid overhead and triangle
1435  // duplication as well).
1436  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
1437  "QuadSurfMesh::find_closest_elems requires synchronize(ELEM_LOCATE_E).")
1438 
1439  // get grid sizes
1440  const size_type ni = elem_grid_->get_ni()-1;
1441  const size_type nj = elem_grid_->get_nj()-1;
1442  const size_type nk = elem_grid_->get_nk()-1;
1443 
1444  // Convert to grid coordinates.
1445  index_type bi, bj, bk, ei, ej, ek;
1446  elem_grid_->unsafe_locate(bi, bj, bk, p);
1447 
1448  // Clamp to closest point on the grid.
1449  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
1450  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
1451  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
1452 
1453  ei = bi; ej = bj; ek = bk;
1454 
1455  double dmin = DBL_MAX;
1456 
1457  bool found = true;
1458  do
1459  {
1460  found = true;
1461  /// This looks incorrect - but it is correct
1462  /// We need to do a full shell without any elements that are closer
1463  /// to make sure there no closer elements
1464  for (index_type i = bi; i <= ei; i++)
1465  {
1466  if (i < 0 || i > ni) continue;
1467  for (index_type j = bj; j <= ej; j++)
1468  {
1469  if (j < 0 || j > nj) continue;
1470  for (index_type k = bk; k <= ek; k++)
1471  {
1472  if (k < 0 || k > nk) continue;
1473  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
1474  {
1475  if (elem_grid_->min_distance_squared(p, i, j, k) < dmin)
1476  {
1477  found = false;
1478  typename SearchGridT<index_type>::iterator it, eit;
1479  elem_grid_->lookup_ijk(it,eit, i, j, k);
1480 
1481  while (it != eit)
1482  {
1483  Core::Geometry::Point rtmp;
1484  const index_type idx = (*it) * 4;
1485  est_closest_point_on_quad(rtmp, p,
1486  points_[faces_[idx ]],
1487  points_[faces_[idx+1]],
1488  points_[faces_[idx+2]],
1489  points_[faces_[idx+3]]);
1490  const double dtmp = (p - rtmp).length2();
1491 
1492  if (dtmp < dmin - epsilon2_)
1493  {
1494  elems.clear();
1495  result = rtmp;
1496  elems.push_back(*it);
1497  dmin = dtmp;
1498  }
1499  else if (dtmp < dmin + epsilon2_)
1500  {
1501  elems.push_back(typename ARRAY::value_type(*it));
1502  }
1503  ++it;
1504  }
1505  }
1506  }
1507  }
1508  }
1509  }
1510  bi--;ei++;
1511  bj--;ej++;
1512  bk--;ek++;
1513  }
1514  while ((!found)||(dmin == DBL_MAX)) ;
1515 
1516 
1517  // As we computed an estimate, we use the Newton's method in the basis functions
1518  // compute a more refined solution. This function may slow down computation.
1519  // This piece of code will calculate the coordinates in the local element framework
1520  // (the newton's method of finding a minimum), then it will project this back
1521  // THIS CODE SHOULD BE FURTHER OPTIMIZED
1522 
1523  if (elems.size() == 1)
1524  {
1525  // if the number of faces is more then one the point we found is located
1526  // on the node or on the edge, which means the estimate is correct.
1527  std::vector<double> coords;
1528  ElemData ed(*this,elems[0]);
1529  basis_.get_coords(coords,result,ed);
1530  result = basis_.interpolate(coords,ed);
1531  dmin = (result-p).length2();
1532  }
1533 
1534  pdist = sqrt(dmin);
1535  return (true);
1536  }
1537 
1538 
1539  double get_epsilon() const
1540  { return (epsilon_); }
1541 
1542  ///////////////////////////////////////////////////
1543  // STATIC VARIABLES AND FUNCTIONS
1544 
1545  /// Export this class using the old Pio system
1546  virtual void io(Piostream&);
1547 
1548  /// This ID is created as soon as this class will be instantiated
1549  static PersistentTypeID quadsurfmesh_typeid;
1550 
1551  /// Core functionality for getting the name of a templated mesh class
1552  static const std::string type_name(int n = -1);
1553  virtual std::string dynamic_type_name() const { return quadsurfmesh_typeid.type; }
1554 
1555  /// Type description, used for finding names of the mesh class for
1556  /// dynamic compilation purposes. Some of this should be obsolete
1557  virtual const TypeDescription *get_type_description() const;
1558  static const TypeDescription* node_type_description();
1559  static const TypeDescription* edge_type_description();
1560  static const TypeDescription* face_type_description();
1561  static const TypeDescription* cell_type_description();
1562  static const TypeDescription* elem_type_description()
1563  { return face_type_description(); }
1564 
1565  /// This function returns a maker for Pio.
1566  static Persistent *maker() { return new QuadSurfMesh<Basis>(); }
1567  /// This function returns a handle for the virtual interface.
1568  static MeshHandle mesh_maker() { return boost::make_shared<QuadSurfMesh<Basis>>(); }
1569 
1570 
1571  //////////////////////////////////////////////////////////////////
1572  // Mesh specific functions (these are not implemented in every mesh)
1573 
1574  template <class INDEX>
1575  bool inside(INDEX idx, const Core::Geometry::Point &p) const
1576  {
1577 
1578  typename Node::array_type nodes;
1579  get_nodes_from_elem(nodes,idx);
1580 
1581  Core::Geometry::BBox bbox;
1582  bbox.extend(points_[nodes[0]]);
1583  bbox.extend(points_[nodes[1]]);
1584  bbox.extend(points_[nodes[2]]);
1585  bbox.extend(points_[nodes[3]]);
1586  bbox.extend(epsilon_);
1587 
1588  if (bbox.inside(p))
1589  {
1590  StackVector<double,2> coords;
1591  ElemData ed(*this, idx);
1592  if(basis_.get_coords(coords, p, ed)) return (true);
1593  }
1594 
1595  return (false);
1596  }
1597 
1598  // Extra functionality needed by this specific geometry.
1599  typename Node::index_type add_find_point(const Core::Geometry::Point &p,
1600  double err = 1.0e-3);
1601  typename Elem::index_type add_quad(typename Node::index_type a,
1602  typename Node::index_type b,
1603  typename Node::index_type c,
1604  typename Node::index_type d);
1605  typename Elem::index_type add_quad(const Core::Geometry::Point &p0, const Core::Geometry::Point &p1,
1606  const Core::Geometry::Point &p2, const Core::Geometry::Point &p3);
1607 
1608 protected:
1609 
1610  //////////////////////////////////////////////////////////////
1611  // These functions are templates and are used to define the
1612  // dynamic compilation interface and the virtual interface
1613  // as they both use different datatypes as indices and arrays
1614  // the following functions have been templated and are inlined
1615  // at the places where they are needed.
1616  //
1617  // Secondly these templates allow for the use of the stack vector
1618  // as well as the STL vector. When an algorithm supports non linear
1619  // functions an STL vector is a better choice, in the other cases
1620  // often a StackVector is enough (The latter improves performance).
1621 
1622 
1623  template<class ARRAY, class INDEX>
1624  inline void get_nodes_from_edge(ARRAY& array, INDEX i) const
1625  {
1626  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1627  "QuadSurfMesh: Must call synchronize EDGES_E on QuadSurfMesh first");
1628 
1629  static int table[4][2] = { {0, 1}, {1, 2}, {2, 3}, {3, 0} };
1630 
1631  const index_type idx = edges_[i][0];
1632  const index_type off = idx % 4;
1633  const index_type node = idx - off;
1634  array.resize(2);
1635  array[0] = static_cast<typename ARRAY::value_type>(faces_[node + table[off][0]]);
1636  array[1] = static_cast<typename ARRAY::value_type>(faces_[node + table[off][1]]);
1637  }
1638 
1639 
1640  template<class ARRAY, class INDEX>
1641  inline void get_nodes_from_face(ARRAY& array, INDEX idx) const
1642  {
1643  array.resize(4);
1644  array[0] = static_cast<typename ARRAY::value_type>(faces_[idx * 4 + 0]);
1645  array[1] = static_cast<typename ARRAY::value_type>(faces_[idx * 4 + 1]);
1646  array[2] = static_cast<typename ARRAY::value_type>(faces_[idx * 4 + 2]);
1647  array[3] = static_cast<typename ARRAY::value_type>(faces_[idx * 4 + 3]);
1648  order_face_nodes(array[0],array[1],array[2],array[3]);
1649  }
1650 
1651  template<class ARRAY, class INDEX>
1652  inline void get_nodes_from_elem(ARRAY& array, INDEX idx) const
1653  {
1654  get_nodes_from_face(array,idx);
1655  }
1656 
1657  template<class ARRAY, class INDEX>
1658  inline void get_edges_from_face(ARRAY& array, INDEX idx) const
1659  {
1660  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1661  "QuadSurfMesh: Must call synchronize EDGES_E on QuadSurfMesh first");
1662 
1663  array.clear();
1664  //array.reserve(4);
1665  index_type he;
1666  int i = 0;
1667  he = static_cast<typename ARRAY::value_type>(halfedge_to_edge_[idx * 4 ]);
1668  if (he >=0) array[i++] = he;
1669  he = static_cast<typename ARRAY::value_type>(halfedge_to_edge_[idx * 4 + 1]);
1670  if (he >=0) array[i++] = he;
1671  he = static_cast<typename ARRAY::value_type>(halfedge_to_edge_[idx * 4 + 2]);
1672  if (he >=0) array[i++] = he;
1673  he = static_cast<typename ARRAY::value_type>(halfedge_to_edge_[idx * 4 + 3]);
1674  if (he >=0) array[i++] = he;
1675  }
1676 
1677 
1678  template<class ARRAY, class INDEX>
1679  inline void get_edges_from_elem(ARRAY& array, INDEX idx) const
1680  {
1681  get_edges_from_face(array,idx);
1682  }
1683 
1684 
1685  template <class ARRAY, class INDEX>
1686  inline void get_faces_from_node(ARRAY& array, INDEX idx) const
1687  {
1688  ASSERTMSG(synchronized_ & Mesh::NODE_NEIGHBORS_E,
1689  "QuadSurfMesh: Must call synchronize NODE_NEIGHBORS_E on QuadSurfMesh first");
1690 
1691  array.resize(node_neighbors_[idx].size());
1692  for (size_t i = 0; i < node_neighbors_[idx].size(); ++i)
1693  array[i] = static_cast<typename ARRAY::value_type>(node_neighbors_[idx][i]);
1694  }
1695 
1696 
1697  template <class ARRAY, class INDEX>
1698  inline void get_faces_from_edge(ARRAY& array, INDEX idx) const
1699  {
1700  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1701  "Must call synchronize EDGES_E on QuadSurfMesh first");
1702 
1703  const std::vector<index_type>& faces = edges_[idx];
1704 
1705  // clear array
1706  array.clear();
1707  // conservative estimate
1708  array.reserve(faces.size());
1709 
1710  size_type fs = faces.size();
1711  for (index_type i=0; i<fs; i++)
1712  {
1713  array.push_back((faces[i]>>2));
1714  }
1715  }
1716 
1717 
1718  template<class ARRAY, class INDEX>
1719  inline void get_edges_from_node(ARRAY& array, INDEX idx) const
1720  {
1721  ASSERTMSG(synchronized_ & Mesh::NODE_NEIGHBORS_E,
1722  "Must call synchronize NODE_NEIGHBORS_E on QuadSurfMesh first");
1723 
1724  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1725  "Must call synchronize EDGES_E on QuadSurfMesh first");
1726  // Get the table of faces that are connected to the two nodes
1727  Core::Thread::Guard nn(synchronize_lock_.get());
1728 
1729  const std::vector<typename Face::index_type>& faces = node_neighbors_[idx];
1730  array.clear();
1731 
1732  typename ARRAY::value_type edge;
1733  for (index_type i=0; i<static_cast<size_type>(faces.size()); i++)
1734  {
1735  for (index_type j=0; j<4; j++)
1736  {
1737  edge = halfedge_to_edge_[faces[i]*4+j];
1738  size_t k=0;
1739  for (; k<array.size(); k++)
1740  if (array[k] == edge) break;
1741  if (k == array.size()) array.push_back(edge);
1742  }
1743  }
1744  }
1745 
1746  template <class ARRAY, class INDEX>
1747  inline void set_nodes_by_elem(ARRAY &array, INDEX idx)
1748  {
1749  for (index_type n = 0; n < 4; ++n)
1750  faces_[idx * 4 + n] = static_cast<index_type>(array[n]);
1751  }
1752 
1753  /// This function has been rewritten to allow for non manifold surfaces to be
1754  /// handled ok.
1755  template <class INDEX1, class INDEX2>
1756  bool get_elem_neighbor(INDEX1 &neighbor, INDEX1 elem, INDEX2 delem) const
1757  {
1758  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1759  "Must call synchronize EDGES_E on QuadSurfMesh first");
1760 
1761  const std::vector<index_type>& faces = edges_[delem];
1762 
1763  size_type fs = faces.size();
1764  for (index_type i=0; i<fs; i++)
1765  {
1766  index_type face = (faces[i]>>2);
1767  if (face != elem)
1768  {
1769  neighbor = face;
1770  return (true);
1771  }
1772  }
1773  return (false);
1774  }
1775 
1776  template <class ARRAY,class INDEX1, class INDEX2>
1777  bool get_elem_neighbors(ARRAY &array, INDEX1 elem, INDEX2 delem) const
1778  {
1779  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1780  "Must call synchronize EDGES_E on QuadSurfMesh first");
1781 
1782  // Find the two nodes that make up the edge
1783  const std::vector<index_type>& faces = edges_[delem];
1784 
1785  // clear array
1786  array.clear();
1787 
1788  size_type fs = faces.size();
1789  for (index_type i=0; i<fs; i++)
1790  {
1791  index_type face = (faces[i]>>2);
1792  if (face != elem)
1793  {
1794  array.push_back(face);
1795  }
1796  }
1797  return (array.size() > 0);
1798  }
1799 
1800  template <class ARRAY, class INDEX>
1801  void get_elem_neighbors(ARRAY &array, INDEX idx) const
1802  {
1803  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1804  "Must call synchronize EDGES_E on QuadSurfMesh first");
1805 
1806  typename Edge::array_type edges;
1807  get_edges_from_face(edges, idx);
1808 
1809  array.clear();
1810  array.reserve(edges.size());
1811 
1812  for (index_type i=0; i<static_cast<index_type>(edges.size()); i++)
1813  {
1814  typename ARRAY::value_type nbor;
1815  if (get_elem_neighbor(nbor, idx, edges[i]))
1816  {
1817  array.push_back(nbor);
1818  }
1819  }
1820  }
1821 
1822 
1823  template <class ARRAY, class INDEX>
1824  void get_node_neighbors(ARRAY &array, INDEX idx) const
1825  {
1826  ASSERTMSG(synchronized_ & Mesh::NODE_NEIGHBORS_E,
1827  "Must call synchronize NODE_NEIGHBORS_E on QuadSurfMesh first");
1828  // clear old contents
1829  array.clear();
1830 
1831  // Get all the neighboring elements
1832  typename Node::array_type nodes;
1833  const std::vector<typename Face::index_type>& faces = node_neighbors_[idx];
1834  // Make a conservative estimate of the number of node neighbors
1835  array.reserve(2*faces.size());
1836 
1837  for (index_type i=0;i<static_cast<index_type>(faces.size());i++)
1838  {
1839  index_type base = faces[i]*4;
1840  for (index_type j=0;j<4;j++)
1841  {
1842  if (faces_[base+j] == idx) continue;
1843  size_t k=0;
1844  for (;k<array.size();k++)
1845  if (static_cast<typename ARRAY::value_type>(faces_[base+j]) == array[k]) break;
1846  if (k==array.size())
1847  array.push_back(static_cast<typename ARRAY::value_type>(faces_[base+j]));
1848  }
1849  }
1850  }
1851 
1852 // That would increase the search speed of node localization
1853 
1854  template <class INDEX>
1855  inline bool locate_node(INDEX &node, const Core::Geometry::Point &p) const
1856  {
1857  typename Node::size_type sz; size(sz);
1858 
1859  /// If there are no nodes we cannot find a closest point
1860  if (sz == 0) return (false);
1861 
1862  /// Check first guess
1863  if (node >= 0 && node < sz)
1864  {
1865  if ((p - points_[node]).length2() < epsilon2_) return (true);
1866  }
1867 
1868  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
1869  "QuadSurfMesh::locate_node requires synchronize(NODE_LOCATE_E).")
1870 
1871  // get grid sizes
1872  const size_type ni = node_grid_->get_ni()-1;
1873  const size_type nj = node_grid_->get_nj()-1;
1874  const size_type nk = node_grid_->get_nk()-1;
1875 
1876  // Convert to grid coordinates.
1877  index_type bi, bj, bk, ei, ej, ek;
1878  node_grid_->unsafe_locate(bi, bj, bk, p);
1879 
1880  // Clamp to closest point on the grid.
1881  if (bi > ni) bi =ni; if (bi < 0) bi = 0;
1882  if (bj > nj) bj =nj; if (bj < 0) bj = 0;
1883  if (bk > nk) bk =nk; if (bk < 0) bk = 0;
1884 
1885  ei = bi;
1886  ej = bj;
1887  ek = bk;
1888 
1889  double dmin = DBL_MAX;
1890  bool found = true;
1891  do
1892  {
1893  found = true;
1894  /// This looks incorrect - but it is correct
1895  /// We need to do a full shell without any elements that are closer
1896  /// to make sure there no closer elements in neighboring searchgrid cells
1897 
1898  for (index_type i = bi; i <= ei; i++)
1899  {
1900  if (i < 0 || i > ni) continue;
1901  for (index_type j = bj; j <= ej; j++)
1902  {
1903  if (j < 0 || j > nj) continue;
1904  for (index_type k = bk; k <= ek; k++)
1905  {
1906  if (k < 0 || k > nk) continue;
1907  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
1908  {
1909  if (node_grid_->min_distance_squared(p, i, j, k) < dmin)
1910  {
1911  found = false;
1912  typename SearchGridT<index_type>::iterator it, eit;
1913  node_grid_->lookup_ijk(it,eit, i, j, k);
1914 
1915  while (it != eit)
1916  {
1917  const Core::Geometry::Point point = points_[*it];
1918  const double dist = (p-point).length2();
1919 
1920  if (dist < dmin)
1921  {
1922  node = INDEX(*it);
1923  dmin = dist;
1924 
1925  if (dist < epsilon2_) return (true);
1926  }
1927  ++it;
1928  }
1929  }
1930  }
1931  }
1932  }
1933  }
1934  bi--;ei++;
1935  bj--;ej++;
1936  bk--;ek++;
1937  }
1938  while ((!found)||(dmin == DBL_MAX)) ;
1939 
1940  return (true);
1941  }
1942 
1943  /// This is currently implemented as an exhaustive search
1944  template <class INDEX>
1945  inline bool locate_edge(INDEX &loc, const Core::Geometry::Point &p) const
1946  {
1947  ASSERTMSG(synchronized_ & EDGES_E,
1948  "QuadSurfMesh::locate_edge requires synchronize(EDGES_E).")
1949 
1950  typename Edge::iterator bi, ei;
1951  typename Node::array_type nodes;
1952 
1953  begin(bi);
1954  end(ei);
1955 
1956  bool found = false;
1957  double mindist = 0.0;
1958  while (bi != ei)
1959  {
1960  get_nodes(nodes,*bi);
1961  const double dist = distance_to_line2(p, points_[nodes[0]],
1962  points_[nodes[1]],epsilon_);
1963  if (!found || dist < mindist)
1964  {
1965  loc = static_cast<INDEX>(*bi);
1966  mindist = dist;
1967  found = true;
1968  }
1969  ++bi;
1970  }
1971  return (found);
1972  }
1973 
1974 
1975  template <class INDEX>
1976  inline bool locate_elem(INDEX &elem, const Core::Geometry::Point &p) const
1977  {
1978  if (basis_.polynomial_order() > 1) return elem_locate(elem, *this, p);
1979 
1980  typename Elem::size_type sz; size(sz);
1981 
1982  /// If there are no nodes we cannot find a closest point
1983  if (sz == 0) return (false);
1984 
1985  /// Check whether the estimate given in idx is the point we are looking for
1986  if ((elem > 0)&&(elem < sz))
1987  {
1988  if (inside(elem,p)) return (true);
1989  }
1990 
1991  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
1992  "QuadSurfMesh::locate_node requires synchronize(ELEM_LOCATE_E).")
1993 
1994  typename SearchGridT<index_type>::iterator it, eit;
1995  if (elem_grid_->lookup(it,eit, p))
1996  {
1997  while (it != eit)
1998  {
1999  if (inside(typename Elem::index_type(*it), p))
2000  {
2001  elem = static_cast<INDEX>(*it);
2002  return (true);
2003  }
2004  ++it;
2005  }
2006  }
2007  return (false);
2008  }
2009 
2010  template <class ARRAY>
2011  inline bool locate_elems(ARRAY &array, const Core::Geometry::BBox &b) const
2012  {
2013 
2014  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
2015  "QuadSurfMesh::locate_elems requires synchronize(ELEM_LOCATE_E).")
2016 
2017  array.clear();
2018  index_type is,js,ks;
2019  index_type ie,je,ke;
2020  elem_grid_->locate_clamp(is,js,ks,b.min());
2021  elem_grid_->locate_clamp(ie,je,ke,b.max());
2022  for (index_type i=is; i<=ie;i++)
2023  for (index_type j=js; j<je;j++)
2024  for (index_type k=ks; k<ke; k++)
2025  {
2026  typename SearchGridT<index_type>::iterator it, eit;
2027  elem_grid_->lookup_ijk(it, eit, i, j, k);
2028  while (it != eit)
2029  {
2030  size_t p=0;
2031  for (;p<array.size();p++) if (array[p] == typename ARRAY::value_type(*it)) break;
2032  if (p == array.size()) array.push_back(typename ARRAY::value_type(*it));
2033  ++it;
2034  }
2035  }
2036 
2037  return (array.size() > 0);
2038  }
2039 
2040 
2041  template <class INDEX, class ARRAY>
2042  inline bool locate_elem(INDEX &elem, ARRAY& coords, const Core::Geometry::Point &p) const
2043  {
2044  if (basis_.polynomial_order() > 1) return elem_locate(elem, *this, p);
2045 
2046  typename Elem::size_type sz; size(sz);
2047 
2048  /// If there are no nodes we cannot find a closest point
2049  if (sz == 0) return (false);
2050 
2051  /// Check whether the estimate given in idx is the point we are looking for
2052  if ((elem > 0)&&(elem < sz))
2053  {
2054  if (inside(elem,p))
2055  {
2056  ElemData ed(*this, elem);
2057  basis_.get_coords(coords, p, ed);
2058  return (true);
2059  }
2060  }
2061 
2062  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
2063  "QuadSurfMesh::locate_node requires synchronize(ELEM_LOCATE_E).")
2064 
2065  typename SearchGridT<index_type>::iterator it, eit;
2066  if (elem_grid_->lookup(it,eit, p))
2067  {
2068  while (it != eit)
2069  {
2070  if (inside(typename Elem::index_type(*it), p))
2071  {
2072  elem = static_cast<INDEX>(*it);
2073  ElemData ed(*this, elem);
2074  basis_.get_coords(coords, p, ed);
2075  return (true);
2076  }
2077  ++it;
2078  }
2079  }
2080  return (false);
2081  }
2082 
2083 
2084  template <class INDEX>
2085  inline void get_node_center(Core::Geometry::Point &p, INDEX idx) const
2086  {
2087  p = points_[idx];
2088  }
2089 
2090  template <class INDEX>
2091  inline void get_edge_center(Core::Geometry::Point &result, INDEX idx) const
2092  {
2093  typename Node::array_type arr;
2094  get_nodes_from_edge(arr, idx);
2095  result = points_[arr[0]];
2096  result += points_[arr[1]];
2097  result *= 0.5;
2098  }
2099 
2100 
2101  template <class INDEX>
2102  inline void get_face_center(Core::Geometry::Point &p, INDEX idx) const
2103  {
2104  // NEED TO OPTIMIZE THIS ONE
2105  typename Node::array_type nodes;
2106  get_nodes_from_face(nodes, idx);
2107  ASSERT(nodes.size() == 4);
2108  typename Node::array_type::iterator nai = nodes.begin();
2109  p = points_[*nai];
2110  ++nai;
2111  while (nai != nodes.end())
2112  {
2113  const Core::Geometry::Point &pp = points_[*nai];
2114  p += pp;
2115  ++nai;
2116  }
2117  p *= (1.0 / 4.0);
2118  }
2119 
2120 
2121  const Core::Geometry::Point &point(typename Node::index_type i) const { return points_[i]; }
2122 
2123  // These require the synchronize_lock_ to be held before calling.
2124  void compute_edges();
2125  void compute_normals();
2126  void compute_node_neighbors();
2127 
2128  void compute_node_grid();
2129  void compute_elem_grid();
2130  void compute_bounding_box();
2131 
2132  /// Used to recompute data for individual cells.
2133  void insert_elem_into_grid(typename Elem::index_type ci);
2134  void remove_elem_from_grid(typename Elem::index_type ci);
2135 
2136  void insert_node_into_grid(typename Node::index_type ci);
2137  void remove_node_from_grid(typename Node::index_type ci);
2138 
2139  template <class NODE>
2140  bool order_face_nodes(NODE& n1,NODE& n2, NODE& n3, NODE& n4) const
2141  {
2142  // Check for degenerate or misformed face
2143  // Opposite faces cannot be equal
2144  if ((n1 == n3)||(n2==n4)) return (false);
2145 
2146  // Face must have three unique identifiers otherwise it was condition
2147  // n1==n3 || n2==n4 would be met.
2148 
2149  if (n1==n2)
2150  {
2151  if (n3==n4) return (false); // this is a line not a face
2152  NODE t;
2153  // shift one position to left
2154  t = n1; n1 = n2; n2 = n3; n3 = n4; n4 = t;
2155  return (true);
2156  }
2157  else if (n2 == n3)
2158  {
2159  if (n1==n4) return (false); // this is a line not a face
2160  NODE t;
2161  // shift two positions to left
2162  t = n1; n1 = n3; n3 = t; t = n2; n2 = n4; n4 = t;
2163  return (true);
2164  }
2165  else if (n3 == n4)
2166  {
2167  NODE t;
2168  // shift one positions to right
2169  t = n4; n4 = n3; n3 = n2; n2 = n1; n1 = t;
2170  return (true);
2171  }
2172  else if (n4 == n1)
2173  {
2174  // proper order
2175  return (true);
2176  }
2177  else
2178  {
2179  if ((n1 < n2)&&(n1 < n3)&&(n1 < n4))
2180  {
2181  // proper order
2182  return (true);
2183  }
2184  else if ((n2 < n3)&&(n2 < n4))
2185  {
2186  NODE t;
2187  // shift one position to left
2188  t = n1; n1 = n2; n2 = n3; n3 = n4; n4 = t;
2189  return (true);
2190  }
2191  else if (n3 < n4)
2192  {
2193  NODE t;
2194  // shift two positions to left
2195  t = n1; n1 = n3; n3 = t; t = n2; n2 = n4; n4 = t;
2196  return (true);
2197  }
2198  else
2199  {
2200  NODE t;
2201  // shift one positions to right
2202  t = n4; n4 = n3; n3 = n2; n2 = n1; n1 = t;
2203  return (true);
2204  }
2205  }
2206  }
2207 
2208 
2209  index_type next(index_type i) { return ((i%4)==3) ? (i-3) : (i+1); }
2210  index_type prev(index_type i) { return ((i%4)==0) ? (i+3) : (i-1); }
2211 
2212  /// array with all the points
2213  std::vector<Core::Geometry::Point> points_;
2214  /// array with the four nodes that make up a face
2215  std::vector<index_type> faces_;
2216 
2217  /// FOR EDGE -> NODES
2218  /// array with information from edge number (unique ones) to the node numbers
2219  /// this one refers to the first node
2220  std::vector<std::vector<index_type> > edges_;
2221  /// FOR FACES -> EDGE NUMBER
2222  /// array with information from halfedge (computed directly from face) to the edge number
2223  std::vector<index_type> halfedge_to_edge_; // halfedge->edge map
2224 
2225  typedef std::vector<std::vector<typename Elem::index_type> > NodeNeighborMap;
2226  NodeNeighborMap node_neighbors_;
2227 
2228  std::vector<Core::Geometry::Vector> normals_; /// normalized per node
2229  boost::shared_ptr<SearchGridT<index_type> > node_grid_; /// Lookup grid for nodes
2230  boost::shared_ptr<SearchGridT<index_type> > elem_grid_; /// Lookup grid for elements
2231 
2232  // Lock and Condition Variable for hand shaking
2233  mutable Core::Thread::Mutex synchronize_lock_;
2234  Core::Thread::ConditionVariable synchronize_cond_;
2235 
2236  // Which tables have been computed
2237  mask_type synchronized_;
2238  // Which tables are currently being computed
2239  mask_type synchronizing_;
2240 
2241  Basis basis_; /// Basis for interpolation
2242 
2243  Core::Geometry::BBox bbox_;
2244  double epsilon_; /// epsilon for calculations 1e-8*diagonal bounding box
2245  double epsilon2_; /// square of epsilon for squared distance comparisons
2246 
2247  /// Pointer to virtual interface
2248  boost::shared_ptr<VMesh> vmesh_; /// Virtual function table
2249 
2250 #ifdef HAVE_HASH_MAP
2251  struct edgehash
2252  {
2253  size_t operator()(const std::pair<index_type, index_type> &a) const
2254  {
2255 #if defined(__ECC) || defined(_MSC_VER)
2256  hash_compare<int> hasher;
2257 #else
2258  hash<int> hasher;
2259 #endif
2260  return hasher(static_cast<int>(hasher(a.first) + a.second));
2261  }
2262 #if defined(__ECC) || defined(_MSC_VER)
2263 
2264  static const size_t bucket_size = 4;
2265  static const size_t min_buckets = 8;
2266 
2267  bool operator()(const std::pair<index_type, index_type> &a, const std::pair<index_type, index_type> &b) const
2268  {
2269  return a.first < b.first || a.first == b.first && a.second < b.second;
2270  }
2271 #endif
2272  };
2273 
2274  struct edgecompare
2275  {
2276  bool operator()(const std::pair<index_type, index_type> &a, const std::pair<index_type, index_type> &b) const
2277  {
2278  return a.first == b.first && a.second == b.second;
2279  }
2280  };
2281 
2282 #else // HAVE_HASH_MAP
2283 
2284  struct edgecompare
2285  {
2286  bool operator()(const std::pair<index_type, index_type> &a, const std::pair<index_type, index_type> &b) const
2287  {
2288  return a.first < b.first || a.first == b.first && a.second < b.second;
2289  }
2290  };
2291 #endif // HAVE_HASH_MAP
2292 
2293 };
2294 
2295 
2296 template <class Basis>
2297 PersistentTypeID
2298 QuadSurfMesh<Basis>::quadsurfmesh_typeid(type_name(-1), "Mesh",
2299  QuadSurfMesh<Basis>::maker);
2300 
2301 template <class Basis>
2302 const std::string
2303 QuadSurfMesh<Basis>::type_name(int n)
2304 {
2305  ASSERT((n >= -1) && n <= 1);
2306  if (n == -1)
2307  {
2308  static const std::string name = TypeNameGenerator::make_template_id(type_name(0), type_name(1));
2309  return name;
2310  }
2311  else if (n == 0)
2312  {
2313  static const std::string nm("QuadSurfMesh");
2314  return nm;
2315  }
2316  else
2317  {
2318  return find_type_name((Basis *)0);
2319  }
2320 }
2321 
2322 
2323 template <class Basis>
2324 QuadSurfMesh<Basis>::QuadSurfMesh()
2325  : points_(0),
2326  faces_(0),
2327  edges_(0),
2328  normals_(0),
2329  synchronize_lock_("QuadSurfMesh lock"),
2330  synchronize_cond_("QuadSurfMesh condition variable"),
2331  synchronized_(Mesh::NODES_E | Mesh::FACES_E | Mesh::CELLS_E),
2332  synchronizing_(0),
2333  epsilon_(0.0),
2334  epsilon2_(0.0)
2335 {
2336  DEBUG_CONSTRUCTOR("QuadSurfMesh")
2337 
2338  /// Initialize the virtual interface when the mesh is created
2339  vmesh_.reset(CreateVQuadSurfMesh(this));
2340 }
2341 
2342 
2343 template <class Basis>
2344 QuadSurfMesh<Basis>::QuadSurfMesh(const QuadSurfMesh &copy)
2345  : Mesh(copy),
2346  points_(0),
2347  faces_(0),
2348  edges_(0),
2349  normals_(0),
2350  synchronize_lock_("QuadSurfMesh lock"),
2351  synchronize_cond_("QuadSurfMesh condition variable"),
2352  synchronized_(Mesh::NODES_E | Mesh::FACES_E | Mesh::CELLS_E),
2353  synchronizing_(0),
2354  epsilon_(0.0),
2355  epsilon2_(0.0)
2356 {
2357  DEBUG_CONSTRUCTOR("QuadSurfMesh")
2358 
2359  /// We need to lock before we can copy these as these
2360  /// structures are generate dynamically when they are
2361  /// needed.
2362  Core::Thread::Guard rlock(copy.synchronize_lock_.get());
2363 
2364  points_ = copy.points_;
2365  faces_ = copy.faces_;
2366 
2367  /// Create a new virtual interface for this copy
2368  /// all pointers have changed hence create a new
2369  /// virtual interface class
2370  vmesh_.reset(CreateVQuadSurfMesh(this));
2371 }
2372 
2373 
2374 template <class Basis>
2375 QuadSurfMesh<Basis>::~QuadSurfMesh()
2376 {
2377  DEBUG_DESTRUCTOR("QuadSurfMesh")
2378 }
2379 
2380 
2381 template <class Basis>
2382 Core::Geometry::BBox
2383 QuadSurfMesh<Basis>::get_bounding_box() const
2384 {
2385  Core::Geometry::BBox result;
2386 
2387  // Compute bounding box
2388  typename Node::iterator ni, nie;
2389  begin(ni);
2390  end(nie);
2391  while (ni != nie)
2392  {
2393  result.extend(points_[*ni]);
2394  ++ni;
2395  }
2396 
2397  return result;
2398 }
2399 
2400 
2401 template <class Basis>
2402 void
2403 QuadSurfMesh<Basis>::get_canonical_transform(Core::Geometry::Transform &t) const
2404 {
2405  t.load_identity();
2406  Core::Geometry::BBox bbox = get_bounding_box();
2407  t.pre_scale(bbox.diagonal());
2408  t.pre_translate(Core::Geometry::Vector(bbox.min()));
2409 }
2410 
2411 template <class Basis>
2412 void
2413 QuadSurfMesh<Basis>::transform(const Core::Geometry::Transform &t)
2414 {
2415  synchronize_lock_.lock();
2416  std::vector<Core::Geometry::Point>::iterator itr = points_.begin();
2417  std::vector<Core::Geometry::Point>::iterator eitr = points_.end();
2418  while (itr != eitr)
2419  {
2420  *itr = t.project(*itr);
2421  ++itr;
2422  }
2423 
2424 
2425 
2426  if (bbox_.valid())
2427  {
2428  bbox_.reset();
2429 
2430  // Compute bounding box
2431  typename Node::iterator ni, nie;
2432  begin(ni);
2433  end(nie);
2434  while (ni != nie)
2435  {
2436  bbox_.extend(point(*ni));
2437  ++ni;
2438  }
2439 
2440  // Compute epsilons associated with the bounding box
2441  epsilon_ = bbox_.diagonal().length()*1e-8;
2442  epsilon2_ = epsilon_*epsilon_;
2443  synchronized_ |= Mesh::BOUNDING_BOX_E;
2444  }
2445 
2446  if (node_grid_) { node_grid_->transform(t); }
2447  if (elem_grid_) { elem_grid_->transform(t); }
2448  synchronize_lock_.unlock();
2449 }
2450 
2451 
2452 template <class Basis>
2453 void
2454 QuadSurfMesh<Basis>::begin(typename QuadSurfMesh::Node::iterator &itr) const
2455 {
2456  ASSERTMSG(synchronized_ & Mesh::NODES_E,
2457  "Must call synchronize NODES_E on QuadSurfMesh first");
2458  itr = 0;
2459 }
2460 
2461 
2462 template <class Basis>
2463 void
2464 QuadSurfMesh<Basis>::end(typename QuadSurfMesh::Node::iterator &itr) const
2465 {
2466  ASSERTMSG(synchronized_ & Mesh::NODES_E,
2467  "Must call synchronize NODES_E on QuadSurfMesh first");
2468  itr = static_cast<size_type>(points_.size());
2469 }
2470 
2471 
2472 template <class Basis>
2473 void
2474 QuadSurfMesh<Basis>::begin(typename QuadSurfMesh::Edge::iterator &itr) const
2475 {
2476  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
2477  "Must call synchronize EDGES_E on QuadSurfMesh first");
2478  itr = 0;
2479 }
2480 
2481 
2482 template <class Basis>
2483 void
2484 QuadSurfMesh<Basis>::end(typename QuadSurfMesh::Edge::iterator &itr) const
2485 {
2486  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
2487  "Must call synchronize EDGES_E on QuadSurfMesh first");
2488  itr = static_cast<size_type>(edges_.size());
2489 }
2490 
2491 
2492 template <class Basis>
2493 void
2494 QuadSurfMesh<Basis>::begin(typename QuadSurfMesh::Face::iterator &itr) const
2495 {
2496  ASSERTMSG(synchronized_ & Mesh::FACES_E,
2497  "Must call synchronize FACES_E on QuadSurfMesh first");
2498  itr = 0;
2499 }
2500 
2501 
2502 template <class Basis>
2503 void
2504 QuadSurfMesh<Basis>::end(typename QuadSurfMesh::Face::iterator &itr) const
2505 {
2506  ASSERTMSG(synchronized_ & Mesh::FACES_E,
2507  "Must call synchronize FACES_E on QuadSurfMesh first");
2508  itr = static_cast<size_type>(faces_.size() / 4);
2509 }
2510 
2511 
2512 template <class Basis>
2513 void
2514 QuadSurfMesh<Basis>::begin(typename QuadSurfMesh::Cell::iterator &itr) const
2515 {
2516  ASSERTMSG(synchronized_ & Mesh::CELLS_E,
2517  "Must call synchronize CELLS_E on QuadSurfMesh first");
2518  itr = 0;
2519 }
2520 
2521 
2522 template <class Basis>
2523 void
2524 QuadSurfMesh<Basis>::end(typename QuadSurfMesh::Cell::iterator &itr) const
2525 {
2526  ASSERTMSG(synchronized_ & Mesh::CELLS_E,
2527  "Must call synchronize CELLS_E on QuadSurfMesh first");
2528  itr = 0;
2529 }
2530 
2531 
2532 template <class Basis>
2533 int
2534 QuadSurfMesh<Basis>::get_weights(const Core::Geometry::Point &p, typename Face::array_type &l,
2535  double *w)
2536 {
2537  typename Face::index_type idx;
2538  if (locate(idx, p))
2539  {
2540  l.resize(1);
2541  l[0] = idx;
2542  w[0] = 1.0;
2543  return 1;
2544  }
2545  return 0;
2546 }
2547 
2548 
2549 template <class Basis>
2550 int
2551 QuadSurfMesh<Basis>::get_weights(const Core::Geometry::Point &p, typename Node::array_type &l,
2552  double *w)
2553 {
2554  typename Face::index_type idx;
2555  if (locate(idx, p))
2556  {
2557  get_nodes(l,idx);
2558  std::vector<double> coords(2);
2559  if (get_coords(coords, p, idx))
2560  {
2561  basis_.get_weights(coords, w);
2562  return basis_.dofs();
2563  }
2564  }
2565  return 0;
2566 }
2567 
2568 
2569 template <class Basis>
2570 void
2571 QuadSurfMesh<Basis>::get_random_point(Core::Geometry::Point &p,
2572  typename Elem::index_type ei,
2573  FieldRNG &rng) const
2574 {
2575  const Core::Geometry::Point &a0 = points_[faces_[ei*4+0]];
2576  const Core::Geometry::Point &a1 = points_[faces_[ei*4+1]];
2577  const Core::Geometry::Point &a2 = points_[faces_[ei*4+2]];
2578 
2579  const Core::Geometry::Point &b0 = points_[faces_[ei*4+2]];
2580  const Core::Geometry::Point &b1 = points_[faces_[ei*4+3]];
2581  const Core::Geometry::Point &b2 = points_[faces_[ei*4+0]];
2582 
2583  const double aarea = Cross(a1 - a0, a2 - a0).length();
2584  const double barea = Cross(b1 - b0, b2 - b0).length();
2585 
2586  if (rng() * (aarea + barea) < aarea)
2587  {
2588  // Fold the quad sample into a triangle.
2589  double u = rng();
2590  double v = rng();
2591  if (u + v > 1.0) { u = 1.0 - u; v = 1.0 - v; }
2592 
2593  // Compute the position of the random point.
2594  p = a0+((a1-a0)*u)+((a2-a0)*v);
2595  }
2596  else
2597  {
2598  // Fold the quad sample into a triangle.
2599  double u = rng();
2600  double v = rng();
2601  if (u + v > 1.0) { u = 1.0 - u; v = 1.0 - v; }
2602 
2603  // Compute the position of the random point.
2604  p = b0+((b1-b0)*u)+((b2-b0)*v);
2605  }
2606 }
2607 
2608 template <class Basis>
2609 bool
2610 QuadSurfMesh<Basis>::synchronize(mask_type sync)
2611 {
2612  // Conversion table
2613  if (sync & (Mesh::DELEMS_E))
2614  { sync |= Mesh::EDGES_E; sync &= ~(Mesh::DELEMS_E); }
2615 
2616  if (sync & Mesh::FIND_CLOSEST_NODE_E)
2617  { sync |= NODE_LOCATE_E; sync &= ~(Mesh::FIND_CLOSEST_NODE_E); }
2618 
2619  if (sync & Mesh::FIND_CLOSEST_ELEM_E)
2620  { sync |= ELEM_LOCATE_E; sync &= ~(Mesh::FIND_CLOSEST_ELEM_E); }
2621 
2622  if (sync & (Mesh::NODE_LOCATE_E|Mesh::ELEM_LOCATE_E)) sync |= Mesh::BOUNDING_BOX_E;
2623  if (sync & Mesh::ELEM_NEIGHBORS_E) sync |= Mesh::NODE_NEIGHBORS_E|Mesh::EDGES_E;
2624 
2625  if (sync & Mesh::NODE_NEIGHBORS_E) sync |= Mesh::EDGES_E;
2626 
2627  // Filter out the only tables available
2628  sync &= (Mesh::EDGES_E|Mesh::NORMALS_E|
2629  Mesh::NODE_NEIGHBORS_E|Mesh::BOUNDING_BOX_E|
2630  Mesh::NODE_LOCATE_E|Mesh::ELEM_LOCATE_E);
2631 
2632  {
2633 
2634  Core::Thread::UniqueLock lock(synchronize_lock_.get());
2635 
2636  // Only sync what hasn't been synched
2637  sync &= (~synchronized_);
2638 
2639 
2640  if (sync == 0)
2641  {
2642  return (true);
2643  }
2644 
2645  if (sync == Mesh::EDGES_E)
2646  {
2647  Synchronize Synchronize(this,sync);
2648  synchronize_lock_.unlock();
2649  Synchronize.run();
2650  synchronize_lock_.lock();
2651  }
2652  else if (sync & Mesh::EDGES_E)
2653  {
2654  mask_type tosync = Mesh::EDGES_E;
2655  Synchronize syncclass(this,tosync);
2656  boost::thread syncthread(syncclass);
2657  }
2658 
2659  if (sync == Mesh::NORMALS_E)
2660  {
2661  Synchronize Synchronize(this,sync);
2662  synchronize_lock_.unlock();
2663  Synchronize.run();
2664  synchronize_lock_.lock();
2665  }
2666  else if (sync & Mesh::NORMALS_E)
2667  {
2668  mask_type tosync = Mesh::NORMALS_E;
2669  Synchronize syncclass(this,tosync);
2670  boost::thread syncthread(syncclass);
2671  }
2672 
2673  if (sync == Mesh::NODE_NEIGHBORS_E)
2674  {
2675  Synchronize Synchronize(this,sync);
2676  synchronize_lock_.unlock();
2677  Synchronize.run();
2678  synchronize_lock_.lock();
2679  }
2680  else if (sync & Mesh::NODE_NEIGHBORS_E)
2681  {
2682  mask_type tosync = Mesh::NODE_NEIGHBORS_E;
2683  Synchronize syncclass(this,tosync);
2684  boost::thread syncthread(syncclass);
2685  }
2686 
2687  if (sync == Mesh::BOUNDING_BOX_E)
2688  {
2689  Synchronize Synchronize(this,sync);
2690  synchronize_lock_.unlock();
2691  Synchronize.run();
2692  synchronize_lock_.lock();
2693  }
2694  else if (sync & Mesh::BOUNDING_BOX_E)
2695  {
2696  mask_type tosync = Mesh::BOUNDING_BOX_E;
2697  Synchronize syncclass(this,tosync);
2698  boost::thread syncthread(syncclass);
2699  }
2700 
2701  if (sync == Mesh::NODE_LOCATE_E)
2702  {
2703  Synchronize Synchronize(this,sync);
2704  synchronize_lock_.unlock();
2705  Synchronize.run();
2706  synchronize_lock_.lock();
2707  }
2708  else if (sync & Mesh::NODE_LOCATE_E)
2709  {
2710  mask_type tosync = Mesh::NODE_LOCATE_E;
2711  Synchronize syncclass(this,tosync);
2712  boost::thread syncthread(syncclass);
2713  }
2714 
2715  if (sync == Mesh::ELEM_LOCATE_E)
2716  {
2717  Synchronize Synchronize(this,sync);
2718  synchronize_lock_.unlock();
2719  Synchronize.run();
2720  synchronize_lock_.lock();
2721  }
2722  else if (sync & Mesh::ELEM_LOCATE_E)
2723  {
2724  mask_type tosync = Mesh::ELEM_LOCATE_E;
2725  Synchronize syncclass(this,tosync);
2726  boost::thread syncthread(syncclass);
2727  }
2728 
2729  // Wait until threads are done
2730  while ((synchronized_ & sync) != sync)
2731  {
2732  synchronize_cond_.wait(lock);
2733  }
2734 
2735  }
2736 
2737  return (true);
2738 }
2739 
2740 template <class Basis>
2741 bool
2742 QuadSurfMesh<Basis>::unsynchronize(mask_type /*sync*/)
2743 {
2744  return (true);
2745 }
2746 
2747 
2748 template <class Basis>
2749 bool
2750 QuadSurfMesh<Basis>::clear_synchronization()
2751 {
2752  synchronize_lock_.lock();
2753  // Undo marking the synchronization
2754  synchronized_ = Mesh::NODES_E | Mesh::ELEMS_E | Mesh::FACES_E;
2755 
2756  if (synchronizing_)
2757  ASSERTFAIL("clear_synchronization is in parallel to a synchronization call from another thread.");
2758  // Free memory where possible
2759 
2760  edges_.clear();
2761  normals_.clear();
2762  halfedge_to_edge_.clear();
2763  node_neighbors_.clear();
2764 
2765  node_grid_.reset();
2766  elem_grid_.reset();
2767 
2768  synchronize_lock_.unlock();
2769  return (true);
2770 }
2771 
2772 template <class Basis>
2773 void
2774 QuadSurfMesh<Basis>::compute_normals()
2775 {
2776  normals_.resize(points_.size()); // 1 per node
2777 
2778  // build table of faces that touch each node
2779  std::vector<std::vector<typename Face::index_type> > node_in_faces(points_.size());
2780  /// face normals (not normalized) so that magnitude is also the area.
2781  std::vector<Core::Geometry::Vector> face_normals(faces_.size());
2782  // Computing normal per face.
2783  typename Node::array_type nodes(4);
2784  typename Face::iterator iter, iter_end;
2785  begin(iter);
2786  end(iter_end);
2787  while (iter != iter_end)
2788  {
2789  get_nodes(nodes, *iter);
2790 
2791  Core::Geometry::Point p0, p1, p2, p3;
2792  get_point(p0, nodes[0]);
2793  get_point(p1, nodes[1]);
2794  get_point(p2, nodes[2]);
2795  get_point(p3, nodes[3]);
2796 
2797  // build table of faces that touch each node
2798  node_in_faces[nodes[0]].push_back(*iter);
2799  node_in_faces[nodes[1]].push_back(*iter);
2800  node_in_faces[nodes[2]].push_back(*iter);
2801  node_in_faces[nodes[3]].push_back(*iter);
2802 
2803  Core::Geometry::Vector v0 = p1 - p0;
2804  Core::Geometry::Vector v1 = p2 - p1;
2805  Core::Geometry::Vector n = Cross(v0, v1);
2806  face_normals[*iter] = n;
2807 
2808  ++iter;
2809  }
2810  //Averaging the normals.
2811  typename std::vector<std::vector<typename Face::index_type> >::iterator nif_iter =
2812  node_in_faces.begin();
2813  index_type i = 0;
2814  while (nif_iter != node_in_faces.end())
2815  {
2816  const std::vector<typename Face::index_type> &v = *nif_iter;
2817  typename std::vector<typename Face::index_type>::const_iterator fiter =
2818  v.begin();
2819  Core::Geometry::Vector ave(0.L,0.L,0.L);
2820  while(fiter != v.end())
2821  {
2822  ave += face_normals[*fiter];
2823  ++fiter;
2824  }
2825  ave.safe_normalize();
2826  normals_[i] = ave; ++i;
2827  ++nif_iter;
2828  }
2829 
2830  synchronize_lock_.lock();
2831  synchronized_ |= Mesh::NORMALS_E;
2832  synchronize_lock_.unlock();
2833 }
2834 
2835 
2836 template <class Basis>
2837 typename QuadSurfMesh<Basis>::Node::index_type
2838 QuadSurfMesh<Basis>::add_find_point(const Core::Geometry::Point &p, double err)
2839 {
2840  typename Node::index_type i;
2841  if (locate(i, p) && (points_[i] - p).length2() < err)
2842  {
2843  return i;
2844  }
2845  else
2846  {
2847  points_.push_back(p);
2848  return static_cast<typename Node::index_type>
2849  (static_cast<size_type>(points_.size()) - 1);
2850  }
2851 }
2852 
2853 
2854 template <class Basis>
2855 typename QuadSurfMesh<Basis>::Elem::index_type
2856 QuadSurfMesh<Basis>::add_quad(typename Node::index_type a,
2857  typename Node::index_type b,
2858  typename Node::index_type c,
2859  typename Node::index_type d)
2860 {
2861  ASSERTMSG(order_face_nodes(a,b,c,d), "add_quad: element that is being added is invalid");
2862  faces_.push_back(a);
2863  faces_.push_back(b);
2864  faces_.push_back(c);
2865  faces_.push_back(d);
2866  synchronize_lock_.lock();
2867  synchronized_ &= ~Mesh::ELEM_NEIGHBORS_E;
2868  synchronized_ &= ~Mesh::NODE_NEIGHBORS_E;
2869  synchronized_ &= ~Mesh::NORMALS_E;
2870  synchronize_lock_.unlock();
2871  return static_cast<typename Elem::index_type>((static_cast<index_type>(faces_.size()) - 1) >> 2);
2872 }
2873 
2874 
2875 #ifdef HAVE_HASH_MAP
2876 
2877 struct edgehash
2878 {
2879  size_t operator()(const std::pair<index_type, index_type> &a) const
2880  {
2881 #if defined(__ECC) || defined(_MSC_VER)
2882  hash_compare<int> hasher;
2883 #else
2884  hash<int> hasher;
2885 #endif
2886  return hasher(static_cast<index_type>(hasher(a.first)) + a.second);
2887  }
2888 #if defined(__ECC) || defined(_MSC_VER)
2889 
2890  static const size_t bucket_size = 4;
2891  static const size_t min_buckets = 8;
2892 
2893  bool operator()(const std::pair<index_type, index_type> &a, const std::pair<index_type, index_type> &b) const
2894  {
2895  return a.first < b.first || a.first == b.first && a.second < b.second;
2896  }
2897 #endif
2898 };
2899 
2900 struct edgecompare
2901 {
2902  bool operator()(const std::pair<index_type, index_type> &a, const std::pair<index_type, index_type> &b) const
2903  {
2904  return a.first == b.first && a.second == b.second;
2905  }
2906 };
2907 
2908 #else
2909 
2910 struct edgecompare
2911 {
2912  bool operator()(const std::pair<index_type, index_type> &a, const std::pair<index_type, index_type> &b) const
2913  {
2914  return a.first < b.first || a.first == b.first && a.second < b.second;
2915  }
2916 };
2917 
2918 #endif
2919 
2920 #ifdef HAVE_HASH_MAP
2921 
2922 #if defined(__ECC) || defined(_MSC_VER)
2923 typedef hash_map<std::pair<index_type, index_type>, index_type, edgehash> EdgeMapType;
2924 #else
2925 typedef hash_map<std::pair<index_type, index_type>, index_type, edgehash, edgecompare> EdgeMapType;
2926 #endif
2927 
2928 #else
2929 
2930 typedef std::map<std::pair<index_type, index_type>, index_type, edgecompare> EdgeMapType;
2931 
2932 #endif
2933 
2934 #ifdef HAVE_HASH_MAP
2935 
2936 #if defined(__ECC) || defined(_MSC_VER)
2937 typedef hash_map<std::pair<index_type, index_type>, std::vector<index_type>, edgehash> EdgeMapType2;
2938 #else
2939 typedef hash_map<std::pair<index_type, index_type>, std::vector<index_type>, edgehash, edgecompare> EdgeMapType2;
2940 #endif
2941 
2942 #else
2943 
2944 typedef std::map<std::pair<index_type, index_type>, std::vector<index_type>, edgecompare> EdgeMapType2;
2945 
2946 #endif
2947 
2948 template <class Basis>
2949 void
2950 QuadSurfMesh<Basis>::compute_edges()
2951 {
2952  EdgeMapType2 edge_map;
2953 
2954  size_type num_faces = static_cast<index_type>(faces_.size()>>2);
2955  for (index_type i=0; i < num_faces; i++)
2956  {
2957  index_type n0 = faces_[4*i];
2958  index_type n1 = faces_[4*i+1];
2959  index_type n2 = faces_[4*i+2];
2960  index_type n3 = faces_[4*i+3];
2961 
2962 
2963  if (n0 > n1) edge_map[std::pair<index_type, index_type>(n1,n0)].push_back(i<<2);
2964  else if (n1 > n0) edge_map[std::pair<index_type, index_type>(n0,n1)].push_back(i<<2);
2965 
2966  if (n1 > n2) edge_map[std::pair<index_type, index_type>(n2,n1)].push_back((i<<2)+1);
2967  else if (n2 > n1) edge_map[std::pair<index_type, index_type>(n1,n2)].push_back((i<<2)+1);
2968 
2969  if (n2 > n3) edge_map[std::pair<index_type, index_type>(n3,n2)].push_back((i<<2)+2);
2970  else if (n3 > n2) edge_map[std::pair<index_type, index_type>(n2,n3)].push_back((i<<2)+2);
2971 
2972  if (n3 > n0) edge_map[std::pair<index_type, index_type>(n0,n3)].push_back((i<<2)+3);
2973  else if (n0 > n3) edge_map[std::pair<index_type, index_type>(n3,n0)].push_back((i<<2)+3);
2974  }
2975 
2976  typename EdgeMapType2::iterator itr;
2977  edges_.clear();
2978  edges_.resize(edge_map.size());
2979  halfedge_to_edge_.resize(faces_.size(),-1);
2980 
2981  size_t k=0;
2982  for (itr = edge_map.begin(); itr != edge_map.end(); ++itr)
2983  {
2984  const std::vector<index_type >& hedges = (*itr).second;
2985  for (size_t j=0; j<hedges.size(); j++)
2986  {
2987  index_type h = hedges[j];
2988  edges_[k].push_back(h);
2989  halfedge_to_edge_[h] = k;
2990  }
2991  k++;
2992  }
2993 
2994  synchronize_lock_.lock();
2995  synchronized_ |= Mesh::EDGES_E;
2996  synchronize_lock_.unlock();
2997 }
2998 
2999 
3000 template <class Basis>
3001 void
3002 QuadSurfMesh<Basis>::compute_node_neighbors()
3003 {
3004  node_neighbors_.clear();
3005  node_neighbors_.resize(points_.size());
3006  index_type i, num_elems = static_cast<index_type>(faces_.size());
3007  for (i = 0; i < num_elems; i++)
3008  {
3009  node_neighbors_[faces_[i]].push_back(i/4);
3010  }
3011 
3012  synchronize_lock_.lock();
3013  synchronized_ |= Mesh::NODE_NEIGHBORS_E;
3014  synchronize_lock_.unlock();
3015 
3016 }
3017 
3018 
3019 template <class Basis>
3020 void
3021 QuadSurfMesh<Basis>::insert_elem_into_grid(typename Elem::index_type ci)
3022 {
3023  /// @todo: This can crash if you insert a new cell outside of the grid.
3024  // Need to recompute grid at that point.
3025  const index_type idx = ci*4;
3026  Core::Geometry::BBox box;
3027  box.extend(points_[faces_[idx]]);
3028  box.extend(points_[faces_[idx+1]]);
3029  box.extend(points_[faces_[idx+2]]);
3030  box.extend(points_[faces_[idx+3]]);
3031  box.extend(epsilon_);
3032  elem_grid_->insert(ci, box);
3033 }
3034 
3035 
3036 template <class Basis>
3037 void
3038 QuadSurfMesh<Basis>::remove_elem_from_grid(typename Elem::index_type ci)
3039 {
3040  const index_type idx = ci*4;
3041  Core::Geometry::BBox box;
3042  box.extend(points_[faces_[idx]]);
3043  box.extend(points_[faces_[idx+1]]);
3044  box.extend(points_[faces_[idx+2]]);
3045  box.extend(points_[faces_[idx+3]]);
3046  box.extend(epsilon_);
3047  elem_grid_->remove(ci, box);
3048 }
3049 
3050 
3051 template <class Basis>
3052 void
3053 QuadSurfMesh<Basis>::insert_node_into_grid(typename Node::index_type ni)
3054 {
3055  /// @todo: This can crash if you insert a new cell outside of the grid.
3056  // Need to recompute grid at that point.
3057  node_grid_->insert(ni,points_[ni]);
3058 }
3059 
3060 
3061 template <class Basis>
3062 void
3063 QuadSurfMesh<Basis>::remove_node_from_grid(typename Node::index_type ni)
3064 {
3065  node_grid_->remove(ni,points_[ni]);
3066 }
3067 
3068 
3069 template <class Basis>
3070 void
3071 QuadSurfMesh<Basis>::compute_node_grid()
3072 {
3073  if (bbox_.valid())
3074  {
3075  // Cubed root of number of elems to get a subdivision ballpark.
3076 
3077  typename Node::size_type esz; size(esz);
3078 
3079  const size_type s = 3*static_cast<size_type>
3080  ((ceil(pow(static_cast<double>(esz) , (1.0/3.0))))/2.0 + 1.0);
3081 
3082  Core::Geometry::Vector diag = bbox_.diagonal();
3083  double trace = (diag.x()+diag.y()+diag.z());
3084  size_type sx = static_cast<size_type>(ceil(0.5+diag.x()/trace*s));
3085  size_type sy = static_cast<size_type>(ceil(0.5+diag.y()/trace*s));
3086  size_type sz = static_cast<size_type>(ceil(0.5+diag.z()/trace*s));
3087 
3088  Core::Geometry::BBox b = bbox_;
3089  b.extend(10*epsilon_);
3090  node_grid_.reset(new SearchGridT<index_type>(sx, sy, sz, b.min(), b.max()));
3091 
3092  typename Node::iterator ni, nie;
3093  begin(ni); end(nie);
3094  while(ni != nie)
3095  {
3096  insert_node_into_grid(*ni);
3097  ++ni;
3098  }
3099  }
3100 
3101  synchronize_lock_.lock();
3102  synchronized_ |= Mesh::NODE_LOCATE_E;
3103  synchronize_lock_.unlock();
3104 }
3105 
3106 template <class Basis>
3107 void
3108 QuadSurfMesh<Basis>::compute_elem_grid()
3109 {
3110  if (bbox_.valid())
3111  {
3112  // Cubed root of number of elems to get a subdivision ballpark.
3113 
3114  typename Elem::size_type esz; size(esz);
3115 
3116  const size_type s = 3*static_cast<size_type>
3117  ((ceil(pow(static_cast<double>(esz) , (1.0/3.0))))/2.0 + 1.0);
3118 
3119  Core::Geometry::Vector diag = bbox_.diagonal();
3120  double trace = (diag.x()+diag.y()+diag.z());
3121  size_type sx = static_cast<size_type>(ceil(0.5+diag.x()/trace*s));
3122  size_type sy = static_cast<size_type>(ceil(0.5+diag.y()/trace*s));
3123  size_type sz = static_cast<size_type>(ceil(0.5+diag.z()/trace*s));
3124 
3125  Core::Geometry::BBox b = bbox_;
3126  b.extend(10*epsilon_);
3127  elem_grid_.reset(new SearchGridT<index_type>(sx, sy, sz, b.min(), b.max()));
3128 
3129  typename Elem::iterator ci, cie;
3130  begin(ci); end(cie);
3131  while(ci != cie)
3132  {
3133  insert_elem_into_grid(*ci);
3134  ++ci;
3135  }
3136  }
3137 
3138  synchronize_lock_.lock();
3139  synchronized_ |= Mesh::ELEM_LOCATE_E;
3140  synchronize_lock_.unlock();
3141 }
3142 
3143 
3144 template <class Basis>
3145 void
3146 QuadSurfMesh<Basis>::compute_bounding_box()
3147 {
3148  bbox_.reset();
3149 
3150  // Compute bounding box
3151  typename Node::iterator ni, nie;
3152  begin(ni);
3153  end(nie);
3154  while (ni != nie)
3155  {
3156  bbox_.extend(points_[*ni]);
3157  ++ni;
3158  }
3159 
3160  // Compute epsilons associated with the bounding box
3161  epsilon_ = bbox_.diagonal().length()*1e-8;
3162  epsilon2_ = epsilon_*epsilon_;
3163 
3164  synchronize_lock_.lock();
3165  synchronized_ |= Mesh::BOUNDING_BOX_E;
3166  synchronize_lock_.unlock();
3167 }
3168 
3169 template <class Basis>
3170 typename QuadSurfMesh<Basis>::Node::index_type
3171 QuadSurfMesh<Basis>::add_point(const Core::Geometry::Point &p)
3172 {
3173  points_.push_back(p);
3174  return static_cast<typename Node::index_type>(static_cast<index_type>(points_.size() - 1));
3175 }
3176 
3177 
3178 template <class Basis>
3179 typename QuadSurfMesh<Basis>::Elem::index_type
3180 QuadSurfMesh<Basis>::add_quad(const Core::Geometry::Point &p0, const Core::Geometry::Point &p1,
3181  const Core::Geometry::Point &p2, const Core::Geometry::Point &p3)
3182 {
3183  return add_quad(add_find_point(p0), add_find_point(p1),
3184  add_find_point(p2), add_find_point(p3));
3185 }
3186 
3187 
3188 #define QUADSURFMESH_VERSION 4
3189 template <class Basis>
3190 void
3191 QuadSurfMesh<Basis>::io(Piostream &stream)
3192 {
3193  const int version = stream.begin_class(type_name(-1), QUADSURFMESH_VERSION);
3194 
3195  Mesh::io(stream);
3196 
3197  Pio(stream, points_);
3198  Pio_index(stream, faces_);
3199 
3200  // In case the face is degenerate
3201  // move the degerenaracy to the end
3202  // this way the visualization works fine
3203  if (version != 1)
3204  {
3205  if (stream.reading())
3206  {
3207  for (size_t i=0; i < faces_.size(); i += 4)
3208  if(!( order_face_nodes(faces_[i],faces_[i+1],faces_[i+2],faces_[i+3])))
3209  std::cerr << "Detected an invalid quadrilateral face\n";
3210  }
3211  }
3212 
3213  if (version == 1)
3214  {
3215  // We no longer save out this table, we actually no longer compute it
3216  // Hence when saving we save an empty table, which may break the old
3217  // SCIRun version (a version that was not working anyway)
3218  // When loading an old file we just ignore the table once it is loaded.
3219  std::vector<under_type> dummy;
3220  Pio(stream,dummy);
3221  }
3222 
3223  if (version >= 3)
3224  {
3225  basis_.io(stream);
3226  }
3227 
3228  stream.end_class();
3229 
3230  if (stream.reading())
3231  {
3232  synchronized_ = Mesh::NODES_E | Mesh::FACES_E | Mesh::CELLS_E;
3233  }
3234 
3235  if (stream.reading())
3236  vmesh_.reset(CreateVQuadSurfMesh(this));
3237 }
3238 
3239 
3240 template <class Basis>
3241 void
3242 QuadSurfMesh<Basis>::size(typename QuadSurfMesh::Node::size_type &s) const
3243 {
3244  typename Node::iterator itr; end(itr);
3245  s = *itr;
3246 }
3247 
3248 
3249 template <class Basis>
3250 void
3251 QuadSurfMesh<Basis>::size(typename QuadSurfMesh::Edge::size_type &s) const
3252 {
3253  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
3254  "Must call synchronize EDGES_E on QuadSurfMesh first");
3255  s = edges_.size();
3256 }
3257 
3258 
3259 template <class Basis>
3260 void
3261 QuadSurfMesh<Basis>::size(typename QuadSurfMesh::Face::size_type &s) const
3262 {
3263  typename Face::iterator itr; end(itr);
3264  s = *itr;
3265 }
3266 
3267 
3268 template <class Basis>
3269 void
3270 QuadSurfMesh<Basis>::size(typename QuadSurfMesh::Cell::size_type &s) const
3271 {
3272  typename Cell::iterator itr; end(itr);
3273  s = *itr;
3274 }
3275 
3276 
3277 template <class Basis>
3278 const TypeDescription*
3279 get_type_description(QuadSurfMesh<Basis> *)
3280 {
3281  static TypeDescription *td = 0;
3282  if (!td)
3283  {
3284  const TypeDescription *sub = get_type_description((Basis*)0);
3285  TypeDescription::td_vec *subs = new TypeDescription::td_vec(1);
3286  (*subs)[0] = sub;
3287  td = new TypeDescription("QuadSurfMesh", subs,
3288  std::string(__FILE__),
3289  "SCIRun",
3290  TypeDescription::MESH_E);
3291  }
3292  return td;
3293 }
3294 
3295 
3296 template <class Basis>
3297 const TypeDescription*
3298 QuadSurfMesh<Basis>::get_type_description() const
3299 {
3300  return SCIRun::get_type_description((QuadSurfMesh<Basis> *)0);
3301 }
3302 
3303 
3304 template <class Basis>
3305 const TypeDescription*
3306 QuadSurfMesh<Basis>::node_type_description()
3307 {
3308  static TypeDescription *td = 0;
3309  if (!td)
3310  {
3311  const TypeDescription *me =
3312  SCIRun::get_type_description((QuadSurfMesh<Basis> *)0);
3313  td = new TypeDescription(me->get_name() + "::Node",
3314  std::string(__FILE__),
3315  "SCIRun",
3316  TypeDescription::MESH_E);
3317  }
3318  return td;
3319 }
3320 
3321 
3322 template <class Basis>
3323 const TypeDescription*
3324 QuadSurfMesh<Basis>::edge_type_description()
3325 {
3326  static TypeDescription *td = 0;
3327  if (!td)
3328  {
3329  const TypeDescription *me =
3330  SCIRun::get_type_description((QuadSurfMesh<Basis> *)0);
3331  td = new TypeDescription(me->get_name() + "::Edge",
3332  std::string(__FILE__),
3333  "SCIRun",
3334  TypeDescription::MESH_E);
3335  }
3336  return td;
3337 }
3338 
3339 
3340 template <class Basis>
3341 const TypeDescription*
3342 QuadSurfMesh<Basis>::face_type_description()
3343 {
3344  static TypeDescription *td = 0;
3345  if (!td)
3346  {
3347  const TypeDescription *me =
3348  SCIRun::get_type_description((QuadSurfMesh<Basis> *)0);
3349  td = new TypeDescription(me->get_name() + "::Face",
3350  std::string(__FILE__),
3351  "SCIRun",
3352  TypeDescription::MESH_E);
3353  }
3354  return td;
3355 }
3356 
3357 
3358 template <class Basis>
3359 const TypeDescription*
3360 QuadSurfMesh<Basis>::cell_type_description()
3361 {
3362  static TypeDescription *td = 0;
3363  if (!td)
3364  {
3365  const TypeDescription *me =
3366  SCIRun::get_type_description((QuadSurfMesh<Basis> *)0);
3367  td = new TypeDescription(me->get_name() + "::Cell",
3368  std::string(__FILE__),
3369  "SCIRun",
3370  TypeDescription::MESH_E);
3371  }
3372  return td;
3373 }
3374 
3375 
3376 
3377 
3378 } // namespace SCIRun
3379 
3380 #endif
SCIRun::QuadSurfMesh::get_faces
void get_faces(typename Face::array_type &array, typename Face::index_type idx) const
Definition: QuadSurfMesh.h:449
SCIRun::QuadSurfMesh::get_delems
void get_delems(typename DElem::array_type &array, typename Face::index_type idx) const
Definition: QuadSurfMesh.h:488
SCIRun::QuadSurfMesh::topology_geometry
virtual int topology_geometry() const
Definition: QuadSurfMesh.h:355
SCIRun::QuadSurfMesh::get_elems
void get_elems(typename Elem::array_type &array, typename Node::index_type idx) const
Definition: QuadSurfMesh.h:469
StackVector.h
Interface to statically allocated std::vector class.
SCIRun::VUnstructuredMesh
Definition: VUnstructuredMesh.h:41
MeshSupport.h
SCIRun::QuadSurfMesh::get_weights
int get_weights(const Core::Geometry::Point &, typename Edge::array_type &, double *)
Definition: QuadSurfMesh.h:601
SCIRun::QuadSurfMesh::epsilon_
double epsilon_
Definition: QuadSurfMesh.h:2244
SCIRun::QuadSurfMesh::basis_order
virtual int basis_order()
Definition: QuadSurfMesh.h:347
SCIRun::QuadSurfMesh::vmesh
virtual VMesh * vmesh()
Access point to virtual interface.
Definition: QuadSurfMesh.h:343
SCIRun::Piostream::reading
bool reading() const
Definition: Persistent.h:164
SCIRun::QuadSurfMesh::Edge::array_type
std::vector< index_type > array_type
Definition: QuadSurfMesh.h:131
SCIRun::QuadSurfMesh::ElemData::edge1_index
index_type edge1_index() const
Definition: QuadSurfMesh.h:204
SCIRun::QuadSurfMesh::get_size
double get_size(typename Edge::index_type idx) const
Definition: QuadSurfMesh.h:530
SCIRun::NodeIterator
Distinct type for node FieldIterator.
Definition: FieldIterator.h:89
SCIRun::QuadSurfMesh::get_center
void get_center(Core::Geometry::Point &result, typename Node::index_type idx) const
get the center point (in object space) of an element
Definition: QuadSurfMesh.h:518
SCIRun::QuadSurfMesh::getFacade
MeshFacadeHandle getFacade() const
Definition: QuadSurfMesh.h:337
SCIRun::QuadSurfMesh::begin
void begin(typename Node::iterator &) const
begin/end iterators
Definition: QuadSurfMesh.h:2454
Point.h
SCIRun::FieldRNG
Definition: FieldRNG.h:37
SCIRun::TypeDescription::get_name
std::string get_name(const std::string &type_sep_start="<", const std::string &type_sep_end="> ") const
Definition: TypeDescription.cc:135
SCIRun::FaceIterator
Distinct type for face Iterator.
Definition: FieldIterator.h:121
SCIRun::QuadSurfMesh::synchronized_
mask_type synchronized_
Definition: QuadSurfMesh.h:2237
SCIRun::QuadSurfMesh::ElemData::node2
const Core::Geometry::Point & node2() const
Definition: QuadSurfMesh.h:230
SCIRun::QuadSurfMesh::Node::iterator
NodeIterator< under_type > iterator
Definition: QuadSurfMesh.h:122
SCIRun::Core::Thread::ConditionVariable
Definition: ConditionVariable.h:45
SCIRun::QuadSurfMesh::remove_node_from_grid
void remove_node_from_grid(typename Node::index_type ci)
Definition: QuadSurfMesh.h:3063
Vector.h
SCIRun::Core::Geometry::Transform::project
Vector project(const Vector &p) const
Definition: Transform.cc:425
SCIRun::QuadSurfMesh::next
index_type next(index_type i)
Definition: QuadSurfMesh.h:2209
SCIRun::QuadSurfMesh::points_
std::vector< Core::Geometry::Point > points_
array with all the points
Definition: QuadSurfMesh.h:2213
SCIRun::QuadSurfMesh::Synchronize::Synchronize
Synchronize(QuadSurfMesh< Basis > *mesh, mask_type sync)
Definition: QuadSurfMesh.h:255
SCIRun::QuadSurfMesh::get_neighbors
void get_neighbors(typename Elem::array_type &array, typename Elem::index_type elem) const
Definition: QuadSurfMesh.h:577
SCIRun::QuadSurfMesh::Edge::index_type
EdgeIndex< under_type > index_type
Definition: QuadSurfMesh.h:128
SCIRun::Mesh::IRREGULAR
Definition: Mesh.h:64
SCIRun::Core::Geometry::BBox::max
Point max() const
Definition: BBox.h:195
SCIRun::QuadSurfMesh::edges_
std::vector< std::vector< index_type > > edges_
Definition: QuadSurfMesh.h:2220
SCIRun::QuadSurfMesh::find_closest_nodes
bool find_closest_nodes(ARRAY1 &distances, ARRAY2 &nodes, const Core::Geometry::Point &p, double maxdist) const
Definition: QuadSurfMesh.h:1027
SCIRun::QuadSurfMesh::get_edge_center
void get_edge_center(Core::Geometry::Point &result, INDEX idx) const
Definition: QuadSurfMesh.h:2091
ASSERTMSG
#define ASSERTMSG(condition, message)
Definition: Exception.h:113
SCIRun::QuadSurfMesh::get_center
void get_center(Core::Geometry::Point &result, typename Face::index_type idx) const
Definition: QuadSurfMesh.h:522
SCIRun::QuadSurfMesh::~QuadSurfMesh
virtual ~QuadSurfMesh()
Destructor.
Definition: QuadSurfMesh.h:2375
SCIRun::QuadSurfMesh::synchronize
virtual bool synchronize(mask_type mask)
Definition: QuadSurfMesh.h:2610
SCIRun::QuadSurfMesh::get_edges_from_node
void get_edges_from_node(ARRAY &array, INDEX idx) const
Definition: QuadSurfMesh.h:1719
SCIRun::QuadSurfMesh::prev
index_type prev(index_type i)
Definition: QuadSurfMesh.h:2210
SCIRun::QuadSurfMesh::get_normal
void get_normal(Core::Geometry::Vector &result, VECTOR &coords, INDEX1 eidx, INDEX2)
Get the normals at the outside of the element.
Definition: QuadSurfMesh.h:625
SCIRun::QuadSurfMesh::get_elem_neighbor
bool get_elem_neighbor(INDEX1 &neighbor, INDEX1 elem, INDEX2 delem) const
Definition: QuadSurfMesh.h:1756
SCIRun::QuadSurfMesh::compute_node_grid
void compute_node_grid()
Definition: QuadSurfMesh.h:3071
SCIRun::QuadSurfMesh::size_type
SCIRun::size_type size_type
Definition: QuadSurfMesh.h:113
SCIRun::QuadSurfMesh::Cell
Definition: QuadSurfMesh.h:141
SCIRun::QuadSurfMesh::get_neighbors
void get_neighbors(std::vector< typename Node::index_type > &array, typename Node::index_type node)
These are more general implementations.
Definition: QuadSurfMesh.h:570
SCIRun::QuadSurfMesh::index_type
SCIRun::index_type index_type
Definition: QuadSurfMesh.h:112
SCIRun::QuadSurfMesh::Synchronize::operator()
void operator()()
Definition: QuadSurfMesh.h:260
SCIRun::Piostream
Definition: Persistent.h:89
QUADSURFMESH_VERSION
#define QUADSURFMESH_VERSION
Definition: QuadSurfMesh.h:3188
SCIRun::CreateVQuadSurfMesh
VMesh * CreateVQuadSurfMesh(MESH *mesh)
Definition: QuadSurfMesh.h:77
SCIRun::Core::Geometry::Transform::pre_scale
void pre_scale(const Vector &)
Definition: Transform.cc:193
SCIRun::QuadSurfMesh::locate
bool locate(typename Edge::index_type &edge, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:585
SCIRun::QuadSurfMesh::synchronize_cond_
Core::Thread::ConditionVariable synchronize_cond_
Definition: QuadSurfMesh.h:2234
SCIRun::QuadSurfMesh::elem_type_description
static const TypeDescription * elem_type_description()
Definition: QuadSurfMesh.h:1562
SCIRun::QuadSurfMesh::get_nodes
void get_nodes(typename Node::array_type &array, typename Face::index_type idx) const
Definition: QuadSurfMesh.h:423
SCIRun::QuadSurfMesh::get_node_center
void get_node_center(Core::Geometry::Point &p, INDEX idx) const
Definition: QuadSurfMesh.h:2085
SCIRun::QuadSurfMesh::find_closest_node
bool find_closest_node(double &pdist, Core::Geometry::Point &result, INDEX &node, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:849
SCIRun::QuadSurfMesh::get_faces
void get_faces(typename Face::array_type &array, typename Node::index_type idx) const
Definition: QuadSurfMesh.h:443
SCIRun::QuadSurfMesh::get_volume
double get_volume(typename Cell::index_type) const
Definition: QuadSurfMesh.h:555
SCIRun::QuadSurfMesh::type_name
static const std::string type_name(int n=-1)
Core functionality for getting the name of a templated mesh class.
Definition: QuadSurfMesh.h:2303
SCIRun::QuadSurfMesh::mesh_maker
static MeshHandle mesh_maker()
This function returns a handle for the virtual interface.
Definition: QuadSurfMesh.h:1568
SCIRun::QuadSurfMesh::handle_type
boost::shared_ptr< QuadSurfMesh< Basis > > handle_type
Definition: QuadSurfMesh.h:116
SCIRun::QuadSurfMesh::get_faces
void get_faces(typename Face::array_type &array, typename Edge::index_type idx) const
Definition: QuadSurfMesh.h:446
SCIRun::QuadSurfMesh::faces_
std::vector< index_type > faces_
array with the four nodes that make up a face
Definition: QuadSurfMesh.h:2215
SCIRun::QuadSurfMesh::add_quad
Elem::index_type add_quad(typename Node::index_type a, typename Node::index_type b, typename Node::index_type c, typename Node::index_type d)
Definition: QuadSurfMesh.h:2856
SCIRun::QuadSurfMesh::get_neighbor
bool get_neighbor(typename Elem::index_type &neighbor, typename Elem::index_type elem, typename DElem::index_type delem) const
Get neighbors of an element or a node.
Definition: QuadSurfMesh.h:564
SCIRun::Persistent
Definition: Persistent.h:187
SCIRun::QuadSurfMesh::derivate
void derivate(const VECTOR1 &coords, INDEX idx, VECTOR2 &J) const
Definition: QuadSurfMesh.h:720
SCIRun::QuadSurfMesh::scaled_jacobian_metric
double scaled_jacobian_metric(INDEX idx) const
Definition: QuadSurfMesh.h:785
SCIRun::QuadSurfMesh::elem_grid_
boost::shared_ptr< SearchGridT< index_type > > elem_grid_
Lookup grid for nodes.
Definition: QuadSurfMesh.h:2230
SCIRun::QuadSurfMesh::Node::index_type
NodeIndex< under_type > index_type
Definition: QuadSurfMesh.h:121
Transform.h
SCIRun::Core::Geometry::BBox::extend
BBox & extend(const Point &p)
Expand the bounding box to include point p.
Definition: BBox.h:98
SCIRun::QuadSurfMesh::get_elems
void get_elems(typename Elem::array_type &array, typename Face::index_type idx) const
Definition: QuadSurfMesh.h:475
SCIRun::QuadSurfMesh::add_node
Node::index_type add_node(const Core::Geometry::Point &p)
Definition: QuadSurfMesh.h:674
SCIRun::QuadSurfMesh::jacobian_metric
double jacobian_metric(INDEX idx) const
Definition: QuadSurfMesh.h:815
SCIRun::QuadSurfMesh::DElem
Edge DElem
Definition: QuadSurfMesh.h:151
SCIRun::QuadSurfMesh::edge_type_description
static const TypeDescription * edge_type_description()
Definition: QuadSurfMesh.h:3324
SCIRun::TypeDescription::MESH_E
Definition: TypeDescription.h:50
SCIRun::Dot
T Dot(const ColumnMatrixGeneric< T > &a, const ColumnMatrixGeneric< T > &b)
Definition: ColumnMatrixFunctions.h:155
TDCSInputData::mesh
FieldHandle mesh()
Definition: BuildTDCSMatrixTests.cc:56
SCIRun::QuadSurfMesh::det_jacobian
double det_jacobian(const VECTOR &coords, INDEX idx) const
Definition: QuadSurfMesh.h:729
SCIRun::Mesh::ELEM_LOCATE_E
Definition: Mesh.h:84
SCIRun::Core::Geometry::Point
Definition: Point.h:49
SCIRun::QuadSurfMesh::ElemData::node1_index
index_type node1_index() const
Definition: QuadSurfMesh.h:182
SCIRun::QuadSurfMesh::Cell::array_type
std::vector< index_type > array_type
Definition: QuadSurfMesh.h:145
SCIRun::Core::Thread::Mutex::get
boost::mutex & get()
Definition: Mutex.h:49
SCIRun::TypeDescription
Definition: TypeDescription.h:45
SCISHARE
#define SCISHARE
Definition: share.h:39
SCIRun::TypeDescription::td_vec
std::vector< const TypeDescription * > td_vec
Definition: TypeDescription.h:56
SCIRun::QuadSurfMesh::edgecompare
Virtual function table.
Definition: QuadSurfMesh.h:2284
SCIRun::EdgeMapType
std::map< std::pair< index_type, index_type >, index_type, edgecompare > EdgeMapType
Definition: QuadSurfMesh.h:2930
SCIRun::Mesh
Definition: Mesh.h:45
SCIRun::QuadSurfMesh::interpolate
void interpolate(Core::Geometry::Point &pt, const VECTOR &coords, INDEX idx) const
Definition: QuadSurfMesh.h:710
SCIRun::Core::Geometry::BBox
Definition: BBox.h:46
SCIRun::QuadSurfMesh::ElemData::node0_index
index_type node0_index() const
Definition: QuadSurfMesh.h:177
SCIRun::MeshHandle
boost::shared_ptr< Mesh > MeshHandle
Definition: DatatypeFwd.h:67
SCIRun::Core::Thread::Mutex
Definition: Mutex.h:43
SCIRun::QuadSurfMesh::get_bounding_box
virtual Core::Geometry::BBox get_bounding_box() const
Get the bounding box of the field.
Definition: QuadSurfMesh.h:2383
SCIRun::QuadSurfMesh::get_length
double get_length(typename Edge::index_type idx) const
More specific names for get_size.
Definition: QuadSurfMesh.h:551
SCIRun::QuadSurfMesh::get_size
double get_size(typename Face::index_type idx) const
Definition: QuadSurfMesh.h:536
SCIRun::QuadSurfMesh::dimensionality
virtual int dimensionality() const
Topological dimension.
Definition: QuadSurfMesh.h:350
SCIRun::CellIterator
Distinct type for cell Iterator.
Definition: FieldIterator.h:134
SCIRun::QuadSurfMesh::compute_normals
void compute_normals()
Definition: QuadSurfMesh.h:2774
VMesh.h
SCIRun::mask_type
unsigned int mask_type
Definition: Types.h:45
SCIRun::QuadSurfMesh::Edge::iterator
EdgeIterator< under_type > iterator
Definition: QuadSurfMesh.h:129
SCIRun::EdgeMapType2
std::map< std::pair< index_type, index_type >, std::vector< index_type >, edgecompare > EdgeMapType2
Definition: QuadSurfMesh.h:2944
FieldRNG.h
SCIRun::QuadSurfMesh::get_basis
Basis & get_basis()
Get the basis class.
Definition: QuadSurfMesh.h:384
SCIRun::Core::Geometry::Vector::asPoint
Point & asPoint() const
Definition: Vector.h:457
SCIRun::QuadSurfMesh::find_closest_nodes
bool find_closest_nodes(ARRAY &nodes, const Core::Geometry::Point &p, double maxdist) const
Definition: QuadSurfMesh.h:964
SCIRun::QuadSurfMesh::has_normals
virtual bool has_normals() const
Has this mesh normals.
Definition: QuadSurfMesh.h:375
ASSERT
#define ASSERT(condition)
Definition: Assert.h:110
SCIRun::QuadSurfMesh::Face::array_type
std::vector< index_type > array_type
Definition: QuadSurfMesh.h:138
SCIRun::QuadSurfMesh::inscribed_circumscribed_radius_metric
double inscribed_circumscribed_radius_metric(INDEX idx) const
Definition: QuadSurfMesh.h:843
SCIRun::QuadSurfMesh::Elem
Face Elem
Definition: QuadSurfMesh.h:150
SCIRun::QuadSurfMesh::ElemData::node2_index
index_type node2_index() const
Definition: QuadSurfMesh.h:187
SCIRun::QuadSurfMesh::get_normal
void get_normal(Core::Geometry::Vector &n, typename Node::index_type i) const
Normals for visualizations.
Definition: QuadSurfMesh.h:616
SCIRun::QuadSurfMesh::ElemData
Definition: QuadSurfMesh.h:159
SCIRun::QuadSurfMesh::get_delems
void get_delems(typename DElem::array_type &array, typename Node::index_type idx) const
Definition: QuadSurfMesh.h:482
SCIRun::Piostream::begin_class
virtual int begin_class(const std::string &name, int current_version)
Definition: Persistent.cc:143
SCIRun::QuadSurfMesh::get_size
double get_size(typename Cell::index_type) const
Definition: QuadSurfMesh.h:547
SCIRun::QuadSurfMesh::halfedge_to_edge_
std::vector< index_type > halfedge_to_edge_
Definition: QuadSurfMesh.h:2223
SCIRun::QuadSurfMesh::face_type_description
static const TypeDescription * face_type_description()
Definition: QuadSurfMesh.h:3342
SCIRun::QuadSurfMesh::io
virtual void io(Piostream &)
Export this class using the old Pio system.
Definition: QuadSurfMesh.h:3191
SCIRun::Core::Geometry::Vector
Definition: Vector.h:63
SCIRun::QuadSurfMesh::synchronize_lock_
Core::Thread::Mutex synchronize_lock_
Lookup grid for elements.
Definition: QuadSurfMesh.h:2233
SCIRun::Core::Thread::UniqueLock
boost::unique_lock< boost::mutex > UniqueLock
Definition: ConditionVariable.h:43
SCIRun::DetMatrix3P
double DetMatrix3P(const VectorOfPoints &p)
Inline templated determinant of matrix.
Definition: Locate.h:120
SCIRun::find_type_name
const string find_type_name(float *)
Definition: TypeName.cc:63
SCIRun::Mesh::EDGES_E
Definition: Mesh.h:72
SCIRun::QuadSurfMesh::node_type_description
static const TypeDescription * node_type_description()
Definition: QuadSurfMesh.h:3306
SCIRun::QuadSurfMesh::locate_node
bool locate_node(INDEX &node, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:1855
SCIRun::QuadSurfMesh::bbox_
Core::Geometry::BBox bbox_
Basis for interpolation.
Definition: QuadSurfMesh.h:2243
SCIRun::QuadSurfMesh::to_index
void to_index(typename Face::index_type &index, index_type i) const
Definition: QuadSurfMesh.h:411
SCIRun::QuadSurfMesh::node_grid_
boost::shared_ptr< SearchGridT< index_type > > node_grid_
normalized per node
Definition: QuadSurfMesh.h:2229
SCIRun::QuadSurfMesh::get_edges
void get_edges(typename Edge::array_type &, typename Cell::index_type) const
Definition: QuadSurfMesh.h:439
SCIRun::DetMatrix3x3
T DetMatrix3x3(const T *p)
Definition: Locate.h:95
SCIRun::QuadSurfMesh::synchronizing_
mask_type synchronizing_
Definition: QuadSurfMesh.h:2239
SCIRun::QuadSurfMesh::Edge::size_type
EdgeIndex< under_type > size_type
Definition: QuadSurfMesh.h:130
SCIRun::QuadSurfMesh::find_closest_elem
bool find_closest_elem(double &pdist, Core::Geometry::Point &result, INDEX &elem, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:1411
SCIRun::QuadSurfMesh::dynamic_type_name
virtual std::string dynamic_type_name() const
Definition: QuadSurfMesh.h:1553
SCIRun::QuadSurfMesh::get_elems
void get_elems(typename Elem::array_type &array, typename Edge::index_type idx) const
Definition: QuadSurfMesh.h:472
SCIRun::QuadSurfMesh::get_faces
void get_faces(typename Face::array_type &, typename Cell::index_type) const
Definition: QuadSurfMesh.h:452
SCIRun::Mesh::BOUNDING_BOX_E
Definition: Mesh.h:87
SCIRun::VMeshShared
Definition: VMeshShared.h:40
SCIRun::QuadSurfMesh::get_nodes
void get_nodes(typename Node::array_type &array, typename Node::index_type idx) const
Get the child elements of the given index.
Definition: QuadSurfMesh.h:417
max
Definition: ParallelLinearAlgebraTests.cc:358
SCIRun::QuadSurfMesh::get_epsilon
double get_epsilon() const
Definition: QuadSurfMesh.h:1539
SCIRun::QuadSurfMesh::Edge
Definition: QuadSurfMesh.h:127
SCIRun::QuadSurfMesh::compute_edges
void compute_edges()
Definition: QuadSurfMesh.h:2950
SCIRun::QuadSurfMesh::get_center
void get_center(Core::Geometry::Point &result, typename Edge::index_type idx) const
Definition: QuadSurfMesh.h:520
SCIRun::SearchGridT
Definition: SearchGridT.h:47
SCIRun::QuadSurfMesh::locate
bool locate(typename Node::index_type &node, const Core::Geometry::Point &p) const
Locate a point in a mesh, find which is the closest node.
Definition: QuadSurfMesh.h:583
name
const char * name[]
Definition: BoostGraphExampleTests.cc:87
SCIRun::QuadSurfMesh::find_closest_node
bool find_closest_node(double &pdist, Core::Geometry::Point &result, INDEX &node, const Core::Geometry::Point &p, double maxdist) const
Definition: QuadSurfMesh.h:856
SCIRun::QuadSurfMesh::ElemData::edge0_index
index_type edge0_index() const
changed the indexing as we now use unique indices
Definition: QuadSurfMesh.h:199
SCIRun::Core::Geometry::est_closest_point_on_quad
void est_closest_point_on_quad(Point &result, const Point &orig, const Point &p0, const Point &p1, const Point &p2, const Point &p3, const double epsilon)
Definition: CompGeom.cc:226
SCIRun::QuadSurfMesh::get_nodes
void get_nodes(typename Node::array_type &array, typename Edge::index_type idx) const
Definition: QuadSurfMesh.h:420
SCIRun::Pio_index
void Pio_index(Piostream &stream, index_type *data, size_type size)
Definition: Persistent.cc:606
SCIRun::QuadSurfMesh::ElemData::edge3_index
index_type edge3_index() const
Definition: QuadSurfMesh.h:214
SCIRun::size_type
long long size_type
Definition: Types.h:40
SCIRun::QuadSurfMesh::insert_node_into_grid
void insert_node_into_grid(typename Node::index_type ci)
Definition: QuadSurfMesh.h:3053
SCIRun::StackVector
Definition: StackVector.h:50
SCIRun::QuadSurfMesh::get_cells
void get_cells(typename Cell::array_type &, typename Cell::index_type) const
Definition: QuadSurfMesh.h:465
SCIRun::QuadSurfMesh::get_weights
int get_weights(const Core::Geometry::Point &p, typename Node::array_type &l, double *w)
Definition: QuadSurfMesh.h:2551
SCIRun::QuadSurfMesh::find_closest_elem
bool find_closest_elem(double &pdist, Core::Geometry::Point &result, ARRAY &coords, INDEX &elem, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:1098
SCIRun::Core::Geometry::Cross
Vector Cross(const Vector &v1, const Vector &v2)
Definition: Vector.h:378
SCIRun::QuadSurfMesh::NodeNeighborMap
std::vector< std::vector< typename Elem::index_type > > NodeNeighborMap
Definition: QuadSurfMesh.h:2225
SCIRun::ScaledDetMatrix3P
double ScaledDetMatrix3P(const VectorOfPoints &p)
Inline templated determinant of matrix.
Definition: Locate.h:132
SCIRun::QuadSurfMesh::inverse_jacobian
double inverse_jacobian(const VECTOR &coords, INDEX idx, double *Ji) const
Definition: QuadSurfMesh.h:762
SCIRun::QuadSurfMesh::pwl_approx_edge
void pwl_approx_edge(std::vector< VECTOR > &coords, INDEX ci, unsigned int which_edge, unsigned int div_per_unit) const
Definition: QuadSurfMesh.h:498
SCIRun::Mesh::NODES_E
Definition: Mesh.h:71
SCIRun::QuadSurfMesh::compute_elem_grid
void compute_elem_grid()
Definition: QuadSurfMesh.h:3108
SCIRun::QuadSurfMesh::get_nodes
void get_nodes(typename Node::array_type &, typename Cell::index_type) const
Definition: QuadSurfMesh.h:426
QuadBicubicHmt.h
SCIRun::Core::Geometry::BBox::diagonal
Vector diagonal() const
Definition: BBox.h:198
SCIRun::QuadSurfMesh::get_cells
void get_cells(typename Cell::array_type &, typename Edge::index_type) const
Definition: QuadSurfMesh.h:459
SCIRun::QuadSurfMesh::to_index
void to_index(typename Node::index_type &index, index_type i) const
Definition: QuadSurfMesh.h:407
SCIRun::QuadSurfMesh::locate_elem
bool locate_elem(INDEX &elem, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:1976
SCIRun::Mesh::ELEM_NEIGHBORS_E
Definition: Mesh.h:81
SCIRun::QuadSurfMesh::ElemData::ElemData
ElemData(const QuadSurfMesh< Basis > &msh, const index_type ind)
Definition: QuadSurfMesh.h:164
Locate.h
SCIRun::Mesh::NORMALS_E
Definition: Mesh.h:79
SCIRun::QuadSurfMesh::get_elem_neighbors
void get_elem_neighbors(ARRAY &array, INDEX idx) const
Definition: QuadSurfMesh.h:1801
SCIRun::QuadSurfMesh::get_delems
void get_delems(typename DElem::array_type &array, typename Edge::index_type idx) const
Definition: QuadSurfMesh.h:485
SCIRun::QuadSurfMesh::pwl_approx_face
void pwl_approx_face(std::vector< std::vector< VECTOR > > &coords, INDEX ci, unsigned int which_face, unsigned int div_per_unit) const
Definition: QuadSurfMesh.h:509
SCIRun::Mesh::CELLS_E
Definition: Mesh.h:74
PersistentSTL.h
Persistent i/o for STL containers.
SCIRun::SearchGridT::iterator
std::vector< INDEX >::iterator iterator
Definition: SearchGridT.h:53
SCIRun::QuadSurfMesh::is_editable
virtual bool is_editable() const
Check whether mesh can be altered by adding nodes or elements.
Definition: QuadSurfMesh.h:369
SCIRun::PersistentTypeID::type
std::string type
Definition: Persistent.h:72
SCIRun::QuadSurfMesh::get_faces_from_node
void get_faces_from_node(ARRAY &array, INDEX idx) const
Definition: QuadSurfMesh.h:1686
SCIRun::QuadSurfMesh::get_nodes_from_face
void get_nodes_from_face(ARRAY &array, INDEX idx) const
Definition: QuadSurfMesh.h:1641
SCIRun::QuadSurfMesh::get_edges_from_face
void get_edges_from_face(ARRAY &array, INDEX idx) const
Definition: QuadSurfMesh.h:1658
CompGeom.h
SCIRun::Core::Geometry::BBox::min
Point min() const
Definition: BBox.h:192
SCIRun::FaceIndex
Distinct type for face index.
Definition: FieldIndex.h:90
SCIRun::Mesh::UNSTRUCTURED
Definition: Mesh.h:62
SCIRun::edgecompare
Definition: QuadSurfMesh.h:2910
SCIRun::QuadSurfMesh::get_weights
int get_weights(const Core::Geometry::Point &, typename Cell::array_type &, double *)
Definition: QuadSurfMesh.h:604
SCIRun::QuadSurfMesh::get_node_neighbors
void get_node_neighbors(ARRAY &array, INDEX idx) const
Definition: QuadSurfMesh.h:1824
SCIRun::QuadSurfMesh::to_index
void to_index(typename Edge::index_type &index, index_type i) const
Definition: QuadSurfMesh.h:409
SCIRun::Core::Geometry::Transform
Definition: Transform.h:53
SCIRun::QuadSurfMesh::Face::index_type
FaceIndex< under_type > index_type
Definition: QuadSurfMesh.h:135
SCIRun::QuadSurfMesh::Face::iterator
FaceIterator< under_type > iterator
Definition: QuadSurfMesh.h:136
SCIRun::Core::Geometry::Vector::x
void x(double)
Definition: Vector.h:175
SCIRun::QuadSurfMesh::find_closest_elems
bool find_closest_elems(double &pdist, Core::Geometry::Point &result, ARRAY &elems, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:1421
SCIRun::QuadSurfMesh::to_index
void to_index(typename Cell::index_type &, index_type) const
Definition: QuadSurfMesh.h:413
SCIRun::QuadSurfMesh::clear_synchronization
bool clear_synchronization()
Definition: QuadSurfMesh.h:2750
SCIRun::QuadSurfMesh::Face
Definition: QuadSurfMesh.h:134
SCIRun::QuadSurfMesh::node_reserve
void node_reserve(size_type s)
Definition: QuadSurfMesh.h:691
SCIRun::QuadSurfMesh::Face::size_type
FaceIndex< under_type > size_type
Definition: QuadSurfMesh.h:137
SCIRun::QuadSurfMesh::resize_nodes
void resize_nodes(size_type s)
Definition: QuadSurfMesh.h:693
SCIRun::QuadSurfMesh::maker
static Persistent * maker()
This function returns a maker for Pio.
Definition: QuadSurfMesh.h:1566
SCIRun::QuadSurfMesh::mask_type
SCIRun::mask_type mask_type
Definition: QuadSurfMesh.h:114
SCIRun::Mesh::FACES_E
Definition: Mesh.h:73
SCIRun::TypeNameGenerator::make_template_id
static const std::string make_template_id(const std::string &templateName, const std::string &templateParam)
Definition: TypeName.h:62
SCIRun::QuadSurfMesh::get_elem_neighbors
bool get_elem_neighbors(ARRAY &array, INDEX1 elem, INDEX2 delem) const
Definition: QuadSurfMesh.h:1777
SCIRun::Core::Geometry::BBox::inside
bool inside(const Point &p) const
Definition: BBox.h:205
SCIRun::QuadSurfMesh::get_edges
void get_edges(typename Edge::array_type &array, typename Face::index_type idx) const
Definition: QuadSurfMesh.h:436
readAllFields.v
v
Definition: readAllFields.py:42
SCIRun::QuadSurfMesh::ElemData::edge2_index
index_type edge2_index() const
Definition: QuadSurfMesh.h:209
SCIRun::QuadSurfMesh::get_center
void get_center(Core::Geometry::Point &, typename Cell::index_type) const
Definition: QuadSurfMesh.h:524
SCIRun::Core::Geometry::Transform::pre_translate
void pre_translate(const Vector &)
Definition: Transform.cc:278
SCIRun::QuadSurfMesh::Cell::size_type
CellIndex< under_type > size_type
Definition: QuadSurfMesh.h:144
SCIRun::Core::Geometry::Vector::normalize
double normalize()
Definition: Vector.h:437
SCIRun::QuadSurfMesh::get_cells
void get_cells(typename Cell::array_type &, typename Face::index_type) const
Definition: QuadSurfMesh.h:462
SCIRun::QuadSurfMesh::compute_node_neighbors
void compute_node_neighbors()
Definition: QuadSurfMesh.h:3002
share.h
SCIRun::QuadSurfMesh::get_elems
void get_elems(typename Face::array_type &, typename Cell::index_type) const
Definition: QuadSurfMesh.h:478
SCIRun::Core::Geometry::Vector::y
void y(double)
Definition: Vector.h:185
SCIRun::QuadSurfMesh::locate
bool locate(typename Elem::index_type &elem, std::vector< double > &coords, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:592
SCIRun::QuadSurfMesh::Cell::iterator
CellIterator< under_type > iterator
Definition: QuadSurfMesh.h:143
SCIRun::Core::Geometry::Vector::safe_normalize
double safe_normalize()
Definition: Vector.h:236
SCIRun::QuadSurfMesh::order_face_nodes
bool order_face_nodes(NODE &n1, NODE &n2, NODE &n3, NODE &n4) const
Definition: QuadSurfMesh.h:2140
SCIRun::QuadSurfMesh::Node::size_type
NodeIndex< under_type > size_type
Definition: QuadSurfMesh.h:123
SCIRun::Core::Geometry::distance_to_line2
double distance_to_line2(const Point &p, const Point &a, const Point &b, const double epsilon)
Definition: CompGeom.cc:53
SCIRun::EdgeIterator
Distinct type for edge Iterator.
Definition: FieldIterator.h:105
SCIRun::QuadSurfMesh::add_find_point
Node::index_type add_find_point(const Core::Geometry::Point &p, double err=1.0e-3)
Definition: QuadSurfMesh.h:2838
SCIRun::QuadSurfMesh::get_faces_from_edge
void get_faces_from_edge(ARRAY &array, INDEX idx) const
Definition: QuadSurfMesh.h:1698
SCIRun::Mesh::ELEMS_E
Definition: Mesh.h:75
ASSERTFAIL
#define ASSERTFAIL(string)
Definition: Assert.h:52
SCIRun::QuadSurfMesh::size
void size(typename Node::size_type &) const
Definition: QuadSurfMesh.h:3242
Mesh.h
SCIRun::QuadSurfMesh
Definition: QuadSurfMesh.h:71
SCIRun::Pio
void Pio(Piostream &stream, Array1< T > &array)
Definition: Array1.h:65
FieldIterator.h
Some convenient simple iterators for fields.
SCIRun::QuadSurfMesh::ElemData::node1
const Core::Geometry::Point & node1() const
Definition: QuadSurfMesh.h:225
SCIRun::QuadSurfMesh::basis_type
Basis basis_type
Definition: QuadSurfMesh.h:117
SCIRun::QuadSurfMesh::ElemData::node0
const Core::Geometry::Point & node0() const
Definition: QuadSurfMesh.h:220
Mutex.h
SCIRun::QuadSurfMesh::node_neighbors_
NodeNeighborMap node_neighbors_
Definition: QuadSurfMesh.h:2226
SCIRun::QuadSurfMesh::get_area
double get_area(typename Face::index_type idx) const
Definition: QuadSurfMesh.h:553
SCIRun::QuadSurfMesh::get_coords
bool get_coords(VECTOR &coords, const Core::Geometry::Point &p, INDEX idx) const
Definition: QuadSurfMesh.h:701
SCIRun::Core::Geometry::Vector::length
double length() const
Definition: Vector.h:365
SCIRun::Mesh::NODE_NEIGHBORS_E
Definition: Mesh.h:80
SCIRun::Piostream::end_class
virtual void end_class()
Definition: Persistent.cc:178
SCIRun::QuadSurfMesh::get_face_center
void get_face_center(Core::Geometry::Point &p, INDEX idx) const
Definition: QuadSurfMesh.h:2102
SCIRun::index_type
long long index_type
Definition: Types.h:39
SCIRun::QuadSurfMesh::insert_elem_into_grid
void insert_elem_into_grid(typename Elem::index_type ci)
Used to recompute data for individual cells.
Definition: QuadSurfMesh.h:3021
ConditionVariable.h
SCIRun::QuadSurfMesh::normals_
std::vector< Core::Geometry::Vector > normals_
Definition: QuadSurfMesh.h:2228
SCIRun::Mesh::NODE_LOCATE_E
Definition: Mesh.h:83
SCIRun::QuadSurfMesh::epsilon2_
double epsilon2_
epsilon for calculations 1e-8*diagonal bounding box
Definition: QuadSurfMesh.h:2245
SCIRun::QuadSurfMesh::get_edges
void get_edges(typename Edge::array_type &array, typename Node::index_type idx)
Definition: QuadSurfMesh.h:430
SCIRun::StackVector::resize
void resize(size_t size, const value_type &val=value_type())
Definition: StackVector.h:61
SCIRun::QuadSurfMesh::add_elem
Elem::index_type add_elem(ARRAY a)
Add a new element to the mesh.
Definition: QuadSurfMesh.h:679
SCIRun::Mesh::DELEMS_E
Definition: Mesh.h:76
SCIRun::QuadSurfMesh::find_closest_elem
bool find_closest_elem(double &pdist, Core::Geometry::Point &result, ARRAY &coords, INDEX &elem, const Core::Geometry::Point &p, double maxdist) const
Definition: QuadSurfMesh.h:1108
SCIRun::QuadSurfMesh::ElemData::index_type
QuadSurfMesh< Basis >::index_type index_type
Definition: QuadSurfMesh.h:162
SearchGridT.h
SCIRun::QuadSurfMesh::Synchronize
Definition: QuadSurfMesh.h:252
SCIRun::QuadSurfMesh::QuadSurfMesh
QuadSurfMesh()
Construct a new mesh.
Definition: QuadSurfMesh.h:2324
BBox.h
SCIRun::Core::Geometry::BBox::reset
void reset()
Definition: BBox.h:95
SCIRun::QuadSurfMesh::ElemData::node3
const Core::Geometry::Point & node3() const
Definition: QuadSurfMesh.h:235
SCIRun::QuadSurfMesh::edgecompare::operator()
bool operator()(const std::pair< index_type, index_type > &a, const std::pair< index_type, index_type > &b) const
Definition: QuadSurfMesh.h:2286
SCIRun::QuadSurfMesh::get_cells
void get_cells(typename Cell::array_type &, typename Node::index_type) const
Definition: QuadSurfMesh.h:456
SCIRun::PersistentTypeID
Definition: Persistent.h:64
SCIRun::QuadSurfMesh::elem_reserve
void elem_reserve(size_type s)
Definition: QuadSurfMesh.h:692
SCIRun::QuadSurfMesh::get_random_point
void get_random_point(Core::Geometry::Point &, typename Elem::index_type, FieldRNG &rng) const
Definition: QuadSurfMesh.h:2571
SCIRun::QuadSurfMesh::get_point
void get_point(Core::Geometry::Point &p, typename Node::index_type i) const
Access the nodes of the mesh.
Definition: QuadSurfMesh.h:608
SCIRun::QuadSurfMesh::vmesh_
boost::shared_ptr< VMesh > vmesh_
square of epsilon for squared distance comparisons
Definition: QuadSurfMesh.h:2248
SCIRun::QuadSurfMesh::get_type_description
virtual const TypeDescription * get_type_description() const
Definition: QuadSurfMesh.h:3298
SCIRun::QuadSurfMesh::set_nodes_by_elem
void set_nodes_by_elem(ARRAY &array, INDEX idx)
Definition: QuadSurfMesh.h:1747
SCIRun::QuadSurfMesh::under_type
SCIRun::index_type under_type
Definition: QuadSurfMesh.h:111
SCIRun::QuadSurfMesh::find_closest_elem_far_from
bool find_closest_elem_far_from(double &pdist, Core::Geometry::Point &result, ARRAY &coords, INDEX &face, const Core::Geometry::Point &p, const Core::Geometry::Point &repelpos, double mindistfromrepel, const Core::Geometry::Vector &repelnormal, double mindotfromrepel, double maxdist) const
Definition: QuadSurfMesh.h:1269
SCIRun::QuadSurfMesh::quadsurfmesh_typeid
static PersistentTypeID quadsurfmesh_typeid
This ID is created as soon as this class will be instantiated.
Definition: QuadSurfMesh.h:1549
SCIRun::QuadSurfMesh::resize_elems
void resize_elems(size_type s)
Definition: QuadSurfMesh.h:694
QuadBilinearLgn.h
SCIRun::EdgeIndex
Distinct type for edge index.
Definition: FieldIndex.h:81
SCIRun::QuadSurfMesh::cell_type_description
static const TypeDescription * cell_type_description()
Definition: QuadSurfMesh.h:3360
SCIRun::QuadSurfMesh::VQuadSurfMesh
friend class VQuadSurfMesh
Make sure the virtual interface has access.
Definition: QuadSurfMesh.h:105
eab.n
int n
Definition: eab.py:9
SCIRun::QuadSurfMesh::get_size
double get_size(typename Node::index_type) const
Get the size of an element (length, area, volume)
Definition: QuadSurfMesh.h:528
SCIRun::Core::Thread::Guard
boost::lock_guard< boost::mutex > Guard
Definition: Mutex.h:55
SCIRun::QuadSurfMesh::Cell::index_type
CellIndex< under_type > index_type
Definition: QuadSurfMesh.h:142
SCIRun::InverseMatrix3x3
T InverseMatrix3x3(const T *p, T *q)
Definition: Locate.h:47
SCIRun::QuadSurfMesh::Node
Index and Iterator types required for Mesh Concept.
Definition: QuadSurfMesh.h:120
SCIRun::QuadSurfMesh::Node::array_type
StackVector< index_type, 4 > array_type
Definition: QuadSurfMesh.h:124
SCIRun::QuadSurfMesh::Synchronize::run
void run()
Definition: QuadSurfMesh.h:265
DEBUG_CONSTRUCTOR
#define DEBUG_CONSTRUCTOR(type)
Definition: Debug.h:64
SCIRun::Mesh::io
void io(Piostream &stream)
Persistent I/O.
Definition: Mesh.cc:387
SCIRun::QuadSurfMesh::locate_elems
bool locate_elems(ARRAY &array, const Core::Geometry::BBox &b) const
Definition: QuadSurfMesh.h:2011
SCIRun::QuadSurfMesh::get_edges_from_elem
void get_edges_from_elem(ARRAY &array, INDEX idx) const
Definition: QuadSurfMesh.h:1679
SCIRun::QuadSurfMesh::jacobian
void jacobian(const VECTOR &coords, INDEX idx, double *J) const
Definition: QuadSurfMesh.h:740
SCIRun::QuadSurfMesh::remove_elem_from_grid
void remove_elem_from_grid(typename Elem::index_type ci)
Definition: QuadSurfMesh.h:3038
SCIRun::QuadSurfMesh::get_neighbors
bool get_neighbors(std::vector< typename Elem::index_type > &array, typename Elem::index_type elem, typename DElem::index_type delem) const
Definition: QuadSurfMesh.h:573
SCIRun::QuadSurfMesh::locate_elem
bool locate_elem(INDEX &elem, ARRAY &coords, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:2042
SCIRun::QuadSurfMesh::transform
virtual void transform(const Core::Geometry::Transform &t)
Core::Geometry::Transform a field (transform all nodes using this transformation matrix) ...
Definition: QuadSurfMesh.h:2413
SCIRun::Core::Geometry::Vector::z
void z(double)
Definition: Vector.h:195
SCIRun::QuadSurfMesh::add_point
Node::index_type add_point(const Core::Geometry::Point &p)
Add a new node to the mesh.
Definition: QuadSurfMesh.h:3171
SCIRun::Core::Datatypes::MeshTraits< VMesh >::MeshFacadeHandle
boost::shared_ptr< MeshFacade< VMesh > > MeshFacadeHandle
Definition: MeshTraits.h:61
SCIRun::VMesh
Definition: VMesh.h:53
SCIRun::QuadSurfMesh::inside
bool inside(INDEX idx, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:1575
SCIRun::QuadSurfMesh::ElemData::node3_index
index_type node3_index() const
Definition: QuadSurfMesh.h:192
SCIRun::QuadSurfMesh::locate_edge
bool locate_edge(INDEX &loc, const Core::Geometry::Point &p) const
This is currently implemented as an exhaustive search.
Definition: QuadSurfMesh.h:1945
SCIRun::QuadSurfMesh::get_delems
void get_delems(typename DElem::array_type &, typename Cell::index_type) const
Definition: QuadSurfMesh.h:491
SCIRun::edgecompare::operator()
bool operator()(const std::pair< index_type, index_type > &a, const std::pair< index_type, index_type > &b) const
Definition: QuadSurfMesh.h:2912
SCIRun::QuadSurfMesh::unsynchronize
virtual bool unsynchronize(mask_type mask)
Definition: QuadSurfMesh.h:2742
VirtualMeshFacade.h
SCIRun::QuadSurfMesh::basis_
Basis basis_
Definition: QuadSurfMesh.h:2241
SCIRun::QuadSurfMesh::clone
virtual QuadSurfMesh * clone() const
Definition: QuadSurfMesh.h:332
SCIRun::QuadSurfMesh::get_nodes_from_elem
void get_nodes_from_elem(ARRAY &array, INDEX idx) const
Definition: QuadSurfMesh.h:1652
SCIRun::QuadSurfMesh::point
const Core::Geometry::Point & point(typename Node::index_type i) const
Definition: QuadSurfMesh.h:2121
SCIRun::QuadSurfMesh::get_canonical_transform
virtual void get_canonical_transform(Core::Geometry::Transform &t) const
Definition: QuadSurfMesh.h:2403
SCIRun::QuadSurfMesh::locate
bool locate(typename Face::index_type &face, const Core::Geometry::Point &p) const
Definition: QuadSurfMesh.h:587
SCIRun::QuadSurfMesh::locate
bool locate(typename Cell::index_type &, const Core::Geometry::Point &) const
Definition: QuadSurfMesh.h:589
SCIRun::QuadSurfMesh::set_point
void set_point(const Core::Geometry::Point &p, typename Node::index_type i)
Definition: QuadSurfMesh.h:610
SCIRun::QuadSurfMesh::get_edges
void get_edges(typename Edge::array_type &array, typename Edge::index_type idx) const
Definition: QuadSurfMesh.h:433
SCIRun::QuadSurfMesh::get_nodes_from_edge
void get_nodes_from_edge(ARRAY &array, INDEX i) const
Definition: QuadSurfMesh.h:1624
eabLatVolData.size
int size
Definition: eabLatVolData.py:2
SCIRun::QuadSurfMesh::end
void end(typename Node::iterator &) const
Definition: QuadSurfMesh.h:2464
QuadBiquadraticLgn.h
DEBUG_DESTRUCTOR
#define DEBUG_DESTRUCTOR(type)
Definition: Debug.h:65
SCIRun::get_type_description
const TypeDescription * get_type_description(Core::Basis::ConstantBasis< T > *)
Definition: Constant.h:209
SCIRun::QuadSurfMesh::compute_bounding_box
void compute_bounding_box()
Definition: QuadSurfMesh.h:3146
SCIRun::elem_locate
bool elem_locate(INDEX &elem, MESH &msh, const Core::Geometry::Point &p)
General case locate, search each elem.
Definition: Mesh.h:188
Generated on Wed Jul 9 2014 14:31:45 for SCIRun by   doxygen 1.8.5