SCIRun  5.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TetVolMesh.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 
30 #ifndef CORE_DATATYPES_TETVOLMESH_H
31 #define CORE_DATATYPES_TETVOLMESH_H 1
32 
33 /// Include what kind of support we want to have
34 /// Need to fix this and couple it sci-defs
35 #include <Core/Datatypes/Legacy/Field/MeshSupport.h>
36 
37 #include <Core/Containers/StackVector.h>
38 #include <Core/Persistent/PersistentSTL.h>
39 
40 #include <Core/GeometryPrimitives/SearchGridT.h>
41 #include <Core/GeometryPrimitives/BBox.h>
42 #include <Core/GeometryPrimitives/CompGeom.h>
43 #include <Core/GeometryPrimitives/Point.h>
44 #include <Core/GeometryPrimitives/Plane.h>
45 #include <Core/GeometryPrimitives/Transform.h>
46 #include <Core/GeometryPrimitives/Vector.h>
47 
48 #include <Core/Basis/Locate.h>
49 #include <Core/Basis/TetLinearLgn.h>
50 #include <Core/Basis/TetQuadraticLgn.h>
51 #include <Core/Basis/TetCubicHmt.h>
52 
53 #include <Core/Datatypes/Legacy/Field/FieldIterator.h>
54 #include <Core/Datatypes/Legacy/Field/FieldRNG.h>
55 #include <Core/Datatypes/Legacy/Field/Mesh.h>
56 #include <Core/Datatypes/Legacy/Field/VMesh.h>
57 #include <Core/Datatypes/Mesh/VirtualMeshFacade.h>
58 #include <Core/Math/MiscMath.h>
59 
60 #include <Core/Utils/Legacy/CheckSum.h>
61 
62 #include <boost/thread.hpp>
63 #include <Core/Thread/Mutex.h>
64 #include <Core/Thread/ConditionVariable.h>
65 
66 #include <set>
67 
68 #include <Core/Datatypes/Legacy/Field/share.h>
69 
70 namespace SCIRun {
71 
72 /////////////////////////////////////////////////////
73 // Declarations for virtual interface
74 
75 /// Functions for creating the virtual interface
76 /// Declare the functions that instantiate the virtual interface
77 template <class Basis> class TetVolMesh;
78 
79 /// make sure any other mesh other than the preinstantiate ones
80 /// returns no virtual interface. Altering this behavior will allow
81 /// for dynamically compiling the interface if needed.
82 template<class MESH>
83 VMesh* CreateVTetVolMesh(MESH*) { return (0); }
84 
85 /// These declarations are needed for a combined dynamic compilation as
86 /// as well as virtual functions solution.
87 /// Declare that these can be found in a library that is already
88 /// precompiled. So dynamic compilation will not instantiate them again.
89 
90 #if (SCIRUN_TETVOL_SUPPORT > 0)
91 
92 SCISHARE VMesh* CreateVTetVolMesh(TetVolMesh<Core::Basis::TetLinearLgn<Core::Geometry::Point> >* mesh);
93 #if (SCIRUN_QUADRATIC_SUPPORT > 0)
94 SCISHARE VMesh* CreateVTetVolMesh(TetVolMesh<Core::Basis::TetQuadraticLgn<Core::Geometry::Point> >* mesh);
95 #endif
96 #if (SCIRUN_CUBIC_SUPPORT > 0)
97 SCISHARE VMesh* CreateVTetVolMesh(TetVolMesh<Core::Basis::TetCubicHmt<Core::Geometry::Point> >* mesh);
98 #endif
99 
100 #endif
101 /////////////////////////////////////////////////////
102 
103 static const int
104 TetVolEdgePerNodeTable[4][6] = {{0,1, 2,0, 3,0},
105  {0,1, 1,2, 3,1},
106  {1,2, 2,0, 3,2},
107  {3,0, 3,1, 3,2}};
108 
109 static const int
110 TetVolFacePerNodeTable[4][9] = {{0,1,2, 0,1,3, 0,3,2},
111  {0,1,2, 1,2,3, 0,1,3},
112  {0,1,2, 1,2,3, 0,3,2},
113  {1,2,3, 0,1,3, 0,3,2}};
114 
115 static const int
116 TetVolFacePerEdgeTable[6][6] = {{0,1,2, 0,1,3}, {0,1,2, 1,2,3},
117  {0,1,2, 0,3,2}, {0,1,3, 0,3,2},
118  {1,2,3, 0,1,3}, {1,2,3, 0,3,2}};
119 
120 static const int
121 TetVolEdgePerFaceTable[4][6] = {{0,1, 1,2, 2,0 },
122  {1,2, 3,1, 3,2 },
123  {0,1, 3,0, 3,1 },
124  {2,0, 3,0, 3,2}};
125 
126 static const int
127 TetVolEdgeTable[6][2] = {{0,1},{1,2},{2,0},{3,0},{3,1},{3,2}};
128 
129 static const int
130 TetVolFaceTable[4][3] = { {0,2,1},{1,2,3},{0,1,3},{0,3,2}};
131 
132 /////////////////////////////////////////////////////
133 // Declarations for TetVolMesh class
134 
135 template <class Basis>
136 class TetVolMesh : public Mesh
137 {
138  /// Make sure the virtual interface has access
139  template<class MESH> friend class VTetVolMesh;
140  template<class MESH> friend class VMeshShared;
141  template<class MESH> friend class VUnstructuredMesh;
142 
143 public:
144  // Types that change depending on 32 or 64bits
145  typedef SCIRun::index_type under_type;
146  typedef SCIRun::index_type index_type;
147  typedef SCIRun::size_type size_type;
148  typedef SCIRun::mask_type mask_type;
149 
150  typedef boost::shared_ptr<TetVolMesh<Basis> > handle_type;
151  typedef Basis basis_type;
152 
153  /// Index and Iterator types required for Mesh Concept.
154  struct Node {
155  typedef NodeIndex<under_type> index_type;
156  typedef NodeIterator<under_type> iterator;
157  typedef NodeIndex<under_type> size_type;
158  typedef StackVector<index_type, 4> array_type;
159  };
160 
161  struct Edge {
162  typedef EdgeIndex<under_type> index_type;
163  typedef EdgeIterator<under_type> iterator;
164  typedef EdgeIndex<under_type> size_type;
165  typedef std::vector<index_type> array_type;
166  };
167 
168  struct Face {
169  typedef FaceIndex<under_type> index_type;
170  typedef FaceIterator<under_type> iterator;
171  typedef FaceIndex<under_type> size_type;
172  typedef std::vector<index_type> array_type;
173  };
174 
175  struct Cell {
176  typedef CellIndex<under_type> index_type;
177  typedef CellIterator<under_type> iterator;
178  typedef CellIndex<under_type> size_type;
179  typedef std::vector<index_type> array_type;
180  };
181 
182  /// Elem refers to the most complex topological object
183  /// DElem refers to object just below Elem in the topological hierarchy
184 
185  typedef Cell Elem;
186  typedef Face DElem;
187 
188  /// Somehow the information of how to interpolate inside an element
189  /// ended up in a separate class, as they need to share information
190  /// this construction was created to transfer data.
191  /// Hopefully in the future this class will disappear again.
192  friend class ElemData;
193 
194  class ElemData
195  {
196  public:
197  typedef typename TetVolMesh<Basis>::index_type index_type;
198 
199  ElemData(const TetVolMesh<Basis>& msh, const index_type ind) :
200  mesh_(msh),
201  index_(ind)
202  {
203  /// Linear and Constant Basis never use edges_
204  if (basis_type::polynomial_order() > 1)
205  {
206  mesh_.get_edges_from_cell(edges_, index_);
207  }
208  }
209 
210  /// the following designed to coordinate with ::get_nodes
211  inline
212  index_type node0_index() const
213  {
214  return mesh_.cells_[index_ * 4];
215  }
216  inline
217  index_type node1_index() const
218  {
219  return mesh_.cells_[index_ * 4 + 1];
220  }
221  inline
222  index_type node2_index() const
223  {
224  return mesh_.cells_[index_ * 4 + 2];
225  }
226  inline
227  index_type node3_index() const
228  {
229  return mesh_.cells_[index_ * 4 + 3];
230  }
231 
232  /// the following designed to coordinate with ::get_edges
233  inline
234  index_type edge0_index() const
235  {
236  return edges_[0];
237  }
238  inline
239  index_type edge1_index() const
240  {
241  return edges_[1];
242  }
243  inline
244  index_type edge2_index() const
245  {
246  return edges_[2];
247  }
248  inline
249  index_type edge3_index() const
250  {
251  return edges_[3];
252  }
253  inline
254  index_type edge4_index() const
255  {
256  return edges_[4];
257  }
258  inline
259  index_type edge5_index() const
260  {
261  return edges_[5];
262  }
263 
264  inline
265  const Core::Geometry::Point &node0() const
266  {
267  return mesh_.points_[node0_index()];
268  }
269  inline
270  const Core::Geometry::Point &node1() const
271  {
272  return mesh_.points_[node1_index()];
273  }
274  inline
275  const Core::Geometry::Point &node2() const
276  {
277  return mesh_.points_[node2_index()];
278  }
279  inline
280  const Core::Geometry::Point &node3() const
281  {
282  return mesh_.points_[node3_index()];
283  }
284 
285  private:
286  /// reference to the mesh
287  const TetVolMesh<Basis> &mesh_;
288  /// copy of element index
289  const index_type index_;
290  /// need edges for quadratic meshes
291  typename Edge::array_type edges_;
292  };
293 
294  friend class Synchronize;
295 
296  class Synchronize //: public Runnable
297  {
298  public:
299  Synchronize(TetVolMesh<Basis>* mesh, mask_type sync) :
300  mesh_(mesh), sync_(sync) {}
301 
302  void operator()()
303  {
304  run();
305  }
306 
307  void run()
308  {
309  mesh_->synchronize_lock_.lock();
310  // block out all the tables that are already synchronized
311  sync_ &= ~(mesh_->synchronized_);
312  // block out all the tables that are already being computed
313  sync_ &= ~(mesh_->synchronizing_);
314  // Now sync_ contains what this thread will synchronize
315  // Denote what this thread will synchronize
316  mesh_->synchronizing_ |= sync_;
317  // Other threads now know what this tread will be doing
318  mesh_->synchronize_lock_.unlock();
319 
320  if (sync_ & Mesh::NODE_NEIGHBORS_E) mesh_->compute_node_neighbors();
321  if (sync_ & Mesh::EDGES_E) mesh_->compute_edges();
322  if (sync_ & Mesh::FACES_E) mesh_->compute_faces();
323  if (sync_ & Mesh::BOUNDING_BOX_E) mesh_->compute_bounding_box();
324 
325  // These depend on the bounding box being synchronized
326  if (sync_ & (Mesh::NODE_LOCATE_E|Mesh::ELEM_LOCATE_E))
327  {
328  {
329  Core::Thread::UniqueLock lock(mesh_->synchronize_lock_.get());
330  while(!(mesh_->synchronized_ & Mesh::BOUNDING_BOX_E))
331  mesh_->synchronize_cond_.wait(lock);
332  }
333  if (sync_ & Mesh::NODE_LOCATE_E) mesh_->compute_node_grid();
334  if (sync_ & Mesh::ELEM_LOCATE_E) mesh_->compute_elem_grid();
335  }
336 
337  mesh_->synchronize_lock_.lock();
338  // Mark the ones that were just synchronized
339  mesh_->synchronized_ |= sync_;
340  // Unmark the the ones that were done
341  mesh_->synchronizing_ &= ~(sync_);
342  /// Tell other threads we are done
343  mesh_->synchronize_cond_.conditionBroadcast();
344  mesh_->synchronize_lock_.unlock();
345  }
346 
347  private:
348  TetVolMesh<Basis>* mesh_;
349  mask_type sync_;
350  };
351 
352  //////////////////////////////////////////////////////////////////
353 
354  /// Construct a new mesh
355  TetVolMesh();
356 
357  /// Copy a mesh, needed for detaching the mesh from a field
358  TetVolMesh(const TetVolMesh &copy);
359 
360  /// Clone function for detaching the mesh and automatically generating
361  /// a new version if needed.
362  virtual TetVolMesh *clone() const { return new TetVolMesh(*this); }
363 
364  /// Destructor
365  virtual ~TetVolMesh();
366 
367  /// Access point to virtual interface
368  virtual VMesh* vmesh() { return vmesh_.get(); }
369 
370  MeshFacadeHandle getFacade() const
371  {
372  return boost::make_shared<Core::Datatypes::VirtualMeshFacade<VMesh>>(vmesh_);
373  }
374 
375  /// This one should go at some point, should be reroute through the
376  /// virtual interface
377  virtual int basis_order() { return (basis_.polynomial_order()); }
378 
379  /// Topological dimension
380  virtual int dimensionality() const { return 3; }
381 
382  /// What kind of mesh is this
383  /// structured = no connectivity data
384  /// regular = no node location data
385  virtual int topology_geometry() const
386  { return (Mesh::UNSTRUCTURED | Mesh::IRREGULAR); }
387 
388  /// Get the bounding box of the field
389  virtual Core::Geometry::BBox get_bounding_box() const;
390 
391  /// Return the transformation that takes a 0-1 space bounding box
392  /// to the current bounding box of this mesh.
393  virtual void get_canonical_transform(Core::Geometry::Transform &t) const;
394 
395  /// Core::Geometry::Transform a field (transform all nodes using this transformation matrix)
396  virtual void transform(const Core::Geometry::Transform &t);
397 
398  /// Check whether mesh can be altered by adding nodes or elements
399  virtual bool is_editable() const { return (true); }
400 
401  /// Has this mesh normals.
402  virtual bool has_normals() const { return (false); }
403 
404  /// Has this mesh face normals
405  virtual bool has_face_normals() const { return (true); }
406 
407  double get_epsilon() const { return (epsilon_); }
408 
409  /// Compute tables for doing topology, these need to be synchronized
410  /// before doing a lot of operations.
411  virtual bool synchronize(mask_type mask);
412  virtual bool unsynchronize(mask_type mask);
413  bool clear_synchronization();
414 
415  /// Get the basis class.
416  Basis& get_basis() { return basis_; }
417 
418  /// begin/end iterators
419  void begin(typename Node::iterator &) const;
420  void begin(typename Edge::iterator &) const;
421  void begin(typename Face::iterator &) const;
422  void begin(typename Cell::iterator &) const;
423 
424  void end(typename Node::iterator &) const;
425  void end(typename Edge::iterator &) const;
426  void end(typename Face::iterator &) const;
427  void end(typename Cell::iterator &) const;
428 
429  /// Get the iteration sizes
430  void size(typename Node::size_type &) const;
431  void size(typename Edge::size_type &) const;
432  void size(typename Face::size_type &) const;
433  void size(typename Cell::size_type &) const;
434 
435  /// These are here to convert indices to unsigned int
436  /// counters. Some how the decision was made to use multi
437  /// dimensional indices in some fields, these functions
438  /// should deal with different pointer types.
439  /// Use the virtual interface to avoid all this non sense.
440  void to_index(typename Node::index_type &index, index_type i) const
441  { index = i; }
442  void to_index(typename Edge::index_type &index, index_type i) const
443  { index = i; }
444  void to_index(typename Face::index_type &index, index_type i) const
445  { index = i; }
446  void to_index(typename Cell::index_type &index, index_type i) const
447  { index = i; }
448 
449  /// Get the child topology elements of the given topology
450  void get_nodes(typename Node::array_type &array,
451  typename Node::index_type idx) const
452  { array.resize(1); array[0]= idx; }
453  void get_nodes(typename Node::array_type &array,
454  typename Edge::index_type idx) const
455  { get_nodes_from_edge(array,idx); }
456  void get_nodes(typename Node::array_type &array,
457  typename Face::index_type idx) const
458  { get_nodes_from_face(array,idx); }
459  void get_nodes(typename Node::array_type &array,
460  typename Cell::index_type idx) const
461  { get_nodes_from_cell(array,idx); }
462 
463  void get_edges(typename Edge::array_type &array,
464  typename Node::index_type idx) const
465  { get_edges_from_node(array,idx); }
466  void get_edges(typename Edge::array_type &array,
467  typename Edge::index_type idx) const
468  { array.resize(1); array[0]= idx; }
469  void get_edges(typename Edge::array_type &array,
470  typename Face::index_type idx) const
471  { get_edges_from_face(array,idx); }
472  void get_edges(typename Edge::array_type &array,
473  typename Cell::index_type idx) const
474  { get_edges_from_cell(array,idx); }
475 
476  void get_faces(typename Face::array_type &array,
477  typename Node::index_type idx) const
478  { get_faces_from_node(array,idx); }
479  void get_faces(typename Face::array_type &array,
480  typename Edge::index_type idx) const
481  { get_faces_from_edge(array,idx); }
482 
483  void get_faces(typename Face::array_type &array,
484  typename Face::index_type idx) const
485  { array.resize(1); array[0]= idx; }
486  void get_faces(typename Face::array_type &array,
487  typename Cell::index_type idx) const
488  { get_faces_from_cell(array,idx); }
489 
490  void get_cells(typename Cell::array_type &array,
491  typename Node::index_type idx) const
492  { get_cells_from_node(array,idx); }
493  void get_cells(typename Cell::array_type &array,
494  typename Edge::index_type idx) const
495  { get_cells_from_edge(array,idx); }
496  void get_cells(typename Cell::array_type &array,
497  typename Face::index_type idx) const
498  { get_cells_from_face(array,idx); }
499  void get_cells(typename Cell::array_type &array,
500  typename Cell::index_type idx) const
501  { array.resize(1); array[0]= idx; }
502 
503  void get_elems(typename Elem::array_type &array,
504  typename Node::index_type idx) const
505  { get_cells_from_node(array,idx); }
506  void get_elems(typename Elem::array_type &array,
507  typename Edge::index_type idx) const
508  { get_cells_from_edge(array,idx); }
509  void get_elems(typename Elem::array_type &array,
510  typename Face::index_type idx) const
511  { get_cells_from_face(array,idx); }
512  void get_elems(typename Elem::array_type &array,
513  typename Cell::index_type idx) const
514  { array.resize(1); array[0]= idx; }
515 
516  void get_delems(typename DElem::array_type &array,
517  typename Node::index_type idx) const
518  { get_faces_from_node(array,idx); }
519  void get_delems(typename DElem::array_type &array,
520  typename Edge::index_type idx) const
521  { get_faces_from_edge(array,idx); }
522  void get_delems(typename DElem::array_type &array,
523  typename Face::index_type idx) const
524  { array.resize(1); array[0]= idx; }
525  void get_delems(typename DElem::array_type &array,
526  typename Cell::index_type idx) const
527  { get_faces_from_cell(array,idx); }
528 
529  /// Generate the list of points that make up a sufficiently accurate
530  /// piecewise linear approximation of an edge.
531  template<class VECTOR, class INDEX>
532  void pwl_approx_edge(std::vector<VECTOR > &coords,
533  INDEX ci,
534  unsigned which_edge,
535  unsigned div_per_unit) const
536  {
537  basis_.approx_edge(which_edge, div_per_unit, coords);
538  }
539 
540  /// Generate the list of points that make up a sufficiently accurate
541  /// piecewise linear approximation of an face.
542  template<class VECTOR, class INDEX>
543  void pwl_approx_face(std::vector<std::vector<VECTOR > > &coords,
544  INDEX ci,
545  unsigned which_face,
546  unsigned div_per_unit) const
547  {
548  basis_.approx_face(which_face, div_per_unit, coords);
549  }
550 
551  /// get the center point (in object space) of an element
552  void get_center(Core::Geometry::Point &result, typename Node::index_type idx) const
553  { get_node_center(result, idx); }
554  void get_center(Core::Geometry::Point &result, typename Edge::index_type idx) const
555  { get_edge_center(result, idx); }
556  void get_center(Core::Geometry::Point &result, typename Face::index_type idx) const
557  { get_face_center(result, idx); }
558  void get_center(Core::Geometry::Point &result, typename Cell::index_type idx) const
559  { get_cell_center(result, idx); }
560 
561  /// Get the size of an element (length, area, volume)
562  double get_size(typename Node::index_type /*idx*/) const
563  { return 0.0; }
564 
565  double get_size(typename Edge::index_type idx) const
566  {
567  typename Node::array_type arr;
568  get_nodes(arr, idx);
569  Core::Geometry::Point p0, p1;
570  get_center(p0, arr[0]);
571  get_center(p1, arr[1]);
572  return ((p1 - p0).length());
573  }
574 
575  double get_size(typename Face::index_type idx) const
576  {
577  typename Node::array_type ra;
578  get_nodes(ra,idx);
579  Core::Geometry::Point p0,p1,p2;
580  get_point(p0,ra[0]);
581  get_point(p1,ra[1]);
582  get_point(p2,ra[2]);
583  return (Cross(p0-p1,p2-p0)).length()*0.5;
584  }
585 
586  double get_size(typename Elem::index_type idx) const
587  {
588  ElemData ed(*this,idx);
589  return (basis_.get_volume(ed));
590  }
591 
592  /// More specific names for get_size
593  double get_length(typename Edge::index_type idx) const
594  { return get_size(idx); };
595  double get_area(typename Face::index_type idx) const
596  { return get_size(idx); };
597  double get_volume(typename Cell::index_type idx) const
598  { return get_size(idx); };
599 
600  /// Get neighbors of an element or a node
601 
602  /// THIS ONE IS FLAWED AS IN 3D SPACE FOR AND ELEMENT TYPE THAT
603  /// IS NOT A VOLUME. HENCE IT WORKS HERE, BUT GENERALLY IT IS FLAWED
604  /// AS IT ASSUMES ONLY ONE NEIGHBOR, WHEREAS FOR ANYTHING ELSE THAN
605  /// A FACE THERE CAN BE MULTIPLE
606  bool get_neighbor(typename Elem::index_type &neighbor,
607  typename Elem::index_type elem,
608  typename DElem::index_type delem) const
609  { return(get_elem_neighbor(neighbor,elem,delem)); }
610 
611  /// These are more general implementations
612  void get_neighbors(std::vector<typename Node::index_type> &array,
613  typename Node::index_type node) const
614  { get_node_neighbors(array,node); }
615  bool get_neighbors(std::vector<typename Elem::index_type> &array,
616  typename Elem::index_type elem,
617  typename DElem::index_type delem) const
618  { return(get_elem_neighbors(array,elem,delem)); }
619  void get_neighbors(typename Elem::array_type &array,
620  typename Elem::index_type elem) const
621  { get_elem_neighbors(array,elem); }
622 
623  /// return false if point is out of range.
624  bool locate(typename Node::index_type &node, const Core::Geometry::Point &p) const
625  { return (locate_node(node,p)); }
626  bool locate(typename Edge::index_type &edge, const Core::Geometry::Point &p) const
627  { return (locate_edge(edge,p)); }
628  bool locate(typename Face::index_type &face, const Core::Geometry::Point &p) const
629  { return (locate_face(face,p)); }
630  bool locate(typename Cell::index_type &cell, const Core::Geometry::Point &p) const
631  { return (locate_elem(cell,p)); }
632 
633  bool locate(typename Elem::index_type &elem,
634  std::vector<double>& coords,
635  const Core::Geometry::Point &p) const
636  { return (locate_elem(elem,coords,p)); }
637 
638  /// These should become obsolete soon, they do not follow the concept
639  /// of the basis functions....
640  int get_weights(const Core::Geometry::Point &p, typename Node::array_type &l, double *w);
641  int get_weights(const Core::Geometry::Point&, typename Edge::array_type&, double*)
642  {ASSERTFAIL("TetVolMesh::get_weights for edges isn't supported"); }
643  int get_weights(const Core::Geometry::Point&, typename Face::array_type&, double*)
644  {ASSERTFAIL("TetVolMesh::get_weights for faces isn't supported"); }
645  int get_weights(const Core::Geometry::Point &p, typename Cell::array_type &l, double *w);
646 
647  /// Access the nodes of the mesh
648  void get_point(Core::Geometry::Point &result, typename Node::index_type index) const
649  { result = points_[index]; }
650  void set_point(const Core::Geometry::Point &point, typename Node::index_type index)
651  { points_[index] = point; }
652  void get_random_point(Core::Geometry::Point &p, typename Elem::index_type i, FieldRNG &r) const;
653 
654  /// Normals for visualizations
655  void get_normal(Core::Geometry::Vector& /*result*/, typename Node::index_type /*index*/) const
656  { ASSERTFAIL("TetVolMesh: This mesh type does not have node normals."); }
657 
658  /// Get the normals at the outside of the element
659  template<class VECTOR, class INDEX1, class INDEX2>
660  void get_normal(Core::Geometry::Vector &result, VECTOR& coords,
661  INDEX1 eidx, INDEX2 fidx)
662  {
663  // Improved algorithm, which should be faster as it is fully
664  // on the stack.
665 
666  // Obtain the inverse jacobian
667  double Ji[9];
668  inverse_jacobian(coords,eidx,Ji);
669 
670  // Get the normal in local coordinates
671  const double un0 = basis_.unit_face_normals[fidx][0];
672  const double un1 = basis_.unit_face_normals[fidx][1];
673  const double un2 = basis_.unit_face_normals[fidx][2];
674  // Do the matrix multiplication: should result in a vector
675  // in the global coordinate space
676  result.x(Ji[0]*un0+Ji[1]*un1+Ji[2]*un2);
677  result.y(Ji[3]*un0+Ji[4]*un1+Ji[5]*un2);
678  result.z(Ji[6]*un0+Ji[7]*un1+Ji[8]*un2);
679 
680  // normalize vector
681  result.normalize();
682  }
683 
684  /// Add a new node to the mesh
685  typename Node::index_type add_point(const Core::Geometry::Point &p);
686  typename Node::index_type add_node(const Core::Geometry::Point &p)
687  { return(add_point(p)); }
688 
689  /// Add a new element to the mesh
690  template<class ARRAY>
691  typename Elem::index_type add_elem(ARRAY a)
692  {
693  ASSERTMSG(a.size() == 4, "Tried to add non-tet element.");
694 
695  return add_tet( static_cast<typename Node::index_type>(a[0]),
696  static_cast<typename Node::index_type>(a[1]),
697  static_cast<typename Node::index_type>(a[2]),
698  static_cast<typename Node::index_type>(a[3]) );
699  }
700 
701  /// Functions to improve memory management. Often one knows how many
702  /// nodes/elements one needs, prereserving memory is often possible.
703  void node_reserve(size_type s) { points_.reserve(static_cast<std::vector<Core::Geometry::Point>::size_type>(s)); }
704  void elem_reserve(size_type s) { cells_.reserve(static_cast<std::vector<index_type>::size_type>(s*4)); }
705  void resize_nodes(size_type s) { points_.resize(static_cast<std::vector<Core::Geometry::Point>::size_type>(s)); }
706  void resize_elems(size_type s) { cells_.resize(static_cast<std::vector<index_type>::size_type>(s*4)); }
707 
708  /// Get the local coordinates for a certain point within an element
709  /// This function uses a couple of newton iterations to find the local
710  /// coordinate of a point
711  template<class VECTOR, class INDEX>
712  bool get_coords(VECTOR &coords, const Core::Geometry::Point &p, INDEX idx) const
713  {
714  ElemData ed(*this, idx);
715  return basis_.get_coords(coords, p, ed);
716  }
717 
718  /// Find the location in the global coordinate system for a local coordinate
719  /// This function is the opposite of get_coords.
720  template<class VECTOR, class INDEX>
721  void interpolate(Core::Geometry::Point &pt, const VECTOR &coords, INDEX idx) const
722  {
723  ElemData ed(*this, idx);
724  pt = basis_.interpolate(coords, ed);
725  }
726 
727  /// Interpolate the derivate of the function, This infact will return the
728  /// jacobian of the local to global coordinate transformation. This function
729  /// is mainly intended for the non linear elements
730  template<class VECTOR1, class INDEX, class VECTOR2>
731  void derivate(const VECTOR1 &coords, INDEX idx, VECTOR2 &J) const
732  {
733  ElemData ed(*this, idx);
734  basis_.derivate(coords, ed, J);
735  }
736 
737  /// Get the determinant of the jacobian, which is the local volume of an element
738  /// and is intended to help with the integration of functions over an element.
739  template<class VECTOR, class INDEX>
740  double det_jacobian(const VECTOR& coords, INDEX idx) const
741  {
742  StackVector<Core::Geometry::Point,3> Jv;
743  ElemData ed(*this,idx);
744  basis_.derivate(coords,ed,Jv);
745  return (DetMatrix3P(Jv));
746  }
747 
748  /// Get the jacobian of the transformation. In case one wants the non inverted
749  /// version of this matrix. This is currentl here for completeness of the
750  /// interface
751  template<class VECTOR, class INDEX>
752  void jacobian(const VECTOR& coords,INDEX idx, double* J) const
753  {
754  StackVector<Core::Geometry::Point,3> Jv;
755  ElemData ed(*this,idx);
756  basis_.derivate(coords,ed,Jv);
757  J[0] = Jv[0].x();
758  J[1] = Jv[0].y();
759  J[2] = Jv[0].z();
760  J[3] = Jv[1].x();
761  J[4] = Jv[1].y();
762  J[5] = Jv[1].z();
763  J[6] = Jv[2].x();
764  J[7] = Jv[2].y();
765  J[8] = Jv[2].z();
766  }
767 
768  /// Get the inverse jacobian of the transformation. This one is needed to
769  /// translate local gradients into global gradients. Hence it is crucial for
770  /// calculating gradients of fields, or constructing finite elements.
771  template<class VECTOR, class INDEX>
772  double inverse_jacobian(const VECTOR& coords, INDEX idx, double* Ji) const
773  {
774  StackVector<Core::Geometry::Point,3> Jv;
775  ElemData ed(*this,idx);
776  basis_.derivate(coords,ed,Jv);
777  return (InverseMatrix3P(Jv,Ji));
778  }
779 
780  template<class INDEX>
781  double scaled_jacobian_metric(INDEX idx) const
782  {
783  StackVector<Core::Geometry::Point,3> Jv;
784  ElemData ed(*this,idx);
785 
786  double temp;
787 
788  basis_.derivate(basis_.unit_center,ed,Jv);
789  double min_jacobian = DetMatrix3P(Jv);
790  INDEX idx2 = idx*4;
791  Core::Geometry::Point p0 = points_[cells_[idx2]];
792  Core::Geometry::Point p1 = points_[cells_[idx2+1]];
793  Core::Geometry::Point p2 = points_[cells_[idx2+2]];
794  Core::Geometry::Point p3 = points_[cells_[idx2+3]];
795 
796  double l0 = (p1-p0).length();
797  double l1 = (p2-p1).length();
798  double l2 = (p0-p2).length();
799  double l3 = (p3-p0).length();
800  double l4 = (p3-p1).length();
801  double l5 = (p3-p2).length();
802 
803  double min_scale = l0*l2*l3;
804  temp = l0*l1*l4;
805  if (temp > min_scale) min_scale = temp;
806  temp = l1*l2*l5;
807  if (temp > min_scale) min_scale = temp;
808  temp = l3*l4*l5;
809  if (temp > min_scale) min_scale = temp;
810  if (min_jacobian > min_scale) min_scale = min_jacobian;
811 
812  return (sqrt(2.0)*min_jacobian/(min_scale));
813  }
814 
815  template<class INDEX>
816  double jacobian_metric(INDEX idx) const
817  {
818  StackVector<Core::Geometry::Point,3> Jv;
819  ElemData ed(*this,idx);
820 
821  double temp;
822 
823  basis_.derivate(basis_.unit_center,ed,Jv);
824  double min_jacobian = DetMatrix3P(Jv);
825 
826  size_t num_vertices = basis_.number_of_vertices();
827  for (size_t j=0;j < num_vertices;j++)
828  {
829  basis_.derivate(basis_.unit_vertices[j],ed,Jv);
830  temp = DetMatrix3P(Jv);
831  if(temp < min_jacobian) min_jacobian = temp;
832  }
833 
834  return (min_jacobian);
835  }
836 
837 
838  template<class INDEX>
839  double inscribed_circumscribed_radius_metric(INDEX idx) const
840  {
841  StackVector<Core::Geometry::Point,3> Vec;
842  ElemData ed(*this,idx);
843 
844  INDEX idx2 = idx*4;
845  Core::Geometry::Point p0 = points_[cells_[idx2]];
846  Core::Geometry::Point p1 = points_[cells_[idx2+1]];
847  Core::Geometry::Point p2 = points_[cells_[idx2+2]];
848  Core::Geometry::Point p3 = points_[cells_[idx2+3]];
849 
850  double l0 = (p1-p0).length();
851  double l1 = (p2-p1).length();
852  double l2 = (p0-p2).length();
853  double l3 = (p3-p0).length();
854  double l4 = (p3-p1).length();
855  double l5 = (p3-p2).length();
856 
857  // area of a triangle using heron's formula
858  // heron's area = sqrt(s(s-a)(s-b)(s-c)) where s=0.5*(a+b+c)
859  double tri_AREA = (l0*l5 + l1*l3 + l2*l4) / 2.0;
860 
861  // volume of the element
862  double tet_VOL = basis_.get_volume(ed);
863 
864  // surface area of element
865  double s_A = 0.5*(l0+l1+l2);
866  double s_B = 0.5*(l1+l4+l5);
867  double s_C = 0.5*(l2+l3+l5);
868  double s_D = 0.5*(l0+l3+l4);
869 
870  double tri_A = sqrt(s_A*(s_A - l0)*(s_A - l1)*(s_A - l2));
871  double tri_B = sqrt(s_B*(s_B - l1)*(s_B - l4)*(s_B - l5));
872  double tri_C = sqrt(s_C*(s_C - l2)*(s_C - l3)*(s_C - l5));
873  double tri_D = sqrt(s_D*(s_D - l0)*(s_D - l3)*(s_D - l4));
874 
875  double tet_SA = tri_A + tri_B + tri_C + tri_D;
876 
877  // calculating inscribed and circumscribed radii
878  double r_in = 3.0 * tet_VOL / tet_SA;
879  double r_cir = sqrt(tri_AREA*(tri_AREA - l0*l5)*(tri_AREA - l1*l3)*(tri_AREA - l2*l4)) / (6.0*tet_VOL);
880 
881  // for an equilateral tetrahedron, the circ radius is 3 times larger than the inscr
882  double ideal_ratio = 0.333333;
883  double normalized_actual_ratio = (r_in / r_cir) / ideal_ratio;
884  return (normalized_actual_ratio);
885  }
886 
887  template <class INDEX>
888  bool find_closest_node(double& pdist, Core::Geometry::Point &result,
889  INDEX &node, const Core::Geometry::Point &p) const
890  {
891  return(find_closest_node(pdist,result,node,p,-1.0));
892  }
893 
894  template <class INDEX>
895  bool find_closest_node(double& pdist, Core::Geometry::Point &result,
896  INDEX &node, const Core::Geometry::Point &p,
897  double maxdist) const
898  {
899  if (maxdist < 0.0) maxdist = DBL_MAX; else maxdist = maxdist*maxdist;
900  typename Node::size_type sz; size(sz);
901 
902  /// If there are no nodes we cannot find the closest one
903  if (sz == 0) return (false);
904 
905  if (node >= 0 && node < sz)
906  {
907  Core::Geometry::Point point = points_[node];
908  double dist = (point-p).length2();
909 
910  if ( dist < epsilon2_ )
911  {
912  result = point;
913  pdist = sqrt(dist);
914  return (true);
915  }
916  }
917 
918  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
919  "TetVolMesh::find_closest_node requires synchronize(NODE_LOCATE_E).");
920 
921  // get grid sizes
922  const size_type ni = node_grid_->get_ni()-1;
923  const size_type nj = node_grid_->get_nj()-1;
924  const size_type nk = node_grid_->get_nk()-1;
925 
926  // Convert to grid coordinates.
927  index_type bi, bj, bk, ei, ej, ek;
928  node_grid_->unsafe_locate(bi, bj, bk, p);
929 
930  // Clamp to closest point on the grid.
931  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
932  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
933  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
934 
935  ei = bi; ej = bj; ek = bk;
936 
937  double dmin = maxdist;
938  bool found = true;
939  bool found_one = false;
940  do
941  {
942  found = true;
943  /// We need to do a full shell without any elements that are closer
944  /// to make sure there no closer elements in neighboring searchgrid cells
945 
946  for (index_type i = bi; i <= ei; i++)
947  {
948  if (i < 0 || i > ni) continue;
949  for (index_type j = bj; j <= ej; j++)
950  {
951  if (j < 0 || j > nj) continue;
952  for (index_type k = bk; k <= ek; k++)
953  {
954  if (k < 0 || k > nk) continue;
955  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
956  {
957  if (node_grid_->min_distance_squared(p, i, j, k) < dmin)
958  {
959  typename SearchGridT<index_type>::iterator it, eit;
960  found = false;
961  node_grid_->lookup_ijk(it,eit, i, j, k);
962 
963  while (it != eit)
964  {
965  const Core::Geometry::Point point = points_[*it];
966  const double dist = (p-point).length2();
967 
968  if (dist < dmin)
969  {
970  found_one = true;
971  result = point;
972  node = INDEX(*it);
973  dmin = dist;
974  /// If we are closer than eps^2 we found a node close enough
975  if (dmin < epsilon2_)
976  {
977  pdist = sqrt(dmin);
978  return (true);
979  }
980  }
981  ++it;
982  }
983  }
984  }
985  }
986  }
987  }
988  bi--;ei++;
989  bj--;ej++;
990  bk--;ek++;
991  }
992  while (!found) ;
993 
994  if (!found_one) return (false);
995  pdist = sqrt(dmin);
996  return (true);
997  }
998 
999  template <class ARRAY>
1000  bool find_closest_nodes(ARRAY &nodes, const Core::Geometry::Point &p, double maxdist) const
1001  {
1002  nodes.clear();
1003 
1004  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
1005  "TetVolMesh::find_closest_node requires synchronize(NODE_LOCATE_E).")
1006 
1007  // get grid sizes
1008  const size_type ni = node_grid_->get_ni()-1;
1009  const size_type nj = node_grid_->get_nj()-1;
1010  const size_type nk = node_grid_->get_nk()-1;
1011 
1012  // Convert to grid coordinates.
1013  index_type bi, bj, bk, ei, ej, ek;
1014 
1015  Core::Geometry::Point max = p+Core::Geometry::Vector(maxdist,maxdist,maxdist);
1016  Core::Geometry::Point min = p+Core::Geometry::Vector(-maxdist,-maxdist,-maxdist);
1017 
1018  node_grid_->unsafe_locate(bi, bj, bk, min);
1019  node_grid_->unsafe_locate(ei, ej, ek, max);
1020 
1021  // Clamp to closest point on the grid.
1022  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
1023  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
1024  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
1025 
1026  if (ei > ni) ei = ni; if (ei < 0) ei = 0;
1027  if (ej > nj) ej = nj; if (ej < 0) ej = 0;
1028  if (ek > nk) ek = nk; if (ek < 0) ek = 0;
1029 
1030  double maxdist2 = maxdist*maxdist;
1031 
1032  for (index_type i = bi; i <= ei; i++)
1033  {
1034  for (index_type j = bj; j <= ej; j++)
1035  {
1036  for (index_type k = bk; k <= ek; k++)
1037  {
1038  if (node_grid_->min_distance_squared(p, i, j, k) < maxdist2)
1039  {
1040  typename SearchGridT<index_type>::iterator it, eit;
1041  node_grid_->lookup_ijk(it, eit, i, j, k);
1042 
1043  while (it != eit)
1044  {
1045  const Core::Geometry::Point point = points_[*it];
1046  const double dist = (p-point).length2();
1047 
1048  if (dist < maxdist2)
1049  {
1050  nodes.push_back(*it);
1051  }
1052  ++it;
1053  }
1054  }
1055  }
1056  }
1057  }
1058 
1059  return(nodes.size() > 0);
1060  }
1061 
1062 
1063 
1064  template <class ARRAY1, class ARRAY2>
1065  bool find_closest_nodes(ARRAY1 &distances, ARRAY2 &nodes, const Core::Geometry::Point &p, double maxdist) const
1066  {
1067  nodes.clear();
1068  distances.clear();
1069 
1070  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
1071  "TetVolMesh::find_closest_node requires synchronize(NODE_LOCATE_E).")
1072 
1073  // get grid sizes
1074  const size_type ni = node_grid_->get_ni()-1;
1075  const size_type nj = node_grid_->get_nj()-1;
1076  const size_type nk = node_grid_->get_nk()-1;
1077 
1078  // Convert to grid coordinates.
1079  index_type bi, bj, bk, ei, ej, ek;
1080 
1081  Core::Geometry::Point max = p+Core::Geometry::Vector(maxdist,maxdist,maxdist);
1082  Core::Geometry::Point min = p+Core::Geometry::Vector(-maxdist,-maxdist,-maxdist);
1083 
1084  node_grid_->unsafe_locate(bi, bj, bk, min);
1085  node_grid_->unsafe_locate(ei, ej, ek, max);
1086 
1087  // Clamp to closest point on the grid.
1088  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
1089  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
1090  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
1091 
1092  if (ei > ni) ei = ni; if (ei < 0) ei = 0;
1093  if (ej > nj) ej = nj; if (ej < 0) ej = 0;
1094  if (ek > nk) ek = nk; if (ek < 0) ek = 0;
1095 
1096  double maxdist2 = maxdist*maxdist;
1097 
1098  for (index_type i = bi; i <= ei; i++)
1099  {
1100  for (index_type j = bj; j <= ej; j++)
1101  {
1102  for (index_type k = bk; k <= ek; k++)
1103  {
1104  if (node_grid_->min_distance_squared(p, i, j, k) < maxdist2)
1105  {
1106  typename SearchGridT<index_type>::iterator it, eit;
1107  node_grid_->lookup_ijk(it, eit, i, j, k);
1108 
1109  while (it != eit)
1110  {
1111  const Core::Geometry::Point point = points_[*it];
1112  const double dist = (p-point).length2();
1113 
1114  if (dist < maxdist2)
1115  {
1116  nodes.push_back(*it);
1117  distances.push_back(dist);
1118  }
1119  ++it;
1120  }
1121  }
1122  }
1123  }
1124  }
1125 
1126  return(nodes.size() > 0);
1127  }
1128 
1129 
1130  /// Find the closest element to a point
1131  template <class INDEX, class ARRAY>
1132  bool find_closest_elem(double& pdist,
1133  Core::Geometry::Point &result,
1134  ARRAY &coords,
1135  INDEX &elem,
1136  const Core::Geometry::Point &p) const
1137  {
1138  return (find_closest_elem(pdist,result,coords,elem,p,-1.0));
1139  }
1140 
1141 
1142  /// Find the closest element to a point
1143  template <class INDEX, class ARRAY>
1144  bool find_closest_elem(double& pdist,
1145  Core::Geometry::Point &result,
1146  ARRAY &coords,
1147  INDEX &elem,
1148  const Core::Geometry::Point &p,
1149  double maxdist) const
1150  {
1151  if (maxdist < 0.0) maxdist = DBL_MAX; else maxdist = maxdist*maxdist;
1152 
1153  typename Elem::size_type sz; size(sz);
1154 
1155  /// If there are no nodes we cannot find a closest point
1156  if (sz == 0) return (false);
1157 
1158  /// Check whether the estimate given in idx is the point we are looking for
1159  if ((elem > 0)&&(elem < sz))
1160  {
1161  if (inside(elem,p))
1162  {
1163  ElemData ed(*this, elem);
1164  basis_.get_coords(coords, p, ed);
1165  result = p;
1166  pdist = 0.0;
1167  return (true);
1168  }
1169  }
1170 
1171  ASSERTMSG(synchronized_ & Mesh::FACES_E,
1172  "TetVolMesh: need to synchronize FACES_E first");
1173  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
1174  "TetVolMesh: need to synchronize ELEM_LOCATE_E first");
1175 
1176  // First check are we inside an element
1177  SearchGridT<index_type>::iterator it, eit;
1178  if (elem_grid_->lookup(it, eit, p))
1179  {
1180  while (it != eit)
1181  {
1182  if (inside(typename Elem::index_type(*it), p))
1183  {
1184  pdist = 0.0;
1185  result = p;
1186  elem = static_cast<INDEX>(*it);
1187  ElemData ed(*this, elem);
1188  basis_.get_coords(coords, p, ed);
1189  return (true);
1190  }
1191  ++it;
1192  }
1193  }
1194 
1195  // If not start searching for the closest outer boundary
1196  // get grid sizes
1197  const size_type ni = elem_grid_->get_ni()-1;
1198  const size_type nj = elem_grid_->get_nj()-1;
1199  const size_type nk = elem_grid_->get_nk()-1;
1200 
1201  // Convert to grid coordinates.
1202  index_type bi, ei, bj, ej, bk, ek;
1203  elem_grid_->unsafe_locate(bi, bj, bk, p);
1204 
1205  // Clamp to closest point on the grid.
1206  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
1207  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
1208  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
1209 
1210  ei = bi; ej = bj; ek = bk;
1211 
1212  double dmin = maxdist;
1213  bool found = true;
1214  bool found_one = false;
1215 
1216  do
1217  {
1218  found = true;
1219  /// We need to do a full shell without any elements that are closer
1220  /// to make sure there no closer elements in neighboring searchgrid cells
1221  for (index_type i = bi; i <= ei; i++)
1222  {
1223  if (i < 0 || i > ni) continue;
1224  for (index_type j = bj; j <= ej; j++)
1225  {
1226  if (j < 0 || j > nj) continue;
1227  for (index_type k = bk; k <= ek; k++)
1228  {
1229  if (k < 0 || k > nk) continue;
1230  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
1231  {
1232  if (elem_grid_->min_distance_squared(p, i, j, k) < dmin)
1233  {
1234  found = false;
1235  typename SearchGridT<index_type>::iterator it, eit;
1236  elem_grid_->lookup_ijk(it,eit, i, j, k);
1237 
1238  while (it != eit)
1239  {
1240  Core::Geometry::Point r;
1241  index_type cidx = (*it);
1242  index_type idx = cidx*4;
1243  unsigned char b = boundary_faces_[cidx];
1244  if (b)
1245  {
1246  if (b & 0x1)
1247  {
1248  closest_point_on_tri(r, p,
1249  points_[cells_[idx ]],
1250  points_[cells_[idx+2]],
1251  points_[cells_[idx+1]]);
1252  const double dtmp = (p - r).length2();
1253  if (dtmp < dmin)
1254  {
1255  found_one = true;
1256  result = r;
1257  elem = INDEX(cidx);
1258  dmin = dtmp;
1259 
1260  if (dmin < epsilon2_)
1261  {
1262  pdist = sqrt(dmin);
1263  ElemData ed(*this,elem);
1264  basis_.get_coords(coords,result,ed);
1265  return (true);
1266  }
1267  }
1268  }
1269  if (b & 0x2)
1270  {
1271  closest_point_on_tri(r, p,
1272  points_[cells_[idx+1]],
1273  points_[cells_[idx+2]],
1274  points_[cells_[idx+3]]);
1275  const double dtmp = (p - r).length2();
1276  if (dtmp < dmin)
1277  {
1278  found_one = true;
1279  result = r;
1280  elem = INDEX(cidx);
1281  dmin = dtmp;
1282 
1283  if (dmin < epsilon2_)
1284  {
1285  pdist = sqrt(dmin);
1286  ElemData ed(*this,elem);
1287  basis_.get_coords(coords,result,ed);
1288  return (true);
1289  }
1290  }
1291  }
1292  if (b & 0x4)
1293  {
1294  closest_point_on_tri(r, p,
1295  points_[cells_[idx]],
1296  points_[cells_[idx+1]],
1297  points_[cells_[idx+3]]);
1298  const double dtmp = (p - r).length2();
1299  if (dtmp < dmin)
1300  {
1301  found_one = true;
1302  result = r;
1303  elem = INDEX(cidx);
1304  dmin = dtmp;
1305 
1306  if (dmin < epsilon2_)
1307  {
1308  pdist = sqrt(dmin);
1309  ElemData ed(*this,elem);
1310  basis_.get_coords(coords,result,ed);
1311  return (true);
1312  }
1313  }
1314  }
1315  if (b & 0x8)
1316  {
1317  closest_point_on_tri(r, p,
1318  points_[cells_[idx]],
1319  points_[cells_[idx+3]],
1320  points_[cells_[idx+2]]);
1321  const double dtmp = (p - r).length2();
1322  if (dtmp < dmin)
1323  {
1324  found_one = true;
1325  result = r;
1326  elem = INDEX(cidx);
1327  dmin = dtmp;
1328 
1329  if (dmin < epsilon2_)
1330  {
1331  pdist = sqrt(dmin);
1332  ElemData ed(*this,elem);
1333  basis_.get_coords(coords,result,ed);
1334  return (true);
1335  }
1336  }
1337  }
1338  }
1339  ++it;
1340  }
1341  }
1342  }
1343  }
1344  }
1345  }
1346  bi--;ei++;
1347  bj--;ej++;
1348  bk--;ek++;
1349  }
1350  while (!found) ;
1351 
1352  if (!found_one) return (false);
1353 
1354  ElemData ed(*this,elem);
1355  basis_.get_coords(coords,result,ed);
1356 
1357  pdist = sqrt(dmin);
1358  return (true);
1359  }
1360 
1361  template <class INDEX>
1362  bool find_closest_elem(double& pdist, Core::Geometry::Point &result,
1363  INDEX &elem, const Core::Geometry::Point &p) const
1364  {
1365  StackVector<double,3> coords;
1366  return(find_closest_elem(pdist,result,coords,elem,p,-1.0));
1367  }
1368 
1369  /// Find the closest elements to a point
1370  template<class ARRAY>
1371  bool find_closest_elems(double& /*pdist*/, Core::Geometry::Point& /*result*/,
1372  ARRAY& /*elems*/, const Core::Geometry::Point& /*p*/) const
1373  {
1374  ASSERTFAIL("TetVolMesh: Find closest element has not yet been implemented.");
1375  return (false);
1376  }
1377 
1378  /// Export this class using the old Pio system
1379  virtual void io(Piostream&);
1380 
1381  ///////////////////////////////////////////////////
1382  // STATIC VARIABLES AND FUNCTIONS
1383 
1384  /// This ID is created as soon as this class will be instantiated
1385  static PersistentTypeID tetvolmesh_typeid;
1386 
1387  /// Core functionality for getting the name of a templated mesh class
1388  static const std::string type_name(int n = -1);
1389  virtual std::string dynamic_type_name() const { return tetvolmesh_typeid.type; }
1390 
1391  /// Type description, used for finding names of the mesh class for
1392  /// dynamic compilation purposes. Soem of this should be obsolete
1393  virtual const TypeDescription *get_type_description() const;
1394  static const TypeDescription* node_type_description();
1395  static const TypeDescription* edge_type_description();
1396  static const TypeDescription* face_type_description();
1397  static const TypeDescription* cell_type_description();
1398  static const TypeDescription* elem_type_description()
1399  { return cell_type_description(); }
1400 
1401  /// This function returns a maker for Pio.
1402  static Persistent* maker() { return new TetVolMesh(); }
1403  /// This function returns a handle for the virtual interface.
1404  static MeshHandle mesh_maker() { return boost::make_shared<TetVolMesh>(); }
1405 
1406  //////////////////////////////////////////////////////////////////
1407  // Mesh specific functions (these are not implemented in every mesh)
1408 
1409  // Extra functionality needed by this specific geometry.
1410  void set_nodes(typename Node::array_type &,
1411  typename Cell::index_type);
1412 
1413  /// Deleting cells, make sure mesh is detached
1414  void delete_cells(std::set<index_type> &to_delete);
1415  void delete_nodes(std::set<index_type> &to_delete);
1416 
1417  /// Orient a cell properly
1418  void orient(typename Cell::index_type ci);
1419 
1420  /// MORE UNSORTED FUNCTIONS
1421  /// Used to recompute data for individual cells. Don't use these, they
1422  // are not synchronous. Use create_cell_syncinfo instead.
1423  void create_cell_edges(typename Cell::index_type);
1424  void delete_cell_edges(typename Cell::index_type, bool table_only = false);
1425  void create_cell_faces(typename Cell::index_type);
1426  void delete_cell_faces(typename Cell::index_type, bool table_only = false);
1427  void create_cell_node_neighbors(typename Cell::index_type);
1428  void delete_cell_node_neighbors(typename Cell::index_type);
1429 
1430  void create_cell_syncinfo(typename Cell::index_type ci);
1431  void delete_cell_syncinfo(typename Cell::index_type ci);
1432  // Create the partial sync info needed by insert_node_in_elem
1433  void create_cell_syncinfo_special(typename Cell::index_type ci);
1434  void delete_cell_syncinfo_special(typename Cell::index_type ci);
1435  // May reorient the tet to make the signed volume positive.
1436  typename Elem::index_type add_tet_pos(typename Node::index_type a,
1437  typename Node::index_type b,
1438  typename Node::index_type c,
1439  typename Node::index_type d);
1440  void mod_tet_pos(typename Elem::index_type tet,
1441  typename Node::index_type a,
1442  typename Node::index_type b,
1443  typename Node::index_type c,
1444  typename Node::index_type d);
1445 
1446  /// Functions needed for cubit interface
1447  /// WE SHOULD MAKE THESE GENERAL AND IN EVERY MESHTYPE
1448  template <class Iter, class Functor>
1449  void fill_points(Iter begin, Iter end, Functor fill_ftor);
1450  template <class Iter, class Functor>
1451  void fill_cells(Iter begin, Iter end, Functor fill_ftor);
1452  template <class Iter, class Functor>
1453  void fill_neighbors(Iter begin, Iter end, Functor fill_ftor);
1454  template <class Iter, class Functor>
1455  void fill_data(Iter begin, Iter end, Functor fill_ftor);
1456 
1457  // THIS FUNCTION NEEDS TO BE REVISED AS IT DOES NOT SCALE PROPERLY
1458  // THE TOLERANCE IS NOT RELATIVE, WHICH IS A PROBLEM.........
1459  typename Node::index_type add_find_point(const Core::Geometry::Point &p,
1460  double err = 1.0e-6);
1461 
1462  // Short cut for generating an element....
1463  typename Elem::index_type add_tet(typename Node::index_type a,
1464  typename Node::index_type b,
1465  typename Node::index_type c,
1466  typename Node::index_type d);
1467 
1468  typename Elem::index_type add_tet(const Core::Geometry::Point &p0,
1469  const Core::Geometry::Point &p1,
1470  const Core::Geometry::Point &p2,
1471  const Core::Geometry::Point &p3);
1472 
1473  bool insert_node_in_elem(typename Elem::array_type &tets,
1474  typename Node::index_type &ni,
1475  typename Elem::index_type ci,
1476  const Core::Geometry::Point &p);
1477 
1478  /// must detach, if altering points!
1479  std::vector<Core::Geometry::Point>& get_points() { return points_; }
1480 
1481  int compute_checksum();
1482 
1483 protected:
1484  //////////////////////////////////////////////////////////////
1485  // These functions are templates and are used to define the
1486  // dynamic compilation interface and the virtual interface
1487  // as they both use different datatypes as indices and arrays
1488  // the following functions have been templated and are inlined
1489  // at the places where they are needed.
1490  //
1491  // Secondly these templates allow for the use of the stack vector
1492  // as well as the STL vector. When an algorithm supports non linear
1493  // functions an STL vector is a better choice, in the other cases
1494  // often a StackVector is enough (The latter improves performance).
1495 
1496  template<class ARRAY, class INDEX>
1497  inline void get_nodes_from_edge(ARRAY& array, INDEX idx) const
1498  {
1499  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1500  "TetVolMesh: Must call synchronize EDGES_E first");
1501 
1502  if (edges_[idx].cells_.size() == 0)
1503  { array.clear(); return; }
1504 
1505  array.resize(2);
1506 
1507  index_type cell_edge_index = edges_[idx].cells_[0];
1508  index_type cell_index = (cell_edge_index>>3) << 2;
1509  index_type edge_index = (cell_edge_index)&0x7;
1510 
1511  const int *offset = TetVolEdgeTable[edge_index];
1512  array[0] = static_cast<typename ARRAY::value_type>(cells_[cell_index+offset[0]]);
1513  array[1] = static_cast<typename ARRAY::value_type>(cells_[cell_index+offset[1]]);
1514  }
1515 
1516  // Always returns nodes in counter-clockwise order
1517  template<class ARRAY, class INDEX>
1518  inline void get_nodes_from_face(ARRAY& array, INDEX idx) const
1519  {
1520  ASSERTMSG(synchronized_ & Mesh::FACES_E,
1521  "TetVolMesh: Must call synchronize FACES_E first");
1522 
1523  if (faces_[idx].cells_[0] == MESH_NO_NEIGHBOR)
1524  { array.clear(); return; }
1525 
1526  index_type cell_face_index = faces_[idx].cells_[0];
1527  index_type cell_index = cell_face_index&(~0x3);
1528  index_type face_index = cell_face_index&0x3;
1529 
1530  array.resize(3);
1531 
1532  const int *offset = TetVolFaceTable[face_index];
1533  array[0] = static_cast<typename ARRAY::value_type>(cells_[cell_index+offset[0]]);
1534  array[1] = static_cast<typename ARRAY::value_type>(cells_[cell_index+offset[1]]);
1535  array[2] = static_cast<typename ARRAY::value_type>(cells_[cell_index+offset[2]]);
1536  }
1537 
1538  template<class ARRAY, class INDEX>
1539  inline void get_nodes_from_cell(ARRAY& array, INDEX idx) const
1540  {
1541  array.resize(4);
1542  const index_type off = idx * 4;
1543  array[0] = static_cast<typename ARRAY::value_type>(cells_[off]);
1544  array[1] = static_cast<typename ARRAY::value_type>(cells_[off + 1]);
1545  array[2] = static_cast<typename ARRAY::value_type>(cells_[off + 2]);
1546  array[3] = static_cast<typename ARRAY::value_type>(cells_[off + 3]);
1547  }
1548 
1549  template<class ARRAY, class INDEX>
1550  inline void get_nodes_from_elem(ARRAY& array, INDEX idx) const
1551  {
1552  get_nodes_from_cell(array,idx);
1553  }
1554 
1555  template<class ARRAY, class INDEX>
1556  inline void get_edges_from_face(ARRAY& array, INDEX idx) const
1557  {
1558  ASSERTMSG(synchronized_ & (Mesh::EDGES_E),
1559  "TetVolMesh: Must call synchronize EDGES_E first");
1560  ASSERTMSG(synchronized_ & (Mesh::FACES_E),
1561  "TetVolMesh: Must call synchronize FACES_E first");
1562 
1563  array.clear();
1564  array.reserve(3);
1565  index_type cell_index = (faces_[idx].cells_[0])&(~0x3);
1566  index_type face_index = (faces_[idx].cells_[0])&(0x3);
1567 
1568  const int* offset = TetVolFaceTable[face_index];
1569 
1570  typename Node::index_type n0 = cells_[cell_index+offset[0]];
1571  typename Node::index_type n1 = cells_[cell_index+offset[1]];
1572  typename Node::index_type n2 = cells_[cell_index+offset[2]];
1573 
1574  if (n0 != n1)
1575  {
1576  PEdgeNode e(n0, n1);
1577  array.push_back(static_cast<typename ARRAY::value_type>(
1578  (*(edge_table_.find(e))).second));
1579  }
1580  if (n1 != n2)
1581  {
1582  PEdgeNode e(n1, n2);
1583  array.push_back(static_cast<typename ARRAY::value_type>(
1584  (*(edge_table_.find(e))).second));
1585  }
1586  if (n2 != n0)
1587  {
1588  PEdgeNode e(n2, n0);
1589  array.push_back(static_cast<typename ARRAY::value_type>(
1590  (*(edge_table_.find(e))).second));
1591  }
1592  }
1593 
1594  template<class ARRAY, class INDEX>
1595  inline void get_edges_from_cell(ARRAY& array, INDEX idx) const
1596  {
1597  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1598  "TetVolMesh: Must call synchronize EDGES_E first");
1599 
1600  const index_type off = idx * 4;
1601  PEdgeNode e00(cells_[off + 0], cells_[off + 1]);
1602  PEdgeNode e01(cells_[off + 1], cells_[off + 2]);
1603  PEdgeNode e02(cells_[off + 2], cells_[off + 0]);
1604  PEdgeNode e03(cells_[off + 0], cells_[off + 3]);
1605  PEdgeNode e04(cells_[off + 1], cells_[off + 3]);
1606  PEdgeNode e05(cells_[off + 2], cells_[off + 3]);
1607  typename Node::index_type n1,n2;
1608 
1609  n1 = cells_[off ]; n2 = cells_[off + 1];
1610  size_t i = 0;
1611  typedef typename ARRAY::value_type T;
1612  if (n1 != n2)
1613  {
1614  PEdge e(n1,n2);
1615  array[i++] = (static_cast<T>((*(edge_table_.find(e))).second));
1616  }
1617  n1 = cells_[off + 1]; n2 = cells_[off + 2];
1618  if (n1 != n2)
1619  {
1620  PEdge e(n1,n2);
1621  array[i++] = (static_cast<T>((*(edge_table_.find(e))).second));
1622  }
1623  n1 = cells_[off + 2]; n2 = cells_[off ];
1624  if (n1 != n2)
1625  {
1626  PEdge e(n1,n2);
1627  array[i++] = (static_cast<T>((*(edge_table_.find(e))).second));
1628  }
1629  n1 = cells_[off ]; n2 = cells_[off + 3];
1630  if (n1 != n2)
1631  {
1632  PEdge e(n1,n2);
1633  array[i++] = (static_cast<T>((*(edge_table_.find(e))).second));
1634  }
1635  n1 = cells_[off + 1]; n2 = cells_[off + 3];
1636  if (n1 != n2)
1637  {
1638  PEdge e(n1,n2);
1639  array[i++] = (static_cast<T>((*(edge_table_.find(e))).second));
1640  }
1641  n1 = cells_[off + 2]; n2 = cells_[off + 3];
1642  if (n1 != n2)
1643  {
1644  PEdge e(n1,n2);
1645  array[i++] = (static_cast<T>((*(edge_table_.find(e))).second));
1646  }
1647  }
1648 
1649  template<class ARRAY, class INDEX>
1650  inline void get_edges_from_elem(ARRAY& array, INDEX idx) const
1651  {
1652  get_edges_from_cell(array,idx);
1653  }
1654 
1655  template<class ARRAY, class INDEX>
1656  inline void get_faces_from_cell(ARRAY& array, INDEX idx) const
1657  {
1658  ASSERTMSG(synchronized_ & Mesh::FACES_E,
1659  "TetVolMesh: Must call synchronize FACES_E first");
1660 
1661  array.clear();
1662  array.resize(4);
1663 
1664  const index_type off = idx << 2;
1665  const typename Node::index_type n0 = cells_[off ];
1666  const typename Node::index_type n1 = cells_[off + 1];
1667  const typename Node::index_type n2 = cells_[off + 2];
1668  const typename Node::index_type n3 = cells_[off + 3];
1669 
1670  PFaceNode f0(n3, n2, n1);
1671  PFaceNode f1(n0, n2, n3);
1672  PFaceNode f2(n3, n1, n0);
1673  PFaceNode f3(n0, n1, n2);
1674 
1675  array[0] = static_cast<typename ARRAY::value_type>(
1676  (*(face_table_.find(f0))).second);
1677  array[1] = static_cast<typename ARRAY::value_type>(
1678  (*(face_table_.find(f1))).second);
1679  array[2] = static_cast<typename ARRAY::value_type>(
1680  (*(face_table_.find(f2))).second);
1681  array[3] = static_cast<typename ARRAY::value_type>(
1682  (*(face_table_.find(f3))).second);
1683  }
1684 
1685  template<class ARRAY, class INDEX>
1686  inline void get_cells_from_node(ARRAY& array, INDEX idx) const
1687  {
1688  ASSERTMSG(synchronized_ & Mesh::NODE_NEIGHBORS_E,
1689  "TetVolMesh: Must call synchronize NODE_NEIGHBORS_E 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>(
1694  (node_neighbors_[idx][i])>>2);
1695  }
1696 
1697  template<class ARRAY, class INDEX>
1698  inline void get_edges_from_node(ARRAY& array, INDEX idx) const
1699  {
1700  ASSERTMSG(synchronized_ & (Mesh::NODE_NEIGHBORS_E),
1701  "HexVolMesh: Must call synchronize NODE_NEIGHBORS_E first");
1702  ASSERTMSG(synchronized_ & (Mesh::EDGES_E),
1703  "HexVolMesh: Must call synchronize EDGES_E first");
1704 
1705  // Get all the nodes that share an edge with this node
1706  const std::vector<typename Cell::index_type>& neighbors = node_neighbors_[idx];
1707 
1708  array.clear();
1709  array.reserve(neighbors.size());
1710  // Iterate through all those edges
1711  for (size_t n = 0; n < neighbors.size(); n++)
1712  {
1713  index_type cell_index = neighbors[n]&(~0x3);
1714  index_type node_index = neighbors[n]&0x3;
1715 
1716  const int *offset = TetVolEdgePerNodeTable[node_index];
1717 
1718  PEdgeNode e(cells_[cell_index+offset[0]],cells_[cell_index+offset[1]]);
1719  typename edge_nt::const_iterator iter = edge_table_.find(e);
1720  if (((edges_[iter->second].cells_[0])&(~0x7))==(cell_index<<1) )
1721  array.push_back(typename ARRAY::value_type(iter->second));
1722 
1723  PEdgeNode e1(cells_[cell_index+offset[2]],cells_[cell_index+offset[3]]);
1724  iter = edge_table_.find(e1);
1725  if (((edges_[iter->second].cells_[0])&(~0x7))==(cell_index<<1) )
1726  array.push_back(typename ARRAY::value_type(iter->second));
1727 
1728  PEdgeNode e2(cells_[cell_index+offset[4]],cells_[cell_index+offset[5]]);
1729  iter = edge_table_.find(e2);
1730  if (((edges_[iter->second].cells_[0])&(~0x7))==(cell_index<<1) )
1731  array.push_back(typename ARRAY::value_type(iter->second));
1732  }
1733  }
1734 
1735  template<class ARRAY, class INDEX>
1736  inline void get_faces_from_edge(ARRAY& array, INDEX idx) const
1737  {
1738  ASSERTMSG(synchronized_ & (Mesh::EDGES_E),
1739  "TetVolMesh: Must call synchronize EDGES_E first");
1740  ASSERTMSG(synchronized_ & (Mesh::FACES_E),
1741  "TetVolMesh: Must call synchronize FACES_E first");
1742 
1743  array.clear();
1744 
1745  for (size_t c=0; c<edges_[idx].cells_.size();c++)
1746  {
1747  index_type cell_index = ((edges_[idx].cells_[c])>>3)<<2;
1748  index_type face_index = (edges_[idx].cells_[c])&0x7;
1749 
1750  const int* off = TetVolFacePerEdgeTable[face_index];
1751 
1752  typename Node::index_type n1, n2, n3;
1753  typename face_nt::const_iterator fiter;
1754 
1755  n1 = cells_[cell_index+off[0]]; n2 = cells_[cell_index+off[1]];
1756  n3 = cells_[cell_index+off[2]];
1757 
1758  fiter = face_table_.find(PFaceNode(n1,n2,n3));
1759  if (((faces_[fiter->second].cells_[0])&(~0x3)) == cell_index)
1760  array.push_back(typename ARRAY::value_type(fiter->second));
1761 
1762  n1 = cells_[cell_index+off[3]]; n2 = cells_[cell_index+off[4]];
1763  n3 = cells_[cell_index+off[5]];
1764 
1765  fiter = face_table_.find(PFaceNode(n1,n2,n3));
1766  if (((faces_[fiter->second].cells_[0])&(~0x3)) == cell_index)
1767  array.push_back(typename ARRAY::value_type(fiter->second));
1768  }
1769  }
1770 
1771  template<class ARRAY, class INDEX>
1772  inline void get_faces_from_node(ARRAY& array, INDEX idx) const
1773  {
1774  ASSERTMSG(synchronized_ & (Mesh::NODE_NEIGHBORS_E),
1775  "TetVolMesh: Must call synchronize NODE_NEIGHBORS_E first");
1776  ASSERTMSG(synchronized_ & (Mesh::EDGES_E),
1777  "TetVolMesh: Must call synchronize EDGES_E first");
1778  ASSERTMSG(synchronized_ & (Mesh::FACES_E),
1779  "TetVolMesh: Must call synchronize FACES_E first");
1780 
1781  // Get all the nodes that share an edge with this node
1782  const std::vector<typename Cell::index_type>& neighbors = node_neighbors_[idx];
1783 
1784  array.clear();
1785  array.reserve(neighbors.size());
1786  // Iterate through all those edges
1787  for (size_t n = 0; n < neighbors.size(); n++)
1788  {
1789  index_type cell_index = neighbors[n]&(~0x3);
1790  index_type node_index = neighbors[n]&0x3;
1791 
1792  const int *offset = TetVolFacePerNodeTable[node_index];
1793 
1794  PFaceNode e(cells_[cell_index+offset[0]],
1795  cells_[cell_index+offset[1]],cells_[cell_index+offset[2]]);
1796 
1797  typename face_nt::const_iterator iter = face_table_.find(e);
1798  if (((faces_[iter->second].cells_[0])&(~0x3))==cell_index)
1799  array.push_back(typename ARRAY::value_type(iter->second));
1800 
1801  PFaceNode e1(cells_[cell_index+offset[3]],cells_[cell_index+offset[4]],
1802  cells_[cell_index+offset[5]]);
1803  iter = face_table_.find(e1);
1804  if (((faces_[iter->second].cells_[0])&(~0x3))==cell_index )
1805  array.push_back(typename ARRAY::value_type(iter->second));
1806 
1807  PFaceNode e2(cells_[cell_index+offset[6]],cells_[cell_index+offset[7]],
1808  cells_[cell_index+offset[8]]);
1809  iter = face_table_.find(e2);
1810  if (((faces_[iter->second].cells_[0])&(~0x3))==cell_index )
1811  array.push_back(typename ARRAY::value_type(iter->second));
1812  }
1813  }
1814 
1815  template<class ARRAY, class INDEX>
1816  inline void get_cells_from_edge(ARRAY& array, INDEX idx) const
1817  {
1818  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
1819  "TetVolMesh: Must call synchronize EDGES_E first");
1820  for (size_t i=0; i< edges_[idx].cells_.size(); i++)
1821  array[i] = static_cast<typename ARRAY::value_type>(edges_[idx].cells_[i]);
1822  }
1823 
1824  template<class ARRAY, class INDEX>
1825  inline void get_cells_from_face(ARRAY& array, INDEX idx) const
1826  {
1827  ASSERTMSG(synchronized_ & Mesh::FACES_E,
1828  "TetVolMesh: Must call synchronize FACES_E first");
1829  if (faces_[idx].cells_[1] == MESH_NO_NEIGHBOR)
1830  {
1831  array.resize(1);
1832  array[0] = static_cast<typename ARRAY::value_type>((faces_[idx].cells_[0])>>2);
1833  }
1834  else
1835  {
1836  array.resize(2);
1837  array[0] = static_cast<typename ARRAY::value_type>((faces_[idx].cells_[0])>>2);
1838  array[1] = static_cast<typename ARRAY::value_type>((faces_[idx].cells_[1])>>2);
1839  }
1840  }
1841 
1842  template <class ARRAY, class INDEX>
1843  inline void set_nodes_by_elem(ARRAY &array, INDEX idx)
1844  {
1845  for (index_type n = 0; n < 4; ++n)
1846  cells_[idx * 4 + n] = static_cast<index_type>(array[n]);
1847  }
1848 
1849  template <class INDEX1, class INDEX2>
1850  inline bool get_elem_neighbor(INDEX1 &neighbor, INDEX1 elem, INDEX2 delem) const
1851  {
1852  ASSERTMSG(synchronized_ & Mesh::FACES_E,
1853  "TetVolMesh: Must call synchronize FACES_E on TetVolMesh first");
1854  const PFaceCell &f = faces_[delem];
1855 
1856  if (elem == static_cast<INDEX1>((f.cells_[0])>>2))
1857  {
1858  if (f.cells_[1] == MESH_NO_NEIGHBOR) return (false);
1859  neighbor = static_cast<INDEX1>((f.cells_[1])>>2);
1860  }
1861  else
1862  {
1863  if (f.cells_[0] == MESH_NO_NEIGHBOR) return (false);
1864  neighbor = static_cast<INDEX1>((f.cells_[0])>>2);
1865  }
1866  return (true);
1867  }
1868 
1869  template <class ARRAY, class INDEX1, class INDEX2>
1870  inline bool get_elem_neighbors(ARRAY &array, INDEX1 elem, INDEX2 delem) const
1871  {
1872  ASSERTMSG(synchronized_ & Mesh::FACES_E,
1873  "Must call synchronize FACES_E on TetVolMesh first");
1874 
1875  const PFaceCell &f = faces_[delem];
1876 
1877  if (elem == static_cast<INDEX1>((f.cells_[0])>>2))
1878  {
1879  if (f.cells_[1] == MESH_NO_NEIGHBOR) return (false);
1880  array.resize(1);
1881  array[0] = static_cast<typename ARRAY::value_type>((f.cells_[1])>>2);
1882  }
1883  else
1884  {
1885  if (f.cells_[0] == MESH_NO_NEIGHBOR) return (false);
1886  array.resize(1);
1887  array[0] = static_cast<typename ARRAY::value_type>((f.cells_[0])>>2);
1888  }
1889  return (true);
1890  }
1891 
1892  template <class ARRAY, class INDEX>
1893  inline void get_elem_neighbors(ARRAY &array, INDEX idx) const
1894  {
1895  typename Face::array_type faces;
1896  get_faces_from_cell(faces, idx);
1897 
1898  array.clear();
1899  array.reserve(faces.size());
1900  typename Face::array_type::iterator iter = faces.begin();
1901  while(iter != faces.end())
1902  {
1903  typename ARRAY::value_type nbor;
1904  if (get_elem_neighbor(nbor, idx, *iter))
1905  {
1906  array.push_back(nbor);
1907  }
1908  ++iter;
1909  }
1910  }
1911 
1912  /// We should optimize this function more
1913  template <class ARRAY, class INDEX>
1914  inline void get_node_neighbors(ARRAY &array, INDEX node) const
1915  {
1916  ASSERTMSG(synchronized_ & Mesh::NODE_NEIGHBORS_E,
1917  "Must call synchronize NODE_NEIGHBORS_E on TetVolMesh first.");
1918  size_t sz = node_neighbors_[node].size();
1919 
1920  std::set<index_type> inserted;
1921  for (size_t i = 0; i < sz; i++)
1922  {
1923  const index_type base = ((node_neighbors_[node][i])&(~0x3));
1924  for (index_type c = base; c < base+4; ++c)
1925  {
1926  if (cells_[c] != node) inserted.insert(cells_[c]);
1927  }
1928  }
1929 
1930  array.clear();
1931  array.reserve(inserted.size());
1932  array.insert(array.begin(), inserted.begin(), inserted.end());
1933  }
1934 
1935  template <class INDEX>
1936  inline bool locate_node(INDEX &node, const Core::Geometry::Point &p) const
1937  {
1938  typename Node::size_type sz; size(sz);
1939 
1940  /// If there are no nodes we cannot find a closest point
1941  if (sz == 0) return (false);
1942 
1943  /// Check first guess
1944  if (node >= 0 && node < sz)
1945  {
1946  if ((p - points_[node]).length2() < epsilon2_) return (true);
1947  }
1948 
1949  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
1950  "TetVolMesh::locate_node requires synchronize(NODE_LOCATE_E).")
1951 
1952  // get grid sizes
1953  const size_type ni = node_grid_->get_ni()-1;
1954  const size_type nj = node_grid_->get_nj()-1;
1955  const size_type nk = node_grid_->get_nk()-1;
1956 
1957  // Convert to grid coordinates.
1958  index_type bi, bj, bk, ei, ej, ek;
1959  node_grid_->unsafe_locate(bi, bj, bk, p);
1960 
1961  // Clamp to l;..;,closest point on the grid.
1962  if (bi > ni) bi =ni; if (bi < 0) bi = 0;
1963  if (bj > nj) bj =nj; if (bj < 0) bj = 0;
1964  if (bk > nk) bk =nk; if (bk < 0) bk = 0;
1965 
1966  ei = bi;
1967  ej = bj;
1968  ek = bk;
1969 
1970  double dmin = DBL_MAX;
1971  bool found;
1972  do
1973  {
1974  found = true;
1975  /// We need to do a full shell without any elements that are closer
1976  /// to make sure there no closer elements in neighboring searchgrid cells
1977 
1978  for (index_type i = bi; i <= ei; i++)
1979  {
1980  if (i < 0 || i > ni) continue;
1981  for (index_type j = bj; j <= ej; j++)
1982  {
1983  if (j < 0 || j > nj) continue;
1984  for (index_type k = bk; k <= ek; k++)
1985  {
1986  if (k < 0 || k > nk) continue;
1987  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
1988  {
1989  if (node_grid_->min_distance_squared(p, i, j, k) < dmin)
1990  {
1991  found = false;
1992  SearchGridT<index_type>::iterator it, eit;
1993  node_grid_->lookup_ijk(it, eit, i, j, k);
1994 
1995  while (it != eit)
1996  {
1997  const Core::Geometry::Point point = points_[*it];
1998  const double dist = (p-point).length2();
1999 
2000  if (dist < dmin)
2001  {
2002  node = INDEX(*it);
2003  dmin = dist;
2004 
2005  if (dist < epsilon2_) return (true);
2006  }
2007  ++it;
2008  }
2009  }
2010  }
2011  }
2012  }
2013  }
2014  bi--;ei++;
2015  bj--;ej++;
2016  bk--;ek++;
2017  }
2018  while ((!found)||(dmin == DBL_MAX)) ;
2019 
2020  return (true);
2021  }
2022 
2023  /// @todo: Fix this function, needs to use search grid
2024  template <class INDEX>
2025  inline bool locate_edge(INDEX &edge, const Core::Geometry::Point &p) const
2026  {
2027  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
2028  "Must call synchronize EDGES_E on TetVolMesh first");
2029 
2030  typename Edge::iterator bi, ei;
2031  typename Node::array_type nodes;
2032  begin(bi);
2033  end(ei);
2034 
2035  bool found = false;
2036  double mindist = 0.0;
2037  while (bi != ei)
2038  {
2039  get_nodes(nodes,*bi);
2040  const double dist = distance_to_line2(p, points_[nodes[0]],
2041  points_[nodes[1]],epsilon_);
2042  if (!found || dist < mindist)
2043  {
2044  edge = static_cast<INDEX>(*bi);
2045  mindist = dist;
2046  found = true;
2047  }
2048  ++bi;
2049  }
2050  return (found);
2051  }
2052 
2053  /// @todo: Fix this function, needs to use search grid
2054  template <class INDEX>
2055  inline bool locate_face(INDEX &face, const Core::Geometry::Point &p) const
2056  {
2057  ASSERTMSG(synchronized_ & Mesh::FACES_E,
2058  "Must call synchronize FACES_E on TetVolMesh first");
2059 
2060  bool found = false;
2061  double mindist = DBL_MAX;
2062  typename Face::iterator bi; begin(bi);
2063  typename Face::iterator ei; end(ei);
2064  while (bi != ei)
2065  {
2066  Core::Geometry::Point c;
2067  get_center(c, *bi);
2068  const double dist = (p - c).length2();
2069  if (!found || dist < mindist)
2070  {
2071  mindist = dist;
2072  face = static_cast<INDEX>(*bi);
2073  found = true;
2074  }
2075  ++bi;
2076  }
2077  return (found);
2078  }
2079 
2080  template <class INDEX>
2081  inline bool locate_elem(INDEX &elem, const Core::Geometry::Point &p) const
2082  {
2083  /// @todo: Generate bounding boxes for elements and integrate this into the
2084  // basis function code.
2085  if (basis_.polynomial_order() > 1) return elem_locate(elem, *this, p);
2086 
2087  typename Elem::size_type sz; size(sz);
2088 
2089  /// If there are no nodes we cannot find a closest point
2090  if (sz == 0) return (false);
2091 
2092  /// Check whether the estimate given in idx is the point we are looking for
2093  if ((elem > 0)&&(elem < sz))
2094  {
2095  if (inside(elem,p)) return (true);
2096  }
2097 
2098  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
2099  "TetVolMesh: need to synchronize ELEM_LOCATE_E first");
2100 
2101  typename SearchGridT<index_type>::iterator it, eit;
2102  if (elem_grid_->lookup(it, eit, p))
2103  {
2104  while (it != eit)
2105  {
2106  if (inside(typename Elem::index_type(*it), p))
2107  {
2108  elem = static_cast<INDEX>(*it);
2109  return (true);
2110  }
2111  ++it;
2112  }
2113  }
2114  return (false);
2115  }
2116 
2117 
2118  template <class ARRAY>
2119  inline bool locate_elems(ARRAY &array, const Core::Geometry::BBox &b) const
2120  {
2121 
2122  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
2123  "TetVolMesh::locate_elems requires synchronize(ELEM_LOCATE_E).")
2124 
2125  array.clear();
2126  index_type is,js,ks;
2127  index_type ie,je,ke;
2128  elem_grid_->locate_clamp(is,js,ks,b.min());
2129  elem_grid_->locate_clamp(ie,je,ke,b.max());
2130  for (index_type i=is; i<=ie;i++)
2131  for (index_type j=js; j<je;j++)
2132  for (index_type k=ks; k<ke; k++)
2133  {
2134  typename SearchGridT<index_type>::iterator it, eit;
2135  elem_grid_->lookup_ijk(it, eit, i, j, k);
2136  while (it != eit)
2137  {
2138  size_t p=0;
2139  for (;p<array.size();p++) if (array[p] == typename ARRAY::value_type(*it)) break;
2140  if (p == array.size()) array.push_back(typename ARRAY::value_type(*it));
2141  ++it;
2142  }
2143  }
2144 
2145  return (array.size() > 0);
2146  }
2147 
2148  template <class INDEX, class ARRAY>
2149  inline bool locate_elem(INDEX &elem, ARRAY& coords, const Core::Geometry::Point &p) const
2150  {
2151  /// @todo: Generate bounding boxes for elements and integrate this into the
2152  // basis function code.
2153  if (basis_.polynomial_order() > 1) return elem_locate(elem, *this, p);
2154 
2155  typename Elem::size_type sz; size(sz);
2156 
2157  /// If there are no nodes we cannot find a closest point
2158  if (sz == 0) return (false);
2159 
2160  /// Check whether the estimate given in idx is the point we are looking for
2161  if ((elem > 0)&&(elem < sz))
2162  {
2163  if (inside(elem,p))
2164  {
2165  ElemData ed(*this, elem);
2166  basis_.get_coords(coords, p, ed);
2167  return (true);
2168  }
2169  }
2170 
2171  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
2172  "TetVolMesh: need to synchronize ELEM_LOCATE_E first");
2173 
2174  typename SearchGridT<index_type>::iterator it, eit;
2175  if (elem_grid_->lookup(it, eit, p))
2176  {
2177  while (it != eit)
2178  {
2179  if (inside(typename Elem::index_type(*it), p))
2180  {
2181  elem = static_cast<INDEX>(*it);
2182  ElemData ed(*this, elem);
2183  basis_.get_coords(coords, p, ed);
2184  return (true);
2185  }
2186  ++it;
2187  }
2188  }
2189 
2190  return (false);
2191  }
2192 
2193  template <class INDEX>
2194  inline void get_node_center(Core::Geometry::Point &p, INDEX idx) const
2195  {
2196  p = points_[idx];
2197  }
2198 
2199  template <class INDEX>
2200  inline void get_edge_center(Core::Geometry::Point &p, INDEX idx) const
2201  {
2202  typename Node::array_type arr;
2203  get_nodes_from_edge(arr, idx);
2204  Core::Geometry::Point p1;
2205  get_node_center(p, arr[0]);
2206  get_node_center(p1, arr[1]);
2207  p += p1;
2208  p *= 0.5;
2209  }
2210 
2211  template <class INDEX>
2212  inline void get_face_center(Core::Geometry::Point &p, INDEX idx) const
2213  {
2214  const double s = 1.0/3.0;
2215  typename Node::array_type arr;
2216  get_nodes_from_face(arr, idx);
2217  p = points_[arr[0]];
2218  const Core::Geometry::Point &p1 = points_[arr[1]];
2219  const Core::Geometry::Point &p2 = points_[arr[2]];
2220 
2221  p += p1;
2222  p += p2;
2223  p *= s;
2224  }
2225 
2226  template <class INDEX>
2227  inline void get_cell_center(Core::Geometry::Point &p, INDEX idx) const
2228  {
2229  const double s = 0.25;
2230  const Core::Geometry::Point &p0 = points_[cells_[idx * 4 + 0]];
2231  const Core::Geometry::Point &p1 = points_[cells_[idx * 4 + 1]];
2232  const Core::Geometry::Point &p2 = points_[cells_[idx * 4 + 2]];
2233  const Core::Geometry::Point &p3 = points_[cells_[idx * 4 + 3]];
2234 
2235  p = Core::Geometry::Point((p0 + p1 + p2 + p3) * s);
2236  }
2237 
2238  //////////////////////////////////////////////////////////////
2239  // Internal functions that the mesh depends on
2240 
2241 protected:
2242 
2243  void compute_node_neighbors();
2244  void compute_edges();
2245  void compute_faces();
2246  void compute_node_grid();
2247  void compute_elem_grid();
2248  void compute_bounding_box();
2249 
2250  void insert_elem_into_grid(typename Elem::index_type ci);
2251  void remove_elem_from_grid(typename Elem::index_type ci);
2252  void insert_node_into_grid(typename Node::index_type ci);
2253  void remove_node_from_grid(typename Node::index_type ci);
2254 
2255  const Core::Geometry::Point &point(typename Node::index_type i) { return points_[i]; }
2256 
2257  template<class INDEX>
2258  bool inside(INDEX idx, const Core::Geometry::Point &p) const
2259  {
2260  const Core::Geometry::Point &p0 = points_[cells_[idx*4+0]];
2261  const Core::Geometry::Point &p1 = points_[cells_[idx*4+1]];
2262  const Core::Geometry::Point &p2 = points_[cells_[idx*4+2]];
2263  const Core::Geometry::Point &p3 = points_[cells_[idx*4+3]];
2264 
2265  const double x0 = p0.x();
2266  const double y0 = p0.y();
2267  const double z0 = p0.z();
2268  const double x1 = p1.x();
2269  const double y1 = p1.y();
2270  const double z1 = p1.z();
2271  const double x2 = p2.x();
2272  const double y2 = p2.y();
2273  const double z2 = p2.z();
2274  const double x3 = p3.x();
2275  const double y3 = p3.y();
2276  const double z3 = p3.z();
2277 
2278  const double a0 = + x1*(y2*z3-y3*z2) + x2*(y3*z1-y1*z3) + x3*(y1*z2-y2*z1);
2279  const double a1 = - x2*(y3*z0-y0*z3) - x3*(y0*z2-y2*z0) - x0*(y2*z3-y3*z2);
2280  const double a2 = + x3*(y0*z1-y1*z0) + x0*(y1*z3-y3*z1) + x1*(y3*z0-y0*z3);
2281  const double a3 = - x0*(y1*z2-y2*z1) - x1*(y2*z0-y0*z2) - x2*(y0*z1-y1*z0);
2282  const double iV6 = 1.0 / (a0+a1+a2+a3);
2283 
2284  const double b0 = - (y2*z3-y3*z2) - (y3*z1-y1*z3) - (y1*z2-y2*z1);
2285  const double c0 = + (x2*z3-x3*z2) + (x3*z1-x1*z3) + (x1*z2-x2*z1);
2286  const double d0 = - (x2*y3-x3*y2) - (x3*y1-x1*y3) - (x1*y2-x2*y1);
2287  const double s0 = iV6 * (a0 + b0*p.x() + c0*p.y() + d0*p.z());
2288  if (s0 < -1e-7) return (false);
2289 
2290  const double b1 = + (y3*z0-y0*z3) + (y0*z2-y2*z0) + (y2*z3-y3*z2);
2291  const double c1 = - (x3*z0-x0*z3) - (x0*z2-x2*z0) - (x2*z3-x3*z2);
2292  const double d1 = + (x3*y0-x0*y3) + (x0*y2-x2*y0) + (x2*y3-x3*y2);
2293  const double s1 = iV6 * (a1 + b1*p.x() + c1*p.y() + d1*p.z());
2294  if (s1 < -1e-7) return (false);
2295 
2296  const double b2 = - (y0*z1-y1*z0) - (y1*z3-y3*z1) - (y3*z0-y0*z3);
2297  const double c2 = + (x0*z1-x1*z0) + (x1*z3-x3*z1) + (x3*z0-x0*z3);
2298  const double d2 = - (x0*y1-x1*y0) - (x1*y3-x3*y1) - (x3*y0-x0*y3);
2299  const double s2 = iV6 * (a2 + b2*p.x() + c2*p.y() + d2*p.z());
2300  if (s2 < -1e-7) return (false);
2301 
2302  const double b3 = +(y1*z2-y2*z1) + (y2*z0-y0*z2) + (y0*z1-y1*z0);
2303  const double c3 = -(x1*z2-x2*z1) - (x2*z0-x0*z2) - (x0*z1-x1*z0);
2304  const double d3 = +(x1*y2-x2*y1) + (x2*y0-x0*y2) + (x0*y1-x1*y0);
2305  const double s3 = iV6 * (a3 + b3*p.x() + c3*p.y() + d3*p.z());
2306  if (s3 < -1e-7) return (false);
2307 
2308  return (true);
2309  }
2310 
2311  /// all the nodes.
2312  std::vector<Core::Geometry::Point> points_;
2313 
2314  /// each 4 indicies make up a tet
2315  std::vector<under_type> cells_;
2316 
2317  /// Face information.
2318  class PFaceCell {
2319  public:
2320  PFaceCell()
2321  {
2322  cells_[0] = MESH_NO_NEIGHBOR;
2323  cells_[1] = MESH_NO_NEIGHBOR;
2324  }
2325 
2326  bool shared() const { return ((cells_[0] != MESH_NO_NEIGHBOR) &&
2327  (cells_[1] != MESH_NO_NEIGHBOR)); }
2328  index_type cells_[2];
2329  };
2330 
2331  class PFaceNode {
2332  public:
2333  // The order of nodes_ corresponds with cells_[0] for CW/CCW purposes.
2334  typename Node::index_type nodes_[3]; /// 3 nodes makes a face.
2335 
2336  PFaceNode() {
2337  nodes_[0] = MESH_NO_NEIGHBOR;
2338  nodes_[1] = MESH_NO_NEIGHBOR;
2339  nodes_[2] = MESH_NO_NEIGHBOR;
2340  }
2341 
2342  // snodes_ must be sorted. See Hash Function below.
2343  PFaceNode(typename Node::index_type n1, typename Node::index_type n2,
2344  typename Node::index_type n3)
2345  {
2346  if (n2 < n1 && n2 < n3)
2347  {
2348  typename Node::index_type t;
2349  if (n3 < n1)
2350  {
2351  t = n1; n1 = n2; n2 = n3; n3 = t;
2352  }
2353  else
2354  {
2355  t = n1; n1 = n2; n2 = t;
2356  }
2357  }
2358  else if (n3 < n1 && n3 < n2)
2359  {
2360  typename Node::index_type t;
2361  if (n1 < n2)
2362  {
2363  t = n1; n1 = n3; n3 = n2; n2 = t;
2364  }
2365  else
2366  {
2367  t = n1; n1 = n3; n3 = t;
2368  }
2369  }
2370  else if (n3 < n2)
2371  {
2372  typename Node::index_type t;
2373  t = n2; n2 = n3; n3 = t;
2374  }
2375  nodes_[0] = n1;
2376  nodes_[1] = n2;
2377  nodes_[2] = n3;
2378  }
2379 
2380  /// true if both have the same nodes (order does not matter)
2381  bool operator==(const PFaceNode &f) const
2382  {
2383  return ((nodes_[0] == f.nodes_[0]) && (nodes_[1] == f.nodes_[1]) &&
2384  (nodes_[2] == f.nodes_[2]));
2385  }
2386 
2387  /// Compares each node. When a non equal node is found the <
2388  /// operator is applied.
2389  bool operator<(const PFaceNode &f) const
2390  {
2391  if (nodes_[0] == f.nodes_[0])
2392  if (nodes_[1] == f.nodes_[1])
2393  return (nodes_[2] < f.nodes_[2]);
2394  else
2395  return (nodes_[1] < f.nodes_[1]);
2396  else
2397  return (nodes_[0] < f.nodes_[0]);
2398  }
2399  };
2400 
2401  class PFace : public PFaceNode, public PFaceCell
2402  {
2403  public:
2404  PFace(typename Node::index_type n1, typename Node::index_type n2,
2405  typename Node::index_type n3) :
2406  PFaceNode(n1,n2,n3) {}
2407  };
2408 
2409  /// Edge information.
2410  class PEdgeNode {
2411  public:
2412  typename Node::index_type nodes_[2];
2413 
2414  PEdgeNode()
2415  {
2416  nodes_[0] = MESH_NO_NEIGHBOR;
2417  nodes_[1] = MESH_NO_NEIGHBOR;
2418  }
2419  // node_[0] must be smaller than node_[1]. See Hash Function below.
2420  PEdgeNode(typename Node::index_type n1,
2421  typename Node::index_type n2)
2422  {
2423  if (n1 < n2)
2424  {
2425  nodes_[0] = n1;
2426  nodes_[1] = n2;
2427  }
2428  else
2429  {
2430  nodes_[0] = n2;
2431  nodes_[1] = n1;
2432  }
2433  }
2434 
2435  /// true if both have the same nodes (order does not matter)
2436  bool operator==(const PEdgeNode &e) const
2437  {
2438  return ((nodes_[0] == e.nodes_[0]) && (nodes_[1] == e.nodes_[1]));
2439  }
2440 
2441  /// Compares each node. When a non equal node is found the <
2442  /// operator is applied.
2443  bool operator<(const PEdgeNode &e) const
2444  {
2445  if (nodes_[0] == e.nodes_[0])
2446  return (nodes_[1] < e.nodes_[1]);
2447  else
2448  return (nodes_[0] < e.nodes_[0]);
2449  }
2450  };
2451 
2452  class PEdgeCell {
2453  public:
2454  /// list of all the cells this edge is in.
2455  std::vector<index_type> cells_;
2456 
2457  bool shared() const { return cells_.size() > 1; }
2458  };
2459 
2460  /// Edge information.
2461  class PEdge : public PEdgeNode, public PEdgeCell
2462  {
2463  public:
2464  PEdge(typename Node::index_type n1,typename Node::index_type n2) :
2465  PEdgeNode(n1,n2) {}
2466  };
2467 
2468  /// hash the egde's node_indecies such that edges with the same nodes
2469  /// hash to the same value. nodes are sorted on edge construction.
2470  static const int sz_int = sizeof(int) * 8; // in bits
2471  struct FaceHash
2472  {
2473  // These are needed by the hash_map particularly
2474  // ANSI C++ allows us to initialize these variables in the
2475  // declaration. However there may be compilers which will complain
2476  // about it.
2477  static const size_t bucket_size = 4;
2478  static const size_t min_buckets = 8;
2479 
2480  /// These are for our own use (making the hash function).
2481  static const int sz_third_int = (int)(sz_int / 3);
2482  static const int up_mask = (~((int)0) << sz_third_int << sz_third_int);
2483  static const int mid_mask = up_mask ^ (~((int)0) << sz_third_int);
2484  static const int low_mask = ~(up_mask | mid_mask);
2485 
2486  /// This is the hash function
2487  template <class PFACE>
2488  size_t operator()(const PFACE &f) const
2489  {
2490  return ((up_mask & (static_cast<int>(f.nodes_[0]) << sz_third_int << sz_third_int)) |
2491  (mid_mask & (static_cast<int>(f.nodes_[1]) << sz_third_int)) |
2492  (low_mask & static_cast<int>(f.nodes_[2])));
2493  }
2494  /// This should return less than rather than equal to.
2495 
2496  template <class PFACE>
2497  bool operator()(const PFACE &f1, const PFACE& f2) const
2498  {
2499  return f1 < f2;
2500  }
2501  };
2502 
2503  /// hash the egde's node_indecies such that edges with the same nodes
2504  /// hash to the same value. nodes are sorted on edge construction.
2505  struct EdgeHash {
2506  /// These are needed by the hash_map particularly
2507  // ANSI C++ allows us to initialize these variables in the
2508  // declaration. However there may be compilers which will complain
2509  // about it.
2510  static const size_t bucket_size = 4;
2511  static const size_t min_buckets = 8;
2512 
2513  /// These are for our own use (making the hash function.
2514  static const int sz_int = sizeof(int) * 8; // in bits
2515  static const int sz_half_int = sizeof(int) << 2; // in bits
2516  static const int up_mask = ((~((int)0)) << sz_half_int);
2517  static const int low_mask = (~((int)0) ^ up_mask);
2518 
2519  /// This is the hash function
2520  template <class PEDGE>
2521  size_t operator()(const PEDGE &e) const
2522  {
2523  return (static_cast<int>(e.nodes_[0]) << sz_half_int) |
2524  (low_mask & static_cast<int>(e.nodes_[1]));
2525  }
2526 
2527  /// This should return less than rather than equal to.
2528  template <class PEDGE>
2529  bool operator()(const PEDGE &e1, const PEDGE& e2) const
2530  {
2531  return e1 < e2;
2532  }
2533  };
2534 
2535 /// Define the hash_map type, as this is not yet an approved type under Windows
2536 /// it is located in stdext
2537 
2538 #ifdef HAVE_HASH_MAP
2539  #if defined(_WIN32)
2540  /// hash_map is in stdext namespace
2541  typedef stdext::hash_map<PFace, typename Face::index_type, FaceHash> face_ht;
2542  typedef stdext::hash_map<PFaceNode, typename Face::index_type, FaceHash> face_nt;
2543  typedef stdext::hash_map<PEdge, typename Edge::index_type, EdgeHash> edge_ht;
2544  typedef stdext::hash_map<PEdgeNode, typename Edge::index_type, EdgeHash> edge_nt;
2545  #else
2546  /// hash_map is in stdext namespace
2547  typedef hash_map<PFace, typename Face::index_type, FaceHash> face_ht;
2548  typedef hash_map<PFaceNode, typename Face::index_type, FaceHash> face_nt;
2549  typedef hash_map<PEdge, typename Edge::index_type, EdgeHash> edge_ht;
2550  typedef hash_map<PEdgeNode, typename Edge::index_type, EdgeHash> edge_nt;
2551  #endif
2552 #else
2553  typedef std::map<PFace, typename Face::index_type, FaceHash> face_ht;
2554  typedef std::map<PFaceNode, typename Face::index_type, FaceHash> face_nt;
2555  typedef std::map<PEdge, typename Edge::index_type, EdgeHash> edge_ht;
2556  typedef std::map<PEdgeNode, typename Edge::index_type, EdgeHash> edge_nt;
2557 #endif
2558 
2559  typedef std::vector<PFaceCell> face_ct;
2560  typedef std::vector<PEdgeCell> edge_ct;
2561 
2562  // These should not be called outside of the synchronize_lock_.
2563 
2564  /// container for face storage. Must be computed each time
2565  /// nodes or cells change.
2566  face_ct faces_;
2567  face_nt face_table_;
2568  /// container for edge storage. Must be computed each time
2569  /// nodes or cells change.
2570  edge_ct edges_;
2571  edge_nt edge_table_;
2572 
2573  inline void remove_edge(typename Node::index_type n1,
2574  typename Node::index_type n2,
2575  typename Cell::index_type ci,
2576  bool table_only = false);
2577 
2578  inline void hash_edge(typename Node::index_type n1,
2579  typename Node::index_type n2,
2580  index_type combined_index,
2581  edge_ht& table);
2582 
2583  inline void add_edge(typename Node::index_type n1,
2584  typename Node::index_type n2,
2585  index_type combined_index);
2586 
2587  inline void remove_face(typename Node::index_type n1,
2588  typename Node::index_type n2,
2589  typename Node::index_type n3,
2590  typename Cell::index_type ci,
2591  bool table_only = false);
2592  inline void hash_face(typename Node::index_type n1,
2593  typename Node::index_type n2,
2594  typename Node::index_type n3,
2595  index_type ci, face_ht& table);
2596  inline void add_face(typename Node::index_type n1,
2597  typename Node::index_type n2,
2598  typename Node::index_type n3,
2599  index_type combined_index);
2600 
2601  std::vector<std::vector<typename Cell::index_type> > node_neighbors_;
2602  std::vector<unsigned char> boundary_faces_;
2603 
2604  /// This grid is used as an acceleration structure to expedite calls
2605  /// to locate. For each cell in the grid, we store a list of which
2606  /// tets overlap that grid cell -- to find the tet which contains a
2607  /// point, we simply find which grid cell contains that point, and
2608  /// then search just those tets that overlap that grid cell.
2609  boost::shared_ptr<SearchGridT<index_type> > node_grid_;
2610  boost::shared_ptr<SearchGridT<index_type> > elem_grid_;
2611 
2612  // Lock and Condition Variable for hand shaking
2613  mutable Core::Thread::Mutex synchronize_lock_;
2614  Core::Thread::ConditionVariable synchronize_cond_;
2615 
2616  // Which tables have been computed
2617  mask_type synchronized_;
2618  // Which tables are currently being computed
2619  mask_type synchronizing_;
2620 
2621  Core::Geometry::BBox bbox_;
2622  Basis basis_;
2623  double epsilon_;
2624  double epsilon2_;
2625  double epsilon3_;
2626 
2627  /// Pointer to virtual interface
2628  boost::shared_ptr<VMesh> vmesh_;
2629 
2630 public:
2631 
2632  /// Always creates 4 tets, does not handle element boundaries properly.
2633  bool insert_node_in_cell(typename Cell::array_type &tets,
2634  typename Cell::index_type ci,
2635  typename Node::index_type &ni,
2636  const Core::Geometry::Point &p);
2637 
2638  void insert_node_in_face(typename Cell::array_type &tets,
2639  typename Node::index_type ni,
2640  typename Cell::index_type ci,
2641  const PFaceNode &pf);
2642 
2643  void insert_node_in_edge(typename Cell::array_type &tets,
2644  typename Node::index_type ni,
2645  typename Cell::index_type ci,
2646  const PEdgeNode &e);
2647 
2648 }; // end class TetVolMesh
2649 
2650 template <class Basis>
2651 int
2652 TetVolMesh<Basis>::compute_checksum()
2653 {
2654  int sum = 0;
2655  sum += SCIRun::compute_checksum(&points_[0],points_.size());
2656  sum += SCIRun::compute_checksum(&cells_[0],cells_.size());
2657  return (sum);
2658 }
2659 
2660 template <class Basis>
2661 TetVolMesh<Basis>::TetVolMesh() :
2662  points_(0),
2663  cells_(0),
2664  faces_(0),
2665  face_table_(),
2666  edges_(0),
2667  edge_table_(),
2668  synchronize_lock_("TetVolMesh lock"),
2669  synchronize_cond_("TetVolMesh condition variable"),
2670  synchronized_(Mesh::NODES_E | Mesh::CELLS_E),
2671  synchronizing_(0),
2672  epsilon_(0.0),
2673  epsilon2_(0.0),
2674  epsilon3_(0.0)
2675 {
2676  DEBUG_CONSTRUCTOR("TetVolMesh")
2677 
2678  vmesh_.reset(CreateVTetVolMesh(this));
2679 }
2680 
2681 template <class Basis>
2682 TetVolMesh<Basis>::TetVolMesh(const TetVolMesh &copy) :
2683  Mesh(copy),
2684  points_(0),
2685  cells_(0),
2686  faces_(0),
2687  face_table_(),
2688  edges_(0),
2689  edge_table_(),
2690  synchronize_lock_("TetVolMesh lock"),
2691  synchronize_cond_("TetVolMesh condition variable"),
2692  synchronized_(Mesh::NODES_E | Mesh::CELLS_E),
2693  synchronizing_(0),
2694  epsilon_(0.0),
2695  epsilon2_(0.0),
2696  epsilon3_(0.0)
2697 {
2698  DEBUG_CONSTRUCTOR("TetVolMesh")
2699 
2700  /// We need to lock before we can copy these as these
2701  /// structures are generate dynamically when they are
2702  /// needed.
2703  copy.synchronize_lock_.lock();
2704 
2705  points_ = copy.points_;
2706  cells_ = copy.cells_;
2707 
2708  // Epsilon does not require much space, hence copy those
2709  synchronized_ |= copy.synchronized_ & Mesh::BOUNDING_BOX_E;
2710  bbox_ = copy.bbox_;
2711  epsilon_ = copy.epsilon_;
2712  epsilon2_ = copy.epsilon2_;
2713  epsilon3_ = copy.epsilon3_;
2714 
2715  copy.synchronize_lock_.unlock();
2716 
2717  /// Create a new virtual interface for this copy
2718  /// all pointers have changed hence create a new
2719  /// virtual interface class
2720  vmesh_.reset(CreateVTetVolMesh(this));
2721 }
2722 
2723 template <class Basis>
2724 TetVolMesh<Basis>::~TetVolMesh()
2725 {
2726  DEBUG_DESTRUCTOR("TetVolMesh")
2727 }
2728 
2729 template <class Basis>
2730 template <class Iter, class Functor>
2731 void
2732 TetVolMesh<Basis>::fill_points(Iter begin, Iter end, Functor fill_ftor)
2733 {
2734  synchronize_lock_.lock();
2735  Iter iter = begin;
2736  points_.resize(end - begin); // resize to the new size
2737  std::vector<Core::Geometry::Point>::iterator piter = points_.begin();
2738  while (iter != end)
2739  {
2740  *piter = fill_ftor(*iter);
2741  ++piter; ++iter;
2742  }
2743  synchronize_lock_.unlock();
2744 }
2745 
2746 template <class Basis>
2747 template <class Iter, class Functor>
2748 void
2749 TetVolMesh<Basis>::fill_cells(Iter begin, Iter end, Functor fill_ftor)
2750 {
2751  synchronize_lock_.lock();
2752  Iter iter = begin;
2753  cells_.resize((end - begin) * 4); // resize to the new size
2754  std::vector<under_type>::iterator citer = cells_.begin();
2755  while (iter != end)
2756  {
2757  index_type *nodes = fill_ftor(*iter); // returns an array of length 4
2758  *citer = nodes[0];
2759  ++citer;
2760  *citer = nodes[1];
2761  ++citer;
2762  *citer = nodes[2];
2763  ++citer;
2764  *citer = nodes[3];
2765  ++citer; ++iter;
2766  }
2767  synchronize_lock_.unlock();
2768 }
2769 
2770 template <class Basis>
2771 PersistentTypeID
2772 TetVolMesh<Basis>::tetvolmesh_typeid(TetVolMesh<Basis>::type_name(-1), "Mesh",
2773  TetVolMesh<Basis>::maker);
2774 
2775 template <class Basis>
2776 const std::string
2777 TetVolMesh<Basis>::type_name(int n)
2778 {
2779  ASSERT((n >= -1) && n <= 1);
2780  if (n == -1)
2781  {
2782  static const std::string name = TypeNameGenerator::make_template_id(type_name(0), type_name(1));
2783  return name;
2784  }
2785  else if (n == 0)
2786  {
2787  static const std::string nm("TetVolMesh");
2788  return nm;
2789  }
2790  else
2791  {
2792  return find_type_name((Basis *)0);
2793  }
2794 }
2795 
2796 /* To generate a random point inside of a tetrahedron, we generate random
2797  barrycentric coordinates (independent random variables between 0 and
2798  1 that sum to 1) for the point. */
2799 template <class Basis>
2800 void
2801 TetVolMesh<Basis>::get_random_point(Core::Geometry::Point &p,
2802  typename Elem::index_type ei,
2803  FieldRNG &rng) const
2804 {
2805  // get positions of the vertices
2806  const Core::Geometry::Point &p0 = points_[cells_[ei*4+0]];
2807  const Core::Geometry::Point &p1 = points_[cells_[ei*4+1]];
2808  const Core::Geometry::Point &p2 = points_[cells_[ei*4+2]];
2809  const Core::Geometry::Point &p3 = points_[cells_[ei*4+3]];
2810 
2811  double t = rng();
2812  double u = rng();
2813  double v = rng();
2814 
2815  // Fold cube into prism.
2816  if (t + u > 1.0)
2817  {
2818  t = 1.0 - t;
2819  u = 1.0 - u;
2820  }
2821 
2822  // Fold prism into tet.
2823  if (u + v > 1.0)
2824  {
2825  const double tmp = v;
2826  v = 1.0 - t - u;
2827  u = 1.0 - tmp;
2828  }
2829  else if (t + u + v > 1.0)
2830  {
2831  const double tmp = v;
2832  v = t + u + v - 1.0;
2833  t = 1.0 - u - tmp;
2834  }
2835 
2836  // Convert to Barycentric and compute new point.
2837  const double a = 1.0 - t - u - v;
2838  p = Core::Geometry::Point(p0*a + p1*t + p2*u + p3*v);
2839 }
2840 
2841 template <class Basis>
2842 Core::Geometry::BBox
2843 TetVolMesh<Basis>::get_bounding_box() const
2844 {
2845  Core::Geometry::BBox result;
2846  typename Node::iterator ni, nie;
2847  begin(ni);
2848  end(nie);
2849  while (ni != nie)
2850  {
2851  result.extend(points_[*ni]);
2852  ++ni;
2853  }
2854  return (result);
2855 }
2856 
2857 template <class Basis>
2858 void
2859 TetVolMesh<Basis>::get_canonical_transform(Core::Geometry::Transform &t) const
2860 {
2861  t.load_identity();
2862  Core::Geometry::BBox bbox = get_bounding_box();
2863  t.pre_scale(bbox.diagonal());
2864  t.pre_translate(Core::Geometry::Vector(bbox.min()));
2865 }
2866 
2867 template <class Basis>
2868 void
2869 TetVolMesh<Basis>::transform(const Core::Geometry::Transform &t)
2870 {
2871  synchronize_lock_.lock();
2872 
2873  std::vector<Core::Geometry::Point>::iterator itr = points_.begin();
2874  std::vector<Core::Geometry::Point>::iterator eitr = points_.end();
2875  while (itr != eitr)
2876  {
2877  *itr = t.project(*itr);
2878  ++itr;
2879  }
2880 
2881  if (bbox_.valid())
2882  {
2883  bbox_.reset();
2884 
2885  // Compute bounding box
2886  typename Node::iterator ni, nie;
2887  begin(ni);
2888  end(nie);
2889  while (ni != nie)
2890  {
2891  bbox_.extend(point(*ni));
2892  ++ni;
2893  }
2894 
2895  // Compute epsilons associated with the bounding box
2896  epsilon_ = bbox_.diagonal().length()*1e-8;
2897  epsilon2_ = epsilon_*epsilon_;
2898  epsilon3_ = epsilon_*epsilon_*epsilon_;
2899 
2900  synchronized_ |= BOUNDING_BOX_E;
2901  }
2902 
2903  if (node_grid_) { node_grid_->transform(t); }
2904  if (elem_grid_) { elem_grid_->transform(t); }
2905 
2906  synchronize_lock_.unlock();
2907 }
2908 
2909 template <class Basis>
2910 void
2911 TetVolMesh<Basis>::remove_face(typename Node::index_type n1,
2912  typename Node::index_type n2,
2913  typename Node::index_type n3,
2914  typename Cell::index_type ci,
2915  bool /*table_only*/)
2916 {
2917  PFaceNode f(n1, n2, n3);
2918  typename face_nt::iterator iter = face_table_.find(f);
2919 
2920  if (iter == face_table_.end())
2921  {
2922  ASSERTFAIL("this face did not exist in the table");
2923  }
2924  PFaceNode found_face = (*iter).first;
2925  index_type found_idx = (*iter).second;
2926 
2927  index_type* cells = faces_[found_idx].cells_;
2928 
2929  if (cells[1] == MESH_NO_NEIGHBOR)
2930  {
2931  // this face belongs to only one cell
2932  //faces_.erase(faces_.begin() + found_idx);
2933  cells[0] = MESH_NO_NEIGHBOR;
2934  cells[1] = MESH_NO_NEIGHBOR;
2935  face_table_.erase(iter);
2936  }
2937  else
2938  {
2939  if (((faces_[found_idx].cells_[0])>>2) == ci)
2940  {
2941  cells[0] = cells[1];
2942  cells[1] = MESH_NO_NEIGHBOR;
2943  }
2944  else if (((faces_[found_idx].cells_[1])>>2) == ci)
2945  {
2946  cells[1] = MESH_NO_NEIGHBOR;
2947  }
2948  else
2949  {
2950  ASSERTFAIL("remove face: face does exist but is ");
2951  }
2952  }
2953 }
2954 
2955 template <class Basis>
2956 void
2957 TetVolMesh<Basis>::hash_face(typename Node::index_type n1,
2958  typename Node::index_type n2,
2959  typename Node::index_type n3,
2960  index_type combined_index,
2961  face_ht& table)
2962 {
2963  PFace f(n1, n2, n3);
2964 
2965  typename face_ht::iterator iter = table.find(f);
2966  if (iter == table.end())
2967  {
2968  f.cells_[0] = combined_index;
2969  table[f] = 0; // insert for the first time
2970  }
2971  else
2972  {
2973  PFace f = (*iter).first;
2974  if (f.cells_[1] != MESH_NO_NEIGHBOR)
2975  {
2976  std::cerr << "TetVolMesh - This Mesh has problems: Cells #"
2977  << (f.cells_[0]>>2) << ", #" << (f.cells_[1]>>2) << ", and #"
2978  << (combined_index>>2) << " are illegally adjacent." << std::endl;
2979  }
2980  else if ((f.cells_[0]>>2) == (combined_index>>2))
2981  {
2982  std::cerr << "TetVolMesh - This Mesh has problems: Cells #"
2983  << (f.cells_[0]>>2) << " and #" << (combined_index>>2)
2984  << " are the same." << std::endl;
2985  }
2986  else
2987  {
2988  f.cells_[1] = combined_index; // add this cell
2989  table.erase(iter);
2990  table[f] = 0;
2991  }
2992  }
2993 }
2994 
2995 template <class Basis>
2996 void
2997 TetVolMesh<Basis>::compute_faces()
2998 {
2999  typename Cell::iterator ci, cie;
3000  begin(ci); end(cie);
3001  typename Node::array_type arr(4);
3002 
3003  face_ht table;
3004 
3005  while (ci != cie)
3006  {
3007  get_nodes(arr, *ci);
3008  // 4 faces -- each is entered CCW from outside looking in
3009 
3010  index_type cell_index = (*ci)<<2;
3011  hash_face(arr[0], arr[2], arr[1], cell_index, table);
3012  hash_face(arr[1], arr[2], arr[3], cell_index + 1, table);
3013  hash_face(arr[0], arr[1], arr[3], cell_index + 2, table);
3014  hash_face(arr[0], arr[3], arr[2], cell_index + 3, table);
3015  ++ci;
3016  }
3017 
3018  faces_.resize(table.size());
3019 
3020  typename face_ht::iterator ht_iter = table.begin();
3021  typename face_ht::iterator ht_iter_end = table.end();
3022 
3023  boundary_faces_.resize(cells_.size() >> 2);
3024 
3025  index_type uidx = 0;
3026  while (ht_iter != ht_iter_end)
3027  {
3028  const PFace& pface = (*ht_iter).first;
3029  faces_[uidx].cells_[0] = pface.cells_[0];
3030  faces_[uidx].cells_[1] = pface.cells_[1];
3031  face_table_[PFaceNode(pface.nodes_[0],pface.nodes_[1],pface.nodes_[2])] = uidx;
3032 
3033  if (pface.cells_[1] == -1)
3034  {
3035  index_type cell = (pface.cells_[0]) >> 2;
3036  index_type face = (pface.cells_[0]) & 0x3;
3037  boundary_faces_[cell] |= 1 << face;
3038  }
3039  ++ht_iter; uidx++;
3040  }
3041 
3042  synchronize_lock_.lock();
3043  synchronized_ |= Mesh::FACES_E;
3044  synchronize_lock_.unlock();
3045 
3046 
3047 }
3048 
3049 
3050 template <class Basis>
3051 void
3052 TetVolMesh<Basis>::add_face(typename Node::index_type n1,
3053  typename Node::index_type n2,
3054  typename Node::index_type n3,
3055  index_type combined_index)
3056 {
3057 
3058  PFaceNode e(n1,n2,n3);
3059  typename face_nt::iterator nt_iter = face_table_.find(e);
3060 
3061  if (nt_iter == face_table_.end())
3062  {
3063  index_type uidx = static_cast<index_type>(faces_.size());
3064  face_table_[e] = uidx;
3065  PFaceCell c;
3066  faces_.push_back(c);
3067  faces_[uidx].cells_[0] = combined_index;
3068  }
3069  else
3070  {
3071  if (faces_[nt_iter->second].cells_[0] == MESH_NO_NEIGHBOR)
3072  {
3073  ASSERTFAIL("Face is in face_table_, but not in faces_ table");
3074  }
3075  if (faces_[nt_iter->second].cells_[1] != MESH_NO_NEIGHBOR)
3076  {
3077  ASSERTFAIL("Adding a face that is already connected twice");
3078  }
3079 
3080  faces_[nt_iter->second].cells_[1] = combined_index;
3081  }
3082 }
3083 
3084 template <class Basis>
3085 void
3086 TetVolMesh<Basis>::hash_edge(typename Node::index_type n1,
3087  typename Node::index_type n2,
3088  index_type combined_index,
3089  edge_ht& table)
3090 {
3091  if (n1 == n2) return;
3092  PEdge e(n1, n2);
3093  typename edge_ht::iterator iter = table.find(e);
3094 
3095  if (iter == table.end())
3096  {
3097  e.cells_.push_back(combined_index); // add this cell
3098  table[e] = 0; // insert for the first time
3099  }
3100  else
3101  {
3102  PEdge e = (*iter).first;
3103  e.cells_.push_back(combined_index); // add this cell
3104  table.erase(iter);
3105  table[e] = 0;
3106  }
3107 }
3108 
3109 template <class Basis>
3110 void
3111 TetVolMesh<Basis>::compute_edges()
3112 {
3113  typename Cell::iterator ci, cie;
3114  begin(ci); end(cie);
3115  edge_ht table;
3116 
3117  typename Node::array_type arr;
3118  while (ci != cie)
3119  {
3120  get_nodes(arr, *ci);
3121  index_type cell_index = (*ci)<<3;
3122  hash_edge(arr[0], arr[1], cell_index ,table);
3123  hash_edge(arr[1], arr[2], cell_index+1,table);
3124  hash_edge(arr[2], arr[0], cell_index+2,table);
3125  hash_edge(arr[3], arr[0], cell_index+3,table);
3126  hash_edge(arr[3], arr[1], cell_index+4,table);
3127  hash_edge(arr[3], arr[2], cell_index+5,table);
3128  ++ci;
3129  }
3130 
3131  edges_.resize(table.size());
3132 
3133  typename edge_ht::iterator ht_iter = table.begin();
3134  typename edge_ht::iterator ht_iter_end = table.end();
3135 
3136  index_type uidx = 0;
3137  while (ht_iter != ht_iter_end)
3138  {
3139  const PEdge& pedge = (*ht_iter).first;
3140  edges_[uidx].cells_ = pedge.cells_;
3141  edge_table_[PEdgeNode(pedge.nodes_[0],pedge.nodes_[1])] = uidx;
3142  ++ht_iter; uidx++;
3143  }
3144 
3145  synchronize_lock_.lock();
3146  synchronized_ |= Mesh::EDGES_E;
3147  synchronize_lock_.unlock();
3148 }
3149 
3150 template <class Basis>
3151 void
3152 TetVolMesh<Basis>::add_edge(typename Node::index_type n1,
3153  typename Node::index_type n2, index_type combined_index)
3154 {
3155  PEdgeNode e(n1,n2);
3156  typename edge_nt::iterator ht_iter = edge_table_.find(e);
3157  if (ht_iter == edge_table_.end())
3158  {
3159  index_type uidx = static_cast<index_type>(edges_.size());
3160  edge_table_[e] = uidx;
3161  PEdgeCell c;
3162  edges_.push_back(c);
3163  edges_[uidx].cells_.push_back(combined_index);
3164  }
3165  else
3166  {
3167  edges_[ht_iter->second].cells_.push_back(combined_index);
3168  }
3169 }
3170 
3171 template <class Basis>
3172 bool
3173 TetVolMesh<Basis>::synchronize(mask_type sync)
3174 {
3175  // Conversion table
3176  if (sync & (Mesh::ELEM_NEIGHBORS_E|Mesh::DELEMS_E))
3177  { sync |= Mesh::FACES_E; sync &= ~(Mesh::ELEM_NEIGHBORS_E|Mesh::DELEMS_E); }
3178 
3179  if (sync & Mesh::FIND_CLOSEST_NODE_E)
3180  { sync |= NODE_LOCATE_E; sync &= ~(Mesh::FIND_CLOSEST_NODE_E); }
3181 
3182  if (sync & Mesh::FIND_CLOSEST_ELEM_E)
3183  { sync |= ELEM_LOCATE_E|FACES_E; sync &= ~(Mesh::FIND_CLOSEST_ELEM_E); }
3184 
3185  if (sync & (Mesh::NODE_LOCATE_E|Mesh::ELEM_LOCATE_E)) sync |= Mesh::BOUNDING_BOX_E;
3186 
3187  // Filter out the only tables available
3188  sync &= (Mesh::EDGES_E|Mesh::FACES_E|
3189  Mesh::NODE_NEIGHBORS_E|Mesh::BOUNDING_BOX_E|
3190  Mesh::NODE_LOCATE_E|Mesh::ELEM_LOCATE_E);
3191 
3192  Core::Thread::UniqueLock lock(synchronize_lock_.get());
3193 
3194  // Only sync was hasn't been synched
3195  sync &= (~synchronized_);
3196 
3197  if (sync == Mesh::EDGES_E)
3198  {
3199  Synchronize Synchronize(this,sync);
3200  synchronize_lock_.unlock();
3201  Synchronize.run();
3202  synchronize_lock_.lock();
3203  }
3204  else if (sync & Mesh::EDGES_E)
3205  {
3206  mask_type tosync = Mesh::EDGES_E;
3207  Synchronize syncclass(this,tosync);
3208  boost::thread syncthread(syncclass);
3209  }
3210 
3211  if (sync == Mesh::FACES_E)
3212  {
3213  Synchronize Synchronize(this,sync);
3214  synchronize_lock_.unlock();
3215  Synchronize.run();
3216  synchronize_lock_.lock();
3217  }
3218  else if (sync & Mesh::FACES_E)
3219  {
3220  mask_type tosync = Mesh::FACES_E;
3221  Synchronize syncclass(this,tosync);
3222  boost::thread syncthread(syncclass);
3223  }
3224 
3225  if (sync == Mesh::NODE_NEIGHBORS_E)
3226  {
3227  Synchronize Synchronize(this,sync);
3228  synchronize_lock_.unlock();
3229  Synchronize.run();
3230  synchronize_lock_.lock();
3231  }
3232  else if (sync & Mesh::NODE_NEIGHBORS_E)
3233  {
3234  mask_type tosync = Mesh::NODE_NEIGHBORS_E;
3235  Synchronize syncclass(this,tosync);
3236  boost::thread syncthread(syncclass);
3237  }
3238 
3239  if (sync == Mesh::BOUNDING_BOX_E)
3240  {
3241  Synchronize Synchronize(this,sync);
3242  synchronize_lock_.unlock();
3243  Synchronize.run();
3244  synchronize_lock_.lock();
3245  }
3246  else if (sync & Mesh::BOUNDING_BOX_E)
3247  {
3248  mask_type tosync = Mesh::BOUNDING_BOX_E;
3249  Synchronize syncclass(this,tosync);
3250  boost::thread syncthread(syncclass);
3251  }
3252 
3253  if (sync == Mesh::NODE_LOCATE_E)
3254  {
3255  Synchronize Synchronize(this,sync);
3256  synchronize_lock_.unlock();
3257  Synchronize.run();
3258  synchronize_lock_.lock();
3259  }
3260  else if (sync & Mesh::NODE_LOCATE_E)
3261  {
3262  mask_type tosync = Mesh::NODE_LOCATE_E;
3263  Synchronize syncclass(this,tosync);
3264  boost::thread syncthread(syncclass);
3265  }
3266 
3267  if (sync == Mesh::ELEM_LOCATE_E)
3268  {
3269  Synchronize Synchronize(this,sync);
3270  synchronize_lock_.unlock();
3271  Synchronize.run();
3272  synchronize_lock_.lock();
3273  }
3274  else if (sync & Mesh::ELEM_LOCATE_E)
3275  {
3276  mask_type tosync = Mesh::ELEM_LOCATE_E;
3277  Synchronize syncclass(this,tosync);
3278  boost::thread syncthread(syncclass);
3279  }
3280 
3281  // Wait until threads are done
3282  while ((synchronized_ & sync) != sync)
3283  {
3284  synchronize_cond_.wait(lock);
3285  }
3286 
3287  return (true);
3288 }
3289 
3290 template <class Basis>
3291 bool
3292 TetVolMesh<Basis>::unsynchronize(mask_type /*sync*/)
3293 {
3294  return (true);
3295 }
3296 
3297 template <class Basis>
3298 bool
3299 TetVolMesh<Basis>::clear_synchronization()
3300 {
3301  // Undo marking the synchronization
3302  synchronize_lock_.lock();
3303 
3304  synchronized_ = Mesh::NODES_E | Mesh::ELEMS_E | Mesh::CELLS_E;
3305 
3306  // Free memory where possible
3307 
3308  faces_.clear();
3309  face_table_.clear();
3310  edges_.clear();
3311  edge_table_.clear();
3312  node_neighbors_.clear();
3313  boundary_faces_.clear();
3314 
3315  node_grid_.reset();
3316  elem_grid_.reset();
3317 
3318  synchronize_lock_.unlock();
3319 
3320  return (true);
3321 }
3322 
3323 template <class Basis>
3324 void
3325 TetVolMesh<Basis>::begin(typename TetVolMesh::Node::iterator &itr) const
3326 {
3327  ASSERTMSG(synchronized_ & Mesh::NODES_E,
3328  "Must call synchronize NODES_E on TetVolMesh first");
3329  itr = 0;
3330 }
3331 
3332 template <class Basis>
3333 void
3334 TetVolMesh<Basis>::end(typename TetVolMesh::Node::iterator &itr) const
3335 {
3336  ASSERTMSG(synchronized_ & Mesh::NODES_E,
3337  "Must call synchronize NODES_E on TetVolMesh first");
3338  itr = static_cast<typename Node::iterator>(points_.size());
3339 }
3340 
3341 template <class Basis>
3342 void
3343 TetVolMesh<Basis>::size(typename TetVolMesh::Node::size_type &s) const
3344 {
3345  ASSERTMSG(synchronized_ & Mesh::NODES_E,
3346  "Must call synchronize NODES_E on TetVolMesh first");
3347  s = points_.size();
3348 }
3349 
3350 template <class Basis>
3351 void
3352 TetVolMesh<Basis>::begin(typename TetVolMesh::Edge::iterator &itr) const
3353 {
3354  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
3355  "Must call synchronize EDGES_E on TetVolMesh first");
3356  itr = 0;
3357 }
3358 
3359 template <class Basis>
3360 void
3361 TetVolMesh<Basis>::end(typename TetVolMesh::Edge::iterator &itr) const
3362 {
3363  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
3364  "Must call synchronize EDGES_E on TetVolMesh first");
3365  itr = static_cast<index_type>(edges_.size());
3366 }
3367 
3368 template <class Basis>
3369 void
3370 TetVolMesh<Basis>::size(typename TetVolMesh::Edge::size_type &s) const
3371 {
3372  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
3373  "Must call synchronize EDGES_E on TetVolMesh first");
3374  s = static_cast<index_type>(edges_.size());
3375 }
3376 
3377 template <class Basis>
3378 void
3379 TetVolMesh<Basis>::begin(typename TetVolMesh::Face::iterator &itr) const
3380 {
3381  ASSERTMSG(synchronized_ & Mesh::FACES_E,
3382  "Must call synchronize FACES_E on TetVolMesh first");
3383  itr = 0;
3384 }
3385 
3386 template <class Basis>
3387 void
3388 TetVolMesh<Basis>::end(typename TetVolMesh::Face::iterator &itr) const
3389 {
3390  ASSERTMSG(synchronized_ & Mesh::FACES_E,
3391  "Must call synchronize FACES_E on TetVolMesh first");
3392  itr = static_cast<index_type>(faces_.size());
3393 }
3394 
3395 template <class Basis>
3396 void
3397 TetVolMesh<Basis>::size(typename TetVolMesh::Face::size_type &s) const
3398 {
3399  ASSERTMSG(synchronized_ & Mesh::FACES_E,
3400  "Must call synchronize FACES_E on TetVolMesh first");
3401  s = static_cast<index_type>(faces_.size());
3402 }
3403 
3404 template <class Basis>
3405 void
3406 TetVolMesh<Basis>::begin(typename TetVolMesh::Cell::iterator &itr) const
3407 {
3408  ASSERTMSG(synchronized_ & CELLS_E,
3409  "Must call synchronize CELLS_E on TetVolMesh first");
3410  itr = 0;
3411 }
3412 
3413 template <class Basis>
3414 void
3415 TetVolMesh<Basis>::end(typename TetVolMesh::Cell::iterator &itr) const
3416 {
3417  ASSERTMSG(synchronized_ & CELLS_E,
3418  "Must call synchronize CELLS_E on TetVolMesh first");
3419  itr = cells_.size() >> 2;
3420 }
3421 
3422 template <class Basis>
3423 void
3424 TetVolMesh<Basis>::size(typename TetVolMesh::Cell::size_type &s) const
3425 {
3426  ASSERTMSG(synchronized_ & CELLS_E,
3427  "Must call synchronize CELLS_E on TetVolMesh first");
3428  s = static_cast<size_type>(cells_.size() >> 2);
3429 }
3430 
3431 template <class Basis>
3432 void
3433 TetVolMesh<Basis>::create_cell_edges(typename Cell::index_type c)
3434 {
3435  typename Node::array_type arr;
3436  get_nodes(arr, c);
3437  index_type cell_index = (c<<3);
3438  add_edge(arr[0], arr[1], cell_index);
3439  add_edge(arr[1], arr[2], cell_index+1);
3440  add_edge(arr[2], arr[0], cell_index+2);
3441  add_edge(arr[3], arr[0], cell_index+3);
3442  add_edge(arr[3], arr[1], cell_index+4);
3443  add_edge(arr[3], arr[2], cell_index+5);
3444 }
3445 
3446 template <class Basis>
3447 void
3448 TetVolMesh<Basis>::remove_edge(typename Node::index_type n1,
3449  typename Node::index_type n2,
3450  typename Cell::index_type ci,
3451  bool table_only)
3452 {
3453  PEdgeNode e(n1, n2);
3454  typename edge_nt::iterator iter = edge_table_.find(e);
3455 
3456  if (iter == edge_table_.end())
3457  {
3458  ASSERTFAIL("this edge did not exist in the table");
3459  }
3460 
3461  PEdgeNode found_edge = (*iter).first;
3462  index_type found_idx = (*iter).second;
3463 
3464  const std::vector<index_type>& cells = edges_[found_idx].cells_;
3465  if (cells.size() < 2)
3466  {
3467  if ((cells[0] >>3) != ci)
3468  {
3469  ASSERTFAIL("this edge does exist in the table but is not connected to this cell");
3470  }
3471  edge_table_.erase(iter);
3472  if (!table_only) edges_[found_idx].cells_.clear();
3473  }
3474  else
3475  {
3476  typename std::vector<index_type>::iterator citer = edges_[found_idx].cells_.begin();
3477  typename std::vector<index_type>::iterator citer_end = edges_[found_idx].cells_.end();
3478 
3479  index_type cell_idx = ci;
3480  while(citer != citer_end)
3481  {
3482  if (((*citer)>>3)==cell_idx) {
3483  // temporary fix
3484  // should be replaced with algorithm
3485  edges_[found_idx].cells_.erase(citer);
3486  citer = edges_[found_idx].cells_.begin();
3487  citer_end = edges_[found_idx].cells_.end();
3488  }
3489  else {
3490  ++citer;
3491  }
3492  }
3493  }
3494 }
3495 
3496 template <class Basis>
3497 void
3498 TetVolMesh<Basis>::delete_cell_edges(typename Cell::index_type c,
3499  bool table_only)
3500 {
3501  typename Node::array_type arr;
3502  get_nodes(arr, c);
3503  remove_edge(arr[0], arr[1], c, table_only);
3504  remove_edge(arr[1], arr[2], c, table_only);
3505  remove_edge(arr[2], arr[0], c, table_only);
3506  remove_edge(arr[3], arr[0], c, table_only);
3507  remove_edge(arr[3], arr[1], c, table_only);
3508  remove_edge(arr[3], arr[2], c, table_only);
3509 }
3510 
3511 template <class Basis>
3512 void
3513 TetVolMesh<Basis>::create_cell_faces(typename Cell::index_type c)
3514 {
3515  typename Node::array_type arr;
3516  get_nodes(arr, c);
3517  index_type cell_index = c<<2;
3518  add_face(arr[0], arr[2], arr[1], cell_index);
3519  add_face(arr[1], arr[2], arr[3], cell_index+1);
3520  add_face(arr[0], arr[1], arr[3], cell_index+2);
3521  add_face(arr[0], arr[3], arr[2], cell_index+3);
3522 }
3523 
3524 
3525 template <class Basis>
3526 void
3527 TetVolMesh<Basis>::delete_cell_faces(typename Cell::index_type c,
3528  bool table_only)
3529 {
3530  typename Node::array_type arr;
3531  get_nodes(arr, c);
3532  remove_face(arr[0], arr[2], arr[1], c, table_only);
3533  remove_face(arr[1], arr[2], arr[3], c, table_only);
3534  remove_face(arr[0], arr[1], arr[3], c, table_only);
3535  remove_face(arr[0], arr[3], arr[2], c, table_only);
3536 }
3537 
3538 
3539 template <class Basis>
3540 void
3541 TetVolMesh<Basis>::create_cell_node_neighbors(typename Cell::index_type c)
3542 {
3543  for (index_type i = c*4; i < c*4+4; ++i)
3544  {
3545  node_neighbors_[cells_[i]].push_back(i);
3546  }
3547 }
3548 
3549 
3550 template <class Basis>
3551 void
3552 TetVolMesh<Basis>::delete_cell_node_neighbors(typename Cell::index_type c)
3553 {
3554  for (index_type i = c*4; i < c*4+4; ++i)
3555  {
3556  const index_type n = cells_[i];
3557  typename std::vector<typename Cell::index_type>::iterator node_cells_end =
3558  node_neighbors_[n].end();
3559  typename std::vector<typename Cell::index_type>::iterator cell =
3560  node_neighbors_[n].begin();
3561  while (cell != node_cells_end && (*cell) != i) ++cell;
3562 
3563  /// ASSERT that the node_neighbors_ structure contains this cell
3564  ASSERT(cell != node_cells_end);
3565 
3566  node_neighbors_[n].erase(cell);
3567  }
3568 }
3569 
3570 template <class Basis>
3571 void
3572 TetVolMesh<Basis>::create_cell_syncinfo(typename Cell::index_type ci)
3573 {
3574  synchronize_lock_.lock();
3575  if (synchronized_ & Mesh::NODE_NEIGHBORS_E)
3576  create_cell_node_neighbors(ci);
3577  if (synchronized_&Mesh::EDGES_E)
3578  create_cell_edges(ci);
3579  if (synchronized_&Mesh::FACES_E)
3580  create_cell_faces(ci);
3581  if (synchronized_ & Mesh::LOCATE_E)
3582  insert_elem_into_grid(ci);
3583  synchronize_lock_.unlock();
3584 }
3585 
3586 template <class Basis>
3587 void
3588 TetVolMesh<Basis>::delete_cell_syncinfo(typename Cell::index_type ci)
3589 {
3590  synchronize_lock_.lock();
3591  if (synchronized_ & Mesh::NODE_NEIGHBORS_E)
3592  delete_cell_node_neighbors(ci);
3593  if (synchronized_&Mesh::EDGES_E)
3594  delete_cell_edges(ci);
3595  if (synchronized_&Mesh::FACES_E)
3596  delete_cell_faces(ci);
3597  if (synchronized_ & Mesh::LOCATE_E)
3598  remove_elem_from_grid(ci);
3599  synchronize_lock_.unlock();
3600 }
3601 
3602 template <class Basis>
3603 void
3604 TetVolMesh<Basis>::create_cell_syncinfo_special(typename Cell::index_type ci)
3605 {
3606  synchronize_lock_.lock();
3607  if (synchronized_ & Mesh::NODE_NEIGHBORS_E)
3608  create_cell_node_neighbors(ci);
3609 
3610  if (synchronized_&Mesh::EDGES_E)
3611  {
3612  typename Node::array_type arr;
3613  get_nodes(arr, ci);
3614  index_type cell_index = ci<<3;
3615  add_edge(arr[0], arr[1], cell_index);
3616  add_edge(arr[1], arr[2], cell_index+1);
3617  add_edge(arr[2], arr[0], cell_index+2);
3618  add_edge(arr[3], arr[0], cell_index+3);
3619  add_edge(arr[3], arr[1], cell_index+4);
3620  add_edge(arr[3], arr[2], cell_index+5);
3621  }
3622 
3623  if (synchronized_&Mesh::FACES_E)
3624  {
3625  typename Node::array_type arr;
3626  get_nodes(arr, ci);
3627  index_type cell_index = ci<<2;
3628  add_face(arr[0], arr[2], arr[1], cell_index);
3629  add_face(arr[1], arr[2], arr[3], cell_index+1);
3630  add_face(arr[0], arr[1], arr[3], cell_index+2);
3631  add_face(arr[0], arr[3], arr[2], cell_index+3);
3632  }
3633  if (synchronized_ & Mesh::LOCATE_E)
3634  insert_elem_into_grid(ci);
3635  synchronize_lock_.unlock();
3636 }
3637 
3638 template <class Basis>
3639 void
3640 TetVolMesh<Basis>::delete_cell_syncinfo_special(typename Cell::index_type ci)
3641 {
3642  synchronize_lock_.lock();
3643  if (synchronized_ & Mesh::NODE_NEIGHBORS_E)
3644  delete_cell_node_neighbors(ci);
3645  if (synchronized_&Mesh::EDGES_E)
3646  delete_cell_edges(ci, true);
3647  if (synchronized_&Mesh::FACES_E)
3648  delete_cell_faces(ci, true);
3649  if (synchronized_ & Mesh::LOCATE_E)
3650  remove_elem_from_grid(ci);
3651  synchronize_lock_.unlock();
3652 }
3653 
3654 
3655 template <class Basis>
3656 void
3657 TetVolMesh<Basis>::set_nodes(typename Node::array_type &array,
3658  typename Cell::index_type idx)
3659 {
3660  ASSERT(array.size() == 4);
3661 
3662  delete_cell_syncinfo(idx);
3663 
3664  for (index_type n = 0; n < 4; ++n)
3665  cells_[idx * 4 + n] = array[n];
3666 
3667  create_cell_syncinfo(idx);
3668 }
3669 
3670 template <class Basis>
3671 void
3672 TetVolMesh<Basis>::compute_node_neighbors()
3673 {
3674  node_neighbors_.clear();
3675  node_neighbors_.resize(points_.size());
3676  index_type i, num_cells = cells_.size();
3677  for (i = 0; i < num_cells; i++)
3678  {
3679  node_neighbors_[cells_[i]].push_back(i);
3680  }
3681 
3682  synchronize_lock_.lock();
3683  synchronized_ |= Mesh::NODE_NEIGHBORS_E;
3684  synchronize_lock_.unlock();
3685 }
3686 
3687 template <class Basis>
3688 int
3689 TetVolMesh<Basis>::get_weights(const Core::Geometry::Point &p, typename Cell::array_type &l,
3690  double *w)
3691 {
3692  typename Cell::index_type idx;
3693  if (locate(idx, p))
3694  {
3695  l.resize(1);
3696  l[0] = idx;
3697  w[0] = 1.0;
3698  return 1;
3699  }
3700  return 0;
3701 }
3702 
3703 template <class Basis>
3704 int
3705 TetVolMesh<Basis>::get_weights(const Core::Geometry::Point &p, typename Node::array_type &l,
3706  double *w)
3707 {
3708  typename Cell::index_type idx;
3709  if (locate(idx, p))
3710  {
3711  get_nodes(l,idx);
3712  std::vector<double> coords(3);
3713  if (get_coords(coords, p, idx))
3714  {
3715  basis_.get_weights(coords, w);
3716  return basis_.dofs();
3717  }
3718  }
3719  return 0;
3720 }
3721 
3722 template <class Basis>
3723 void
3724 TetVolMesh<Basis>::insert_elem_into_grid(typename Cell::index_type ci)
3725 {
3726  /// @todo: This can crash if you insert a new cell outside of the grid.
3727  // Need to recompute grid at that point.
3728 
3729  const index_type idx = ci*4;
3730  Core::Geometry::BBox box;
3731  box.extend(points_[cells_[idx]]);
3732  box.extend(points_[cells_[idx+1]]);
3733  box.extend(points_[cells_[idx+2]]);
3734  box.extend(points_[cells_[idx+3]]);
3735  box.extend(epsilon_);
3736  elem_grid_->insert(ci, box);
3737 }
3738 
3739 
3740 template <class Basis>
3741 void
3742 TetVolMesh<Basis>::remove_elem_from_grid(typename Cell::index_type ci)
3743 {
3744  const index_type idx = ci*4;
3745  Core::Geometry::BBox box;
3746  box.extend(points_[cells_[idx]]);
3747  box.extend(points_[cells_[idx+1]]);
3748  box.extend(points_[cells_[idx+2]]);
3749  box.extend(points_[cells_[idx+3]]);
3750  box.extend(epsilon_);
3751  elem_grid_->remove(ci, box);
3752 }
3753 
3754 template <class Basis>
3755 void
3756 TetVolMesh<Basis>::insert_node_into_grid(typename Node::index_type ni)
3757 {
3758  /// @todo: This can crash if you insert a new cell outside of the grid.
3759  // Need to recompute grid at that point.
3760  node_grid_->insert(ni,points_[ni]);
3761 }
3762 
3763 template <class Basis>
3764 void
3765 TetVolMesh<Basis>::remove_node_from_grid(typename Node::index_type ni)
3766 {
3767  node_grid_->remove(ni,points_[ni]);
3768 }
3769 
3770 template <class Basis>
3771 void
3772 TetVolMesh<Basis>::compute_elem_grid()
3773 {
3774  if (bbox_.valid())
3775  {
3776  // Cubed root of number of cells to get a subdivision ballpark.
3777 
3778  typename Elem::size_type esz; size(esz);
3779 
3780  const size_type s =
3781  3*static_cast<size_type>((ceil(pow(static_cast<double>(esz) , (1.0/3.0))))/2.0 + 1.0);
3782 
3783  Core::Geometry::Vector diag = bbox_.diagonal();
3784  double trace = (diag.x()+diag.y()+diag.z());
3785  size_type sx = static_cast<size_type>(ceil(0.5+diag.x()/trace*s));
3786  size_type sy = static_cast<size_type>(ceil(0.5+diag.y()/trace*s));
3787  size_type sz = static_cast<size_type>(ceil(0.5+diag.z()/trace*s));
3788 
3789  Core::Geometry::BBox b = bbox_; b.extend(10*epsilon_);
3790  elem_grid_.reset(new SearchGridT<index_type>(sx, sy, sz, b.min(), b.max()));
3791 
3792  typename Elem::iterator ci, cie;
3793  begin(ci); end(cie);
3794  while(ci != cie)
3795  {
3796  insert_elem_into_grid(*ci);
3797  ++ci;
3798  }
3799  }
3800 
3801  synchronize_lock_.lock();
3802  synchronized_ |= Mesh::ELEM_LOCATE_E;
3803  synchronize_lock_.unlock();
3804 }
3805 
3806 template <class Basis>
3807 void
3808 TetVolMesh<Basis>::compute_node_grid()
3809 {
3810  ASSERTMSG(bbox_.valid(),"TetVolMesh BBox not valid");
3811  if (bbox_.valid())
3812  {
3813  // Cubed root of number of cells to get a subdivision ballpark.
3814 
3815  typename Elem::size_type esz; size(esz);
3816 
3817  const size_type s = 3*static_cast<size_type>
3818  ((ceil(pow(static_cast<double>(esz) , (1.0/3.0))))/2.0 + 1.0);
3819 
3820  Core::Geometry::Vector diag = bbox_.diagonal();
3821  double trace = (diag.x()+diag.y()+diag.z());
3822  size_type sx = static_cast<size_type>(ceil(0.5+diag.x()/trace*s));
3823  size_type sy = static_cast<size_type>(ceil(0.5+diag.y()/trace*s));
3824  size_type sz = static_cast<size_type>(ceil(0.5+diag.z()/trace*s));
3825 
3826  Core::Geometry::BBox b = bbox_; b.extend(10*epsilon_);
3827  node_grid_.reset(new SearchGridT<index_type>(sx, sy, sz, b.min(), b.max()));
3828 
3829  typename Node::iterator ni, nie;
3830  begin(ni); end(nie);
3831  while(ni != nie)
3832  {
3833  insert_node_into_grid(*ni);
3834  ++ni;
3835  }
3836  }
3837 
3838  synchronize_lock_.lock();
3839  synchronized_ |= Mesh::NODE_LOCATE_E;
3840  synchronize_lock_.unlock();
3841 }
3842 
3843 template <class Basis>
3844 void
3845 TetVolMesh<Basis>::compute_bounding_box()
3846 {
3847  bbox_.reset();
3848 
3849  // Compute bounding box
3850  typename Node::iterator ni, nie;
3851  begin(ni);
3852  end(nie);
3853 
3854  // make sure it does not fail on empty meshes
3855  if (ni == nie)
3856  {
3857  bbox_.extend(Core::Geometry::Point(0.0,0.0,0.0));
3858  bbox_.extend(Core::Geometry::Point(1.0,1.0,1.0));
3859  }
3860 
3861  while (ni != nie)
3862  {
3863  bbox_.extend(point(*ni));
3864  ++ni;
3865  }
3866 
3867  // Compute epsilons associated with the bounding box
3868  epsilon_ = bbox_.diagonal().length()*1e-8;
3869  epsilon2_ = epsilon_*epsilon_;
3870  epsilon3_ = epsilon_*epsilon_*epsilon_;
3871 
3872  synchronize_lock_.lock();
3873  synchronized_ |= BOUNDING_BOX_E;
3874  synchronize_lock_.unlock();
3875 }
3876 
3877 template <class Basis>
3878 typename TetVolMesh<Basis>::Node::index_type
3879 TetVolMesh<Basis>::add_find_point(const Core::Geometry::Point &p, double /*err*/)
3880 {
3881  typename Node::index_type i;
3882  if (locate(i, p) && (points_[i] - p).length2() < epsilon2_)
3883  {
3884  return i;
3885  }
3886  else
3887  {
3888  points_.push_back(p);
3889  if (synchronized_ & Mesh::NODE_NEIGHBORS_E)
3890  {
3891  synchronize_lock_.lock();
3892  node_neighbors_.push_back(std::vector<typename Cell::index_type>());
3893  synchronize_lock_.unlock();
3894  }
3895  return static_cast<typename Node::index_type>(points_.size() - 1);
3896  }
3897 }
3898 
3899 template <class Basis>
3900 typename TetVolMesh<Basis>::Elem::index_type
3901 TetVolMesh<Basis>::add_tet(typename Node::index_type a,
3902  typename Node::index_type b,
3903  typename Node::index_type c,
3904  typename Node::index_type d)
3905 {
3906  const index_type tet = static_cast<index_type>(cells_.size()) / 4;
3907  cells_.push_back(a);
3908  cells_.push_back(b);
3909  cells_.push_back(c);
3910  cells_.push_back(d);
3911  return tet;
3912 }
3913 
3914 template <class Basis>
3915 typename TetVolMesh<Basis>::Node::index_type
3916 TetVolMesh<Basis>::add_point(const Core::Geometry::Point &p)
3917 {
3918  points_.push_back(p);
3919  return static_cast<typename Node::index_type>(points_.size() - 1);
3920 }
3921 
3922 
3923 template <class Basis>
3924 typename TetVolMesh<Basis>::Elem::index_type
3925 TetVolMesh<Basis>::add_tet(const Core::Geometry::Point &p0, const Core::Geometry::Point &p1,
3926  const Core::Geometry::Point &p2, const Core::Geometry::Point &p3)
3927 {
3928  return add_tet(add_find_point(p0), add_find_point(p1),
3929  add_find_point(p2), add_find_point(p3));
3930 }
3931 
3932 template <class Basis>
3933 typename TetVolMesh<Basis>::Elem::index_type
3934 TetVolMesh<Basis>::add_tet_pos(typename Node::index_type a,
3935  typename Node::index_type b,
3936  typename Node::index_type c,
3937  typename Node::index_type d)
3938 {
3939  const index_type tet = static_cast<index_type>(cells_.size()) / 4;
3940  const Core::Geometry::Point &p0 = point(a);
3941  const Core::Geometry::Point &p1 = point(b);
3942  const Core::Geometry::Point &p2 = point(c);
3943  const Core::Geometry::Point &p3 = point(d);
3944 
3945  if (Dot(Cross(p1-p0,p2-p0),p3-p0) >= 0.0)
3946  {
3947  cells_.push_back(a);
3948  cells_.push_back(b);
3949  }
3950  else
3951  {
3952  cells_.push_back(b);
3953  cells_.push_back(a);
3954  }
3955  cells_.push_back(c);
3956  cells_.push_back(d);
3957  return tet;
3958 }
3959 
3960 template <class Basis>
3961 void
3962 TetVolMesh<Basis>::mod_tet_pos(typename TetVolMesh<Basis>::Elem::index_type ci,
3963  typename Node::index_type a,
3964  typename Node::index_type b,
3965  typename Node::index_type c,
3966  typename Node::index_type d)
3967 {
3968  const Core::Geometry::Point &p0 = point(a);
3969  const Core::Geometry::Point &p1 = point(b);
3970  const Core::Geometry::Point &p2 = point(c);
3971  const Core::Geometry::Point &p3 = point(d);
3972 
3973  if (Dot(Cross(p1-p0,p2-p0),p3-p0) >= 0.0)
3974  {
3975  cells_[ci*4+0] = a;
3976  cells_[ci*4+1] = b;
3977  }
3978  else
3979  {
3980  cells_[ci*4+0] = b;
3981  cells_[ci*4+1] = a;
3982  }
3983  cells_[ci*4+2] = c;
3984  cells_[ci*4+3] = d;
3985 }
3986 
3987 template <class Basis>
3988 void
3989 TetVolMesh<Basis>::delete_cells(std::set<index_type> &to_delete)
3990 {
3991  synchronize_lock_.lock();
3992  std::set<index_type>::reverse_iterator iter = to_delete.rbegin();
3993  while (iter != to_delete.rend())
3994  {
3995  // erase the correct cell
3996  typename TetVolMesh<Basis>::Cell::index_type ci = *iter++;
3997  index_type ind = ci * 4;
3998  std::vector<index_type>::iterator cb = cells_.begin() + ind;
3999  std::vector<index_type>::iterator ce = cb;
4000  ce+=4;
4001  cells_.erase(cb, ce);
4002  }
4003 
4004  synchronized_ &= ~Mesh::LOCATE_E;
4005  synchronized_ &= ~Mesh::NODE_NEIGHBORS_E;
4006  if (synchronized_ & Mesh::FACES_E)
4007  {
4008  synchronized_ &= ~Mesh::FACES_E;
4009  compute_faces();
4010  }
4011  if (synchronized_ & Mesh::EDGES_E)
4012  {
4013  synchronized_ &= ~Mesh::EDGES_E;
4014  compute_edges();
4015  }
4016  synchronize_lock_.unlock();
4017 }
4018 
4019 template <class Basis>
4020 void
4021 TetVolMesh<Basis>::delete_nodes(std::set<index_type> &to_delete)
4022 {
4023  synchronize_lock_.lock();
4024  std::set<index_type>::reverse_iterator iter = to_delete.rbegin();
4025  while (iter != to_delete.rend())
4026  {
4027  typename TetVolMesh::Node::index_type n = *iter++;
4028  std::vector<Core::Geometry::Point>::iterator pit = points_.begin() + n;
4029  points_.erase(pit);
4030  }
4031  synchronized_ &= ~Mesh::LOCATE_E;
4032  synchronized_ &= ~Mesh::NODE_NEIGHBORS_E;
4033  if (synchronized_ & Mesh::FACES_E)
4034  {
4035  synchronized_ &= ~Mesh::FACES_E;
4036  compute_faces();
4037  }
4038  if (synchronized_ & Mesh::EDGES_E)
4039  {
4040  synchronized_ &= ~Mesh::EDGES_E;
4041  compute_edges();
4042  }
4043  synchronize_lock_.unlock();
4044 }
4045 
4046 template <class Basis>
4047 bool
4048 TetVolMesh<Basis>::insert_node_in_cell(typename Cell::array_type &tets,
4049  typename Cell::index_type ci,
4050  typename Node::index_type &pi,
4051  const Core::Geometry::Point &p)
4052 {
4053  if (!inside(ci, p)) return false;
4054 
4055  pi = add_point(p);
4056 
4057  delete_cell_syncinfo_special(ci);
4058 
4059  tets.resize(4, ci);
4060  const index_type index = ci*4;
4061  tets[1] = add_tet_pos(cells_[index+0], cells_[index+3], cells_[index+1], pi);
4062  tets[2] = add_tet_pos(cells_[index+1], cells_[index+3], cells_[index+2], pi);
4063  tets[3] = add_tet_pos(cells_[index+0], cells_[index+2], cells_[index+3], pi);
4064 
4065  mod_tet_pos(ci, cells_[index+0], cells_[index+1], cells_[index+2], pi);
4066 
4067  create_cell_syncinfo_special(ci);
4068  create_cell_syncinfo_special(tets[1]);
4069  create_cell_syncinfo_special(tets[2]);
4070  create_cell_syncinfo_special(tets[3]);
4071 
4072  return true;
4073 }
4074 
4075 template <class Basis>
4076 void
4077 TetVolMesh<Basis>::insert_node_in_face(typename Cell::array_type &tets,
4078  typename Node::index_type pi,
4079  typename Cell::index_type ci,
4080  const PFaceNode &f)
4081 {
4082 
4083  // NOTE: This function assumes that the same operation has or will happen
4084  // to any neighbor cell across the face.
4085 
4086  index_type skip;
4087  for (skip = 0; skip < 4; skip++)
4088  {
4089  typename Node::index_type ni = cells_[ci*4+skip];
4090  if (ni != f.nodes_[0] && ni != f.nodes_[1] && ni != f.nodes_[2]) break;
4091  }
4092 
4093  delete_cell_syncinfo_special(ci);
4094 
4095  const index_type i = ci*4;
4096  tets.push_back(ci);
4097  if (skip == 0)
4098  {
4099  tets.push_back(add_tet_pos(cells_[i+0], cells_[i+3], cells_[i+1], pi));
4100  tets.push_back(add_tet_pos(cells_[i+0], cells_[i+2], cells_[i+3], pi));
4101  mod_tet_pos(ci, cells_[i+0], cells_[i+1], cells_[i+2], pi);
4102  }
4103  else if (skip == 1)
4104  {
4105  tets.push_back(add_tet_pos(cells_[i+0], cells_[i+3], cells_[i+1], pi));
4106  tets.push_back(add_tet_pos(cells_[i+1], cells_[i+3], cells_[i+2], pi));
4107  mod_tet_pos(ci, cells_[i+0], cells_[i+1], cells_[i+2], pi);
4108  }
4109  else if (skip == 2)
4110  {
4111  tets.push_back(add_tet_pos(cells_[i+1], cells_[i+3], cells_[i+2], pi));
4112  tets.push_back(add_tet_pos(cells_[i+0], cells_[i+2], cells_[i+3], pi));
4113  mod_tet_pos(ci, cells_[i+0], cells_[i+1], cells_[i+2], pi);
4114  }
4115  else if (skip == 3)
4116  {
4117  tets.push_back(add_tet_pos(cells_[i+0], cells_[i+3], cells_[i+1], pi));
4118  tets.push_back(add_tet_pos(cells_[i+1], cells_[i+3], cells_[i+2], pi));
4119  mod_tet_pos(ci, cells_[i+0], cells_[i+2], cells_[i+3], pi);
4120  }
4121 
4122  create_cell_syncinfo_special(ci);
4123  create_cell_syncinfo_special(tets[tets.size()-2]);
4124  create_cell_syncinfo_special(tets[tets.size()-1]);
4125 
4126 }
4127 
4128 template <class Basis>
4129 void
4130 TetVolMesh<Basis>::insert_node_in_edge(typename Cell::array_type &tets,
4131  typename Node::index_type pi,
4132  typename Cell::index_type ci,
4133  const PEdgeNode &e)
4134 {
4135  index_type skip1=-1, skip2=-1;
4136 
4137  for (index_type j = 0; j < 4; j++)
4138  {
4139  typename Node::index_type ni = cells_[ci*4+j];
4140  if (ni != e.nodes_[0] && ni != e.nodes_[1])
4141  {
4142  if (skip1 == -1) skip1 = j;
4143  else skip2 = j;
4144  }
4145  }
4146 
4147  delete_cell_syncinfo_special(ci);
4148 
4149  bool pushed = false;
4150  tets.push_back(ci);
4151  const index_type i = ci*4;
4152  if (skip1 != 0 && skip2 != 0)
4153  {
4154  // add 1 3 2 pi
4155  tets.push_back(add_tet_pos(cells_[i+1], cells_[i+3], cells_[i+2], pi));
4156  pushed = true;
4157  }
4158  if (skip1 != 1 && skip2 != 1)
4159  {
4160  // add 0 2 3 pi
4161  if (pushed)
4162  {
4163  mod_tet_pos(ci, cells_[i+0], cells_[i+2], cells_[i+3], pi);
4164  }
4165  else
4166  {
4167  tets.push_back(add_tet_pos(cells_[i+0], cells_[i+2], cells_[i+3], pi));
4168  }
4169  pushed = true;
4170  }
4171  if (skip1 != 2 && skip2 != 2)
4172  {
4173  // add 0 3 1 pi
4174  if (pushed)
4175  {
4176  mod_tet_pos(ci, cells_[i+0], cells_[i+3], cells_[i+1], pi);
4177  }
4178  else
4179  {
4180  tets.push_back(add_tet_pos(cells_[i+0], cells_[i+3], cells_[i+1], pi));
4181  }
4182  pushed = true;
4183  }
4184  if (skip1 != 3 && skip2 != 3)
4185  {
4186  // add 0 1 2 pi
4187  ASSERTMSG(pushed,
4188  "insert_node_in_cell_edge::skip1 or skip2 were invalid.");
4189  mod_tet_pos(ci, cells_[i+0], cells_[i+1], cells_[i+2], pi);
4190  }
4191 
4192  create_cell_syncinfo_special(ci);
4193  create_cell_syncinfo_special(tets[tets.size()-1]);
4194 }
4195 
4196 template <class Basis>
4197 bool
4198 TetVolMesh<Basis>::insert_node_in_elem(typename Elem::array_type &tets,
4199  typename Node::index_type &pi,
4200  typename Elem::index_type ci,
4201  const Core::Geometry::Point &p)
4202 {
4203 
4204  const Core::Geometry::Point &p0 = points_[cells_[ci*4 + 0]];
4205  const Core::Geometry::Point &p1 = points_[cells_[ci*4 + 1]];
4206  const Core::Geometry::Point &p2 = points_[cells_[ci*4 + 2]];
4207  const Core::Geometry::Point &p3 = points_[cells_[ci*4 + 3]];
4208 
4209  // Compute all the new tet areas.
4210  const double aerr = Abs(Dot(Cross(p1 - p0, p2 - p0), p3 - p0)) * 0.01;
4211  const double a0 = Abs(Dot(Cross(p1 - p, p2 - p), p3 - p));
4212  const double a1 = Abs(Dot(Cross(p - p0, p2 - p0), p3 - p0));
4213  const double a2 = Abs(Dot(Cross(p1 - p0, p - p0), p3 - p0));
4214  const double a3 = Abs(Dot(Cross(p1 - p0, p2 - p0), p - p0));
4215 
4216  mask_type mask = 0;
4217  if (a0 >= aerr && a0 >= epsilon3_) { mask |= 1; }
4218  if (a1 >= aerr && a1 >= epsilon3_) { mask |= 2; }
4219  if (a2 >= aerr && a2 >= epsilon3_) { mask |= 4; }
4220  if (a3 >= aerr && a3 >= epsilon3_) { mask |= 8; }
4221 
4222  // If we're completely inside then we do a normal 4 tet insert.
4223  // Test this first because it's most common.
4224  if (mask == 15) { return insert_node_in_cell(tets, ci, pi, p); }
4225 
4226  // If the tet is degenerate then we just return any corner and are done.
4227  else if (mask == 0)
4228  {
4229  tets.clear();
4230  tets.push_back(ci);
4231  pi = cells_[ci*4 + 0];
4232  return true;
4233  }
4234 
4235  // If we're on a corner, we're done. The corner is the point.
4236  else if (mask == 1)
4237  {
4238  tets.clear();
4239  tets.push_back(ci);
4240  pi = cells_[ci*4 + 0];
4241  return true;
4242  }
4243  else if (mask == 2)
4244  {
4245  tets.clear();
4246  tets.push_back(ci);
4247  pi = cells_[ci*4 + 1];
4248  return true;
4249  }
4250  else if (mask == 4)
4251  {
4252  tets.clear();
4253  tets.push_back(ci);
4254  pi = cells_[ci*4 + 2];
4255  return true;
4256  }
4257  else if (mask == 8)
4258  {
4259  tets.clear();
4260  tets.push_back(ci);
4261  pi = cells_[ci*4 + 3];
4262  return true;
4263  }
4264 
4265  // If we're on an edge, we do an edge insert.
4266  else if (mask == 3 || mask == 5 || mask == 6 ||
4267  mask == 9 || mask == 10 || mask == 12)
4268  {
4269  PEdgeNode etmp;
4270  if (mask == 3)
4271  {
4272  /* 0-1 edge */
4273  etmp = PEdgeNode(cells_[ci*4 + 0], cells_[ci*4 + 1]);
4274  }
4275  else if (mask == 5)
4276  {
4277  /* 0-2 edge */
4278  etmp = PEdgeNode(cells_[ci*4 + 0], cells_[ci*4 + 2]);
4279  }
4280  else if (mask == 6)
4281  {
4282  /* 1-2 edge */
4283  etmp = PEdgeNode(cells_[ci*4 + 1], cells_[ci*4 + 2]);
4284  }
4285  else if (mask == 9)
4286  {
4287  /* 0-3 edge */
4288  etmp = PEdgeNode(cells_[ci*4 + 0], cells_[ci*4 + 3]);
4289  }
4290  else if (mask == 10)
4291  {
4292  /* 1-3 edge */
4293  etmp = PEdgeNode(cells_[ci*4 + 1], cells_[ci*4 + 3]);
4294  }
4295  else if (mask == 12)
4296  {
4297  /* 2-3 edge */
4298  etmp = PEdgeNode(cells_[ci*4 + 2], cells_[ci*4 + 3]);
4299  }
4300 
4301  typename edge_nt::iterator iter = edge_table_.find(etmp);
4302  PEdgeNode e = iter->first;
4303  const std::vector<index_type>& cells = edges_[iter->second].cells_;
4304 
4305  pi = add_point(p);
4306  tets.clear();
4307  for (size_t i = 0; i < cells.size(); i++)
4308  {
4309  insert_node_in_edge(tets, pi, (cells[i]>>3), e);
4310  }
4311  }
4312 
4313  // If we're on a face, we do a face insert.
4314  else if (mask == 7 || mask == 11 || mask == 13 || mask == 14)
4315  {
4316  typename Face::index_type fi;
4317  PFaceNode ftmp;
4318  if (mask == 7)
4319  {
4320  /* on 0 1 2 face */
4321  ftmp = PFaceNode(cells_[ci*4 + 0], cells_[ci*4 + 1], cells_[ci*4 + 2]);
4322  }
4323  else if (mask == 11)
4324  {
4325  /* on 0 1 3 face */
4326  ftmp = PFaceNode(cells_[ci*4 + 0], cells_[ci*4 + 1], cells_[ci*4 + 3]);
4327  }
4328  else if (mask == 13)
4329  {
4330  /* on 0 2 3 face */
4331  ftmp = PFaceNode(cells_[ci*4 + 0], cells_[ci*4 + 2], cells_[ci*4 + 3]);
4332  }
4333  else if (mask == 14)
4334  {
4335  /* on 1 2 3 face */
4336  ftmp = PFaceNode(cells_[ci*4 + 1], cells_[ci*4 + 2], cells_[ci*4 + 3]);
4337  }
4338 
4339  typename face_nt::iterator iter = face_table_.find(ftmp);
4340  const PFaceNode& n = iter->first;
4341  const PFaceCell& f = faces_[iter->second];
4342  typename Cell::index_type nbr_tet =
4343  (ci == (f.cells_[0])>>2) ? ((f.cells_[1])>>2) : ((f.cells_[0])>>2);
4344 
4345  pi = add_point(p);
4346  tets.clear();
4347  insert_node_in_face(tets, pi, ci, n);
4348  if (nbr_tet != MESH_NO_NEIGHBOR)
4349  {
4350  insert_node_in_face(tets, pi, nbr_tet, n);
4351  }
4352  }
4353 
4354  return true;
4355 }
4356 
4357 template <class Basis>
4358 void
4359 TetVolMesh<Basis>::orient(typename Cell::index_type ci)
4360 {
4361  const Core::Geometry::Point &p0 = point(cells_[ci*4+0]);
4362  const Core::Geometry::Point &p1 = point(cells_[ci*4+1]);
4363  const Core::Geometry::Point &p2 = point(cells_[ci*4+2]);
4364  const Core::Geometry::Point &p3 = point(cells_[ci*4+3]);
4365 
4366  // Unsigned volumex6 of the tet.
4367  const double sgn = Dot(Cross(p1-p0,p2-p0),p3-p0);
4368  if (sgn < 0.0)
4369  {
4370  typename Node::index_type tmp = cells_[ci*4+0];
4371  cells_[ci*4+0] = cells_[ci*4+1];
4372  cells_[ci*4+1] = tmp;
4373  }
4374 }
4375 
4376 #define TETVOLMESH_VERSION 4
4377 
4378 template <class Basis>
4379 void
4380 TetVolMesh<Basis>::io(Piostream &stream)
4381 {
4382  const int version = stream.begin_class(type_name(-1),
4383  TETVOLMESH_VERSION);
4384  Mesh::io(stream);
4385 
4386  SCIRun::Pio(stream, points_);
4387  SCIRun::Pio_index(stream, cells_);
4388  if (version == 1)
4389  {
4390  std::vector<unsigned int> neighbors;
4391  SCIRun::Pio(stream, neighbors);
4392  }
4393 
4394  if (version >= 3)
4395  {
4396  basis_.io(stream);
4397  }
4398 
4399  stream.end_class();
4400 
4401  if (stream.reading())
4402  {
4403  synchronized_ = NODES_E | CELLS_E;
4404 
4405  vmesh_.reset(CreateVTetVolMesh(this));
4406  }
4407 }
4408 
4409 template <class Basis>
4410 const TypeDescription* get_type_description(TetVolMesh<Basis> *)
4411 {
4412  static TypeDescription *td = 0;
4413  if (!td)
4414  {
4415  const TypeDescription *sub = get_type_description((Basis*)0);
4416  TypeDescription::td_vec *subs = new TypeDescription::td_vec(1);
4417  (*subs)[0] = sub;
4418  td = new TypeDescription("TetVolMesh", subs,
4419  std::string(__FILE__),
4420  "SCIRun",
4421  TypeDescription::MESH_E);
4422  }
4423  return td;
4424 }
4425 
4426 template <class Basis>
4427 const TypeDescription*
4428 TetVolMesh<Basis>::get_type_description() const
4429 {
4430  return SCIRun::get_type_description((TetVolMesh<Basis> *)0);
4431 }
4432 
4433 template <class Basis>
4434 const TypeDescription*
4435 TetVolMesh<Basis>::node_type_description()
4436 {
4437  static TypeDescription *td = 0;
4438  if (!td)
4439  {
4440  const TypeDescription *me =
4441  SCIRun::get_type_description((TetVolMesh<Basis> *)0);
4442  td = new TypeDescription(me->get_name() + "::Node",
4443  std::string(__FILE__),
4444  "SCIRun",
4445  TypeDescription::MESH_E);
4446  }
4447  return td;
4448 }
4449 
4450 template <class Basis>
4451 const TypeDescription*
4452 TetVolMesh<Basis>::edge_type_description()
4453 {
4454  static TypeDescription *td = 0;
4455  if (!td)
4456  {
4457  const TypeDescription *me =
4458  SCIRun::get_type_description((TetVolMesh<Basis> *)0);
4459  td = new TypeDescription(me->get_name() + "::Edge",
4460  std::string(__FILE__),
4461  "SCIRun",
4462  TypeDescription::MESH_E);
4463  }
4464  return td;
4465 }
4466 
4467 template <class Basis>
4468 const TypeDescription*
4469 TetVolMesh<Basis>::face_type_description()
4470 {
4471  static TypeDescription *td = 0;
4472  if (!td)
4473  {
4474  const TypeDescription *me =
4475  SCIRun::get_type_description((TetVolMesh<Basis> *)0);
4476  td = new TypeDescription(me->get_name() + "::Face",
4477  std::string(__FILE__),
4478  "SCIRun",
4479  TypeDescription::MESH_E);
4480  }
4481  return td;
4482 }
4483 
4484 template <class Basis>
4485 const TypeDescription*
4486 TetVolMesh<Basis>::cell_type_description()
4487 {
4488  static TypeDescription *td = 0;
4489  if (!td)
4490  {
4491  const TypeDescription *me =
4492  SCIRun::get_type_description((TetVolMesh<Basis> *)0);
4493  td = new TypeDescription(me->get_name() + "::Cell",
4494  std::string(__FILE__),
4495  "SCIRun",
4496  TypeDescription::MESH_E);
4497  }
4498  return td;
4499 }
4500 
4501 } // namespace SCIRun
4502 
4503 #endif
SCIRun::TetVolMesh::~TetVolMesh
virtual ~TetVolMesh()
Destructor.
Definition: TetVolMesh.h:2724
SCIRun::TetVolMesh::ElemData::edge3_index
index_type edge3_index() const
Definition: TetVolMesh.h:249
SCIRun::TetVolMesh::get_center
void get_center(Core::Geometry::Point &result, typename Cell::index_type idx) const
Definition: TetVolMesh.h:558
SCIRun::TetVolMesh::get_nodes
void get_nodes(typename Node::array_type &array, typename Cell::index_type idx) const
Definition: TetVolMesh.h:459
StackVector.h
Interface to statically allocated std::vector class.
SCIRun::VUnstructuredMesh
Definition: VUnstructuredMesh.h:41
SCIRun::TetVolMesh::insert_node_in_face
void insert_node_in_face(typename Cell::array_type &tets, typename Node::index_type ni, typename Cell::index_type ci, const PFaceNode &pf)
Definition: TetVolMesh.h:4077
SCIRun::TetVolMesh::PEdge::PEdge
PEdge(typename Node::index_type n1, typename Node::index_type n2)
Definition: TetVolMesh.h:2464
SCIRun::TetVolMesh::node_neighbors_
std::vector< std::vector< typename Cell::index_type > > node_neighbors_
Definition: TetVolMesh.h:2601
MeshSupport.h
SCIRun::TetVolMesh::get_faces
void get_faces(typename Face::array_type &array, typename Cell::index_type idx) const
Definition: TetVolMesh.h:486
SCIRun::TetVolMesh::add_tet
Elem::index_type add_tet(typename Node::index_type a, typename Node::index_type b, typename Node::index_type c, typename Node::index_type d)
Definition: TetVolMesh.h:3901
SCIRun::TetVolMesh::ElemData::edge1_index
index_type edge1_index() const
Definition: TetVolMesh.h:239
TetCubicHmt.h
SCIRun::Piostream::reading
bool reading() const
Definition: Persistent.h:164
SCIRun::NodeIterator
Distinct type for node FieldIterator.
Definition: FieldIterator.h:89
SCIRun::TetVolMesh::get_area
double get_area(typename Face::index_type idx) const
Definition: TetVolMesh.h:595
SCIRun::TetVolMesh::compute_bounding_box
void compute_bounding_box()
Definition: TetVolMesh.h:3845
SCIRun::TetVolMesh::get_neighbors
void get_neighbors(typename Elem::array_type &array, typename Elem::index_type elem) const
Definition: TetVolMesh.h:619
SCIRun::TetVolMesh::compute_node_neighbors
void compute_node_neighbors()
Definition: TetVolMesh.h:3672
SCIRun::TetVolMesh::get_cells
void get_cells(typename Cell::array_type &array, typename Edge::index_type idx) const
Definition: TetVolMesh.h:493
Point.h
SCIRun::FieldRNG
Definition: FieldRNG.h:37
SCIRun::TetVolMesh::get_cells
void get_cells(typename Cell::array_type &array, typename Node::index_type idx) const
Definition: TetVolMesh.h:490
SCIRun::TetVolMesh::PEdgeNode::PEdgeNode
PEdgeNode(typename Node::index_type n1, typename Node::index_type n2)
Definition: TetVolMesh.h:2420
SCIRun::TetVolMesh::get_faces
void get_faces(typename Face::array_type &array, typename Edge::index_type idx) const
Definition: TetVolMesh.h:479
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::TetVolMesh::cell_type_description
static const TypeDescription * cell_type_description()
Definition: TetVolMesh.h:4486
SCIRun::TetVolMesh::EdgeHash::operator()
size_t operator()(const PEDGE &e) const
This is the hash function.
Definition: TetVolMesh.h:2521
SCIRun::FaceIterator
Distinct type for face Iterator.
Definition: FieldIterator.h:121
SCIRun::TetVolMesh::insert_node_in_cell
bool insert_node_in_cell(typename Cell::array_type &tets, typename Cell::index_type ci, typename Node::index_type &ni, const Core::Geometry::Point &p)
Always creates 4 tets, does not handle element boundaries properly.
Definition: TetVolMesh.h:4048
SCIRun::Core::Thread::ConditionVariable
Definition: ConditionVariable.h:45
SCIRun::Mesh::LOCATE_E
Definition: Mesh.h:85
Vector.h
SCIRun::Core::Geometry::Transform::project
Vector project(const Vector &p) const
Definition: Transform.cc:425
SCIRun::TetVolMesh::ElemData::edge0_index
index_type edge0_index() const
the following designed to coordinate with ::get_edges
Definition: TetVolMesh.h:234
SCIRun::TetVolMesh::inverse_jacobian
double inverse_jacobian(const VECTOR &coords, INDEX idx, double *Ji) const
Definition: TetVolMesh.h:772
SCIRun::Mesh::IRREGULAR
Definition: Mesh.h:64
SCIRun::TetVolMesh::ElemData
Definition: TetVolMesh.h:194
SCIRun::TetVolMesh::get_nodes
void get_nodes(typename Node::array_type &array, typename Face::index_type idx) const
Definition: TetVolMesh.h:456
SCIRun::TetVolMesh::ElemData::node0_index
index_type node0_index() const
the following designed to coordinate with ::get_nodes
Definition: TetVolMesh.h:212
SCIRun::Core::Geometry::BBox::max
Point max() const
Definition: BBox.h:195
SCIRun::TetVolMesh::get_point
void get_point(Core::Geometry::Point &result, typename Node::index_type index) const
Access the nodes of the mesh.
Definition: TetVolMesh.h:648
SCIRun::TetVolMesh::index_type
SCIRun::index_type index_type
Definition: TetVolMesh.h:146
SCIRun::TetVolMesh::Synchronize
Definition: TetVolMesh.h:296
ASSERTMSG
#define ASSERTMSG(condition, message)
Definition: Exception.h:113
SCIRun::TetVolMesh::vmesh
virtual VMesh * vmesh()
Access point to virtual interface.
Definition: TetVolMesh.h:368
SCIRun::TetVolMesh::epsilon3_
double epsilon3_
Definition: TetVolMesh.h:2625
SCIRun::TetVolMesh::get_elem_neighbor
bool get_elem_neighbor(INDEX1 &neighbor, INDEX1 elem, INDEX2 delem) const
Definition: TetVolMesh.h:1850
SCIRun::Abs
double Abs(double d)
Definition: MiscMath.h:369
SCIRun::TetVolMesh::node_reserve
void node_reserve(size_type s)
Definition: TetVolMesh.h:703
SCIRun::TetVolMesh::basis_order
virtual int basis_order()
Definition: TetVolMesh.h:377
SCIRun::TetVolMesh::get_face_center
void get_face_center(Core::Geometry::Point &p, INDEX idx) const
Definition: TetVolMesh.h:2212
SCIRun::TetVolMesh::Node::index_type
NodeIndex< under_type > index_type
Definition: TetVolMesh.h:155
SCIRun::TetVolMesh::get_edges_from_node
void get_edges_from_node(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1698
SCIRun::TetVolMesh::Face::size_type
FaceIndex< under_type > size_type
Definition: TetVolMesh.h:171
SCIRun::TetVolMesh::Node::size_type
NodeIndex< under_type > size_type
Definition: TetVolMesh.h:157
SCIRun::TetVolMesh::elem_type_description
static const TypeDescription * elem_type_description()
Definition: TetVolMesh.h:1398
SCIRun::TetVolMesh::point
const Core::Geometry::Point & point(typename Node::index_type i)
Definition: TetVolMesh.h:2255
SCIRun::Piostream
Definition: Persistent.h:89
SCIRun::TetVolMesh::PEdge
Edge information.
Definition: TetVolMesh.h:2461
SCIRun::TetVolMesh::locate
bool locate(typename Edge::index_type &edge, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:626
SCIRun::Core::Geometry::Transform::pre_scale
void pre_scale(const Vector &)
Definition: Transform.cc:193
SCIRun::TetVolMesh::locate_face
bool locate_face(INDEX &face, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:2055
TETVOLMESH_VERSION
#define TETVOLMESH_VERSION
Definition: TetVolMesh.h:4376
SCIRun::TetVolMesh::EdgeHash
Definition: TetVolMesh.h:2505
SCIRun::TetVolMesh::ElemData::node2
const Core::Geometry::Point & node2() const
Definition: TetVolMesh.h:275
SCIRun::TetVolMesh::get_node_center
void get_node_center(Core::Geometry::Point &p, INDEX idx) const
Definition: TetVolMesh.h:2194
SCIRun::TetVolMesh::clear_synchronization
bool clear_synchronization()
Definition: TetVolMesh.h:3299
SCIRun::TetVolMesh::PEdgeNode
Edge information.
Definition: TetVolMesh.h:2410
SCIRun::TetVolMesh::PFaceNode::PFaceNode
PFaceNode(typename Node::index_type n1, typename Node::index_type n2, typename Node::index_type n3)
Definition: TetVolMesh.h:2343
SCIRun::TetVolMesh::insert_elem_into_grid
void insert_elem_into_grid(typename Elem::index_type ci)
Definition: TetVolMesh.h:3724
SCIRun::TetVolMesh::EdgeHash::low_mask
static const int low_mask
Definition: TetVolMesh.h:2517
SCIRun::TetVolMesh::ElemData::node3_index
index_type node3_index() const
Definition: TetVolMesh.h:227
TetLinearLgn.h
SCIRun::TetVolMesh::get_weights
int get_weights(const Core::Geometry::Point &, typename Edge::array_type &, double *)
Definition: TetVolMesh.h:641
SCIRun::TetVolMesh::FaceHash::sz_third_int
static const int sz_third_int
These are for our own use (making the hash function).
Definition: TetVolMesh.h:2481
SCIRun::TetVolMesh::ElemData::node2_index
index_type node2_index() const
Definition: TetVolMesh.h:222
SCIRun::compute_checksum
int compute_checksum(T *data, std::size_t length)
Definition: CheckSum.h:38
SCIRun::TetVolMesh::mod_tet_pos
void mod_tet_pos(typename Elem::index_type tet, typename Node::index_type a, typename Node::index_type b, typename Node::index_type c, typename Node::index_type d)
Definition: TetVolMesh.h:3962
SCIRun::TetVolMesh::get_basis
Basis & get_basis()
Get the basis class.
Definition: TetVolMesh.h:416
SCIRun::TetVolMesh::get_normal
void get_normal(Core::Geometry::Vector &result, VECTOR &coords, INDEX1 eidx, INDEX2 fidx)
Get the normals at the outside of the element.
Definition: TetVolMesh.h:660
SCIRun::TetVolMesh::EdgeHash::bucket_size
static const size_t bucket_size
These are needed by the hash_map particularly.
Definition: TetVolMesh.h:2510
SCIRun::InverseMatrix3P
double InverseMatrix3P(const PointVector &p, double *q)
Inline templated inverse matrix.
Definition: Locate.h:72
SCIRun::TetVolMesh::remove_edge
void remove_edge(typename Node::index_type n1, typename Node::index_type n2, typename Cell::index_type ci, bool table_only=false)
Definition: TetVolMesh.h:3448
SCIRun::Core::Geometry::Point::y
void y(const double)
Definition: Point.h:135
SCIRun::TetVolMesh::resize_nodes
void resize_nodes(size_type s)
Definition: TetVolMesh.h:705
SCIRun::TetVolMesh::create_cell_edges
void create_cell_edges(typename Cell::index_type)
Definition: TetVolMesh.h:3433
SCIRun::Persistent
Definition: Persistent.h:187
SCIRun::TetVolMesh::Elem
Cell Elem
Definition: TetVolMesh.h:185
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::TetVolMesh::get_points
std::vector< Core::Geometry::Point > & get_points()
must detach, if altering points!
Definition: TetVolMesh.h:1479
SCIRun::TetVolMesh::delete_cell_node_neighbors
void delete_cell_node_neighbors(typename Cell::index_type)
Definition: TetVolMesh.h:3552
SCIRun::TetVolMesh::PFaceNode::PFaceNode
PFaceNode()
3 nodes makes a face.
Definition: TetVolMesh.h:2336
SCIRun::TypeDescription::MESH_E
Definition: TypeDescription.h:50
SCIRun::TetVolMesh::mask_type
SCIRun::mask_type mask_type
Definition: TetVolMesh.h:148
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::TetVolMesh::face_ct
std::vector< PFaceCell > face_ct
Definition: TetVolMesh.h:2559
SCIRun::Mesh::ELEM_LOCATE_E
Definition: Mesh.h:84
SCIRun::Core::Geometry::Point
Definition: Point.h:49
SCIRun::TetVolMesh::PEdgeNode::operator==
bool operator==(const PEdgeNode &e) const
true if both have the same nodes (order does not matter)
Definition: TetVolMesh.h:2436
SCIRun::TypeDescription
Definition: TypeDescription.h:45
SCIRun::TetVolMesh::get_edges
void get_edges(typename Edge::array_type &array, typename Cell::index_type idx) const
Definition: TetVolMesh.h:472
SCIRun::TetVolMesh::inside
bool inside(INDEX idx, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:2258
SCISHARE
#define SCISHARE
Definition: share.h:39
SCIRun::TypeDescription::td_vec
std::vector< const TypeDescription * > td_vec
Definition: TypeDescription.h:56
SCIRun::TetVolMesh::Cell::index_type
CellIndex< under_type > index_type
Definition: TetVolMesh.h:176
SCIRun::TetVolMesh::FaceHash::up_mask
static const int up_mask
Definition: TetVolMesh.h:2482
SCIRun::TetVolMesh::set_point
void set_point(const Core::Geometry::Point &point, typename Node::index_type index)
Definition: TetVolMesh.h:650
SCIRun::TetVolMesh::face_type_description
static const TypeDescription * face_type_description()
Definition: TetVolMesh.h:4469
SCIRun::TetVolMesh::ElemData::index_type
TetVolMesh< Basis >::index_type index_type
Definition: TetVolMesh.h:197
SCIRun::Mesh
Definition: Mesh.h:45
SCIRun::TetVolMesh::find_closest_nodes
bool find_closest_nodes(ARRAY1 &distances, ARRAY2 &nodes, const Core::Geometry::Point &p, double maxdist) const
Definition: TetVolMesh.h:1065
SCIRun::TetVolMesh::create_cell_syncinfo
void create_cell_syncinfo(typename Cell::index_type ci)
Definition: TetVolMesh.h:3572
SCIRun::TetVolMesh::Edge::size_type
EdgeIndex< under_type > size_type
Definition: TetVolMesh.h:164
SCIRun::Core::Geometry::BBox
Definition: BBox.h:46
SCIRun::MeshHandle
boost::shared_ptr< Mesh > MeshHandle
Definition: DatatypeFwd.h:67
SCIRun::Core::Thread::Mutex
Definition: Mutex.h:43
SCIRun::TetVolMesh::pwl_approx_edge
void pwl_approx_edge(std::vector< VECTOR > &coords, INDEX ci, unsigned which_edge, unsigned div_per_unit) const
Definition: TetVolMesh.h:532
SCIRun::TetVolMesh::PEdgeNode::PEdgeNode
PEdgeNode()
Definition: TetVolMesh.h:2414
SCIRun::TetVolMesh::Node::iterator
NodeIterator< under_type > iterator
Definition: TetVolMesh.h:156
MESH_NO_NEIGHBOR
#define MESH_NO_NEIGHBOR
Definition: MeshTraits.h:67
SCIRun::TetVolMesh::insert_node_into_grid
void insert_node_into_grid(typename Node::index_type ci)
Definition: TetVolMesh.h:3756
SCIRun::CellIterator
Distinct type for cell Iterator.
Definition: FieldIterator.h:134
SCIRun::TetVolMesh::delete_cell_syncinfo
void delete_cell_syncinfo(typename Cell::index_type ci)
Definition: TetVolMesh.h:3588
SCIRun::TetVolMesh::PEdgeCell::shared
bool shared() const
Definition: TetVolMesh.h:2457
VMesh.h
SCIRun::mask_type
unsigned int mask_type
Definition: Types.h:45
SCIRun::TetVolMesh::topology_geometry
virtual int topology_geometry() const
Definition: TetVolMesh.h:385
SCIRun::TetVolMesh::get_volume
double get_volume(typename Cell::index_type idx) const
Definition: TetVolMesh.h:597
SCIRun::TetVolMesh::FaceHash::operator()
bool operator()(const PFACE &f1, const PFACE &f2) const
This should return less than rather than equal to.
Definition: TetVolMesh.h:2497
SCIRun::TetVolMesh::PEdgeNode::operator<
bool operator<(const PEdgeNode &e) const
Definition: TetVolMesh.h:2443
SCIRun::TetVolMesh::add_tet_pos
Elem::index_type add_tet_pos(typename Node::index_type a, typename Node::index_type b, typename Node::index_type c, typename Node::index_type d)
Definition: TetVolMesh.h:3934
SCIRun::TetVolMesh::PFaceNode::operator<
bool operator<(const PFaceNode &f) const
Definition: TetVolMesh.h:2389
SCIRun::TetVolMesh::get_nodes
void get_nodes(typename Node::array_type &array, typename Edge::index_type idx) const
Definition: TetVolMesh.h:453
FieldRNG.h
SCIRun::TetVolMesh::locate_edge
bool locate_edge(INDEX &edge, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:2025
SCIRun::TetVolMesh::FaceHash::low_mask
static const int low_mask
Definition: TetVolMesh.h:2484
SCIRun::TetVolMesh::jacobian
void jacobian(const VECTOR &coords, INDEX idx, double *J) const
Definition: TetVolMesh.h:752
SCIRun::TetVolMesh::get_nodes_from_edge
void get_nodes_from_edge(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1497
SCIRun::CreateVTetVolMesh
VMesh * CreateVTetVolMesh(MESH *)
Definition: TetVolMesh.h:83
ASSERT
#define ASSERT(condition)
Definition: Assert.h:110
SCIRun::TetVolMesh::TetVolMesh
TetVolMesh()
Construct a new mesh.
Definition: TetVolMesh.h:2661
SCIRun::TetVolMesh::get_size
double get_size(typename Node::index_type) const
Get the size of an element (length, area, volume)
Definition: TetVolMesh.h:562
SCIRun::TetVolMesh::begin
void begin(typename Node::iterator &) const
begin/end iterators
Definition: TetVolMesh.h:3325
SCIRun::TetVolMesh::set_nodes
void set_nodes(typename Node::array_type &, typename Cell::index_type)
Definition: TetVolMesh.h:3657
SCIRun::TetVolMesh::get_faces
void get_faces(typename Face::array_type &array, typename Node::index_type idx) const
Definition: TetVolMesh.h:476
SCIRun::TetVolMesh::locate_elems
bool locate_elems(ARRAY &array, const Core::Geometry::BBox &b) const
Definition: TetVolMesh.h:2119
SCIRun::TetVolMesh::create_cell_syncinfo_special
void create_cell_syncinfo_special(typename Cell::index_type ci)
Definition: TetVolMesh.h:3604
SCIRun::TetVolMesh::get_coords
bool get_coords(VECTOR &coords, const Core::Geometry::Point &p, INDEX idx) const
Definition: TetVolMesh.h:712
SCIRun::TetVolMesh::compute_checksum
int compute_checksum()
Definition: TetVolMesh.h:2652
Plane.h
SCIRun::TetVolMesh::handle_type
boost::shared_ptr< TetVolMesh< Basis > > handle_type
Definition: TetVolMesh.h:150
SCIRun::TetVolMesh::faces_
face_ct faces_
Definition: TetVolMesh.h:2566
SCIRun::Piostream::begin_class
virtual int begin_class(const std::string &name, int current_version)
Definition: Persistent.cc:143
SCIRun::TetVolMesh::orient
void orient(typename Cell::index_type ci)
Orient a cell properly.
Definition: TetVolMesh.h:4359
SCIRun::TetVolMesh::is_editable
virtual bool is_editable() const
Check whether mesh can be altered by adding nodes or elements.
Definition: TetVolMesh.h:399
SCIRun::TetVolMesh::inscribed_circumscribed_radius_metric
double inscribed_circumscribed_radius_metric(INDEX idx) const
Definition: TetVolMesh.h:839
SCIRun::TetVolMesh::get_elem_neighbors
void get_elem_neighbors(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1893
SCIRun::TetVolMesh::get_edges
void get_edges(typename Edge::array_type &array, typename Edge::index_type idx) const
Definition: TetVolMesh.h:466
SCIRun::TetVolMesh::get_canonical_transform
virtual void get_canonical_transform(Core::Geometry::Transform &t) const
Definition: TetVolMesh.h:2859
SCIRun::TetVolMesh::get_elems
void get_elems(typename Elem::array_type &array, typename Cell::index_type idx) const
Definition: TetVolMesh.h:512
SCIRun::Core::Geometry::Vector
Definition: Vector.h:63
SCIRun::TetVolMesh::Face::iterator
FaceIterator< under_type > iterator
Definition: TetVolMesh.h:170
SCIRun::TetVolMesh::Cell::array_type
std::vector< index_type > array_type
Definition: TetVolMesh.h:179
SCIRun::TetVolMesh::Edge
Definition: TetVolMesh.h:161
SCIRun::TetVolMesh::get_node_neighbors
void get_node_neighbors(ARRAY &array, INDEX node) const
We should optimize this function more.
Definition: TetVolMesh.h:1914
SCIRun::TetVolMesh::add_face
void add_face(typename Node::index_type n1, typename Node::index_type n2, typename Node::index_type n3, index_type combined_index)
Definition: TetVolMesh.h:3052
SCIRun::TetVolMesh::get_edge_center
void get_edge_center(Core::Geometry::Point &p, INDEX idx) const
Definition: TetVolMesh.h:2200
SCIRun::Core::Thread::UniqueLock
boost::unique_lock< boost::mutex > UniqueLock
Definition: ConditionVariable.h:43
SCIRun::TetVolMesh::to_index
void to_index(typename Node::index_type &index, index_type i) const
Definition: TetVolMesh.h:440
SCIRun::TetVolMesh::ElemData::edge4_index
index_type edge4_index() const
Definition: TetVolMesh.h:254
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::Core::Geometry::Point::z
void z(const double)
Definition: Point.h:145
SCIRun::Mesh::EDGES_E
Definition: Mesh.h:72
SCIRun::TetVolMesh::elem_reserve
void elem_reserve(size_type s)
Definition: TetVolMesh.h:704
SCIRun::TetVolMesh::get_cell_center
void get_cell_center(Core::Geometry::Point &p, INDEX idx) const
Definition: TetVolMesh.h:2227
SCIRun::TetVolMesh::VTetVolMesh
friend class VTetVolMesh
Make sure the virtual interface has access.
Definition: TetVolMesh.h:139
SCIRun::TetVolMesh::get_edges_from_cell
void get_edges_from_cell(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1595
SCIRun::TetVolMesh::remove_node_from_grid
void remove_node_from_grid(typename Node::index_type ci)
Definition: TetVolMesh.h:3765
SCIRun::TetVolMesh::create_cell_faces
void create_cell_faces(typename Cell::index_type)
Definition: TetVolMesh.h:3513
SCIRun::TetVolMesh::edge_ht
std::map< PEdge, typename Edge::index_type, EdgeHash > edge_ht
Definition: TetVolMesh.h:2555
SCIRun::TetVolMesh::get_delems
void get_delems(typename DElem::array_type &array, typename Edge::index_type idx) const
Definition: TetVolMesh.h:519
SCIRun::TetVolMesh::get_edges_from_face
void get_edges_from_face(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1556
SCIRun::TetVolMesh::insert_node_in_elem
bool insert_node_in_elem(typename Elem::array_type &tets, typename Node::index_type &ni, typename Elem::index_type ci, const Core::Geometry::Point &p)
Definition: TetVolMesh.h:4198
SCIRun::TetVolMesh::get_cells_from_node
void get_cells_from_node(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1686
SCIRun::TetVolMesh::find_closest_node
bool find_closest_node(double &pdist, Core::Geometry::Point &result, INDEX &node, const Core::Geometry::Point &p, double maxdist) const
Definition: TetVolMesh.h:895
SCIRun::TetVolMesh::add_node
Node::index_type add_node(const Core::Geometry::Point &p)
Definition: TetVolMesh.h:686
SCIRun::TetVolMesh::edges_
edge_ct edges_
Definition: TetVolMesh.h:2570
SCIRun::TetVolMesh::get_bounding_box
virtual Core::Geometry::BBox get_bounding_box() const
Get the bounding box of the field.
Definition: TetVolMesh.h:2843
SCIRun::TetVolMesh::face_nt
std::map< PFaceNode, typename Face::index_type, FaceHash > face_nt
Definition: TetVolMesh.h:2554
SCIRun::Mesh::BOUNDING_BOX_E
Definition: Mesh.h:87
SCIRun::TetVolMesh::unsynchronize
virtual bool unsynchronize(mask_type mask)
Definition: TetVolMesh.h:3292
SCIRun::TetVolMesh::set_nodes_by_elem
void set_nodes_by_elem(ARRAY &array, INDEX idx)
Definition: TetVolMesh.h:1843
SCIRun::VMeshShared
Definition: VMeshShared.h:40
SCIRun::TetVolMesh::FaceHash::mid_mask
static const int mid_mask
Definition: TetVolMesh.h:2483
SCIRun::TetVolMesh::cells_
std::vector< under_type > cells_
each 4 indicies make up a tet
Definition: TetVolMesh.h:2315
max
Definition: ParallelLinearAlgebraTests.cc:358
SCIRun::TetVolMesh::ElemData::ElemData
ElemData(const TetVolMesh< Basis > &msh, const index_type ind)
Definition: TetVolMesh.h:199
SCIRun::TetVolMesh::FaceHash::min_buckets
static const size_t min_buckets
Definition: TetVolMesh.h:2478
SCIRun::TetVolMesh::EdgeHash::sz_int
static const int sz_int
These are for our own use (making the hash function.
Definition: TetVolMesh.h:2514
SCIRun::SearchGridT
Definition: SearchGridT.h:47
SCIRun::TetVolMesh::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
Find the closest element to a point.
Definition: TetVolMesh.h:1144
SCIRun::TetVolMesh::remove_face
void remove_face(typename Node::index_type n1, typename Node::index_type n2, typename Node::index_type n3, typename Cell::index_type ci, bool table_only=false)
Definition: TetVolMesh.h:2911
name
const char * name[]
Definition: BoostGraphExampleTests.cc:87
SCIRun::Core::Geometry::Point::x
void x(const double)
Definition: Point.h:125
SCIRun::TetVolMesh::get_faces_from_edge
void get_faces_from_edge(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1736
SCIRun::TetVolMesh
Definition: TetVolMesh.h:77
SCIRun::Pio_index
void Pio_index(Piostream &stream, index_type *data, size_type size)
Definition: Persistent.cc:606
SCIRun::size_type
long long size_type
Definition: Types.h:40
SCIRun::TetVolMesh::delete_nodes
void delete_nodes(std::set< index_type > &to_delete)
Definition: TetVolMesh.h:4021
SCIRun::TetVolMesh::size
void size(typename Node::size_type &) const
Get the iteration sizes.
Definition: TetVolMesh.h:3343
SCIRun::TetVolMesh::scaled_jacobian_metric
double scaled_jacobian_metric(INDEX idx) const
Definition: TetVolMesh.h:781
SCIRun::TetVolMesh::get_faces
void get_faces(typename Face::array_type &array, typename Face::index_type idx) const
Definition: TetVolMesh.h:483
TetQuadraticLgn.h
SCIRun::TetVolMesh::add_point
Node::index_type add_point(const Core::Geometry::Point &p)
Add a new node to the mesh.
Definition: TetVolMesh.h:3916
SCIRun::StackVector
Definition: StackVector.h:50
SCIRun::TetVolMesh::resize_elems
void resize_elems(size_type s)
Definition: TetVolMesh.h:706
SCIRun::TetVolMesh::dimensionality
virtual int dimensionality() const
Topological dimension.
Definition: TetVolMesh.h:380
SCIRun::TetVolMesh::get_epsilon
double get_epsilon() const
Definition: TetVolMesh.h:407
SCIRun::Core::Geometry::Cross
Vector Cross(const Vector &v1, const Vector &v2)
Definition: Vector.h:378
SCIRun::TetVolMesh::size_type
SCIRun::size_type size_type
Definition: TetVolMesh.h:147
SCIRun::TetVolMesh::EdgeHash::up_mask
static const int up_mask
Definition: TetVolMesh.h:2516
SCIRun::TetVolMesh::vmesh_
boost::shared_ptr< VMesh > vmesh_
Pointer to virtual interface.
Definition: TetVolMesh.h:2628
SCIRun::Mesh::NODES_E
Definition: Mesh.h:71
SCIRun::TetVolMesh::interpolate
void interpolate(Core::Geometry::Point &pt, const VECTOR &coords, INDEX idx) const
Definition: TetVolMesh.h:721
SCIRun::TetVolMesh::edge_ct
std::vector< PEdgeCell > edge_ct
Definition: TetVolMesh.h:2560
SCIRun::TetVolMesh::Cell::iterator
CellIterator< under_type > iterator
Definition: TetVolMesh.h:177
SCIRun::Core::Geometry::BBox::diagonal
Vector diagonal() const
Definition: BBox.h:198
SCIRun::Mesh::ELEM_NEIGHBORS_E
Definition: Mesh.h:81
SCIRun::TetVolMesh::det_jacobian
double det_jacobian(const VECTOR &coords, INDEX idx) const
Definition: TetVolMesh.h:740
SCIRun::TetVolMesh::FaceHash::operator()
size_t operator()(const PFACE &f) const
This is the hash function.
Definition: TetVolMesh.h:2488
SCIRun::TetVolMesh::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: TetVolMesh.h:606
Locate.h
SCIRun::TetVolMesh::edge_nt
std::map< PEdgeNode, typename Edge::index_type, EdgeHash > edge_nt
Definition: TetVolMesh.h:2556
SCIRun::TetVolMesh::insert_node_in_edge
void insert_node_in_edge(typename Cell::array_type &tets, typename Node::index_type ni, typename Cell::index_type ci, const PEdgeNode &e)
Definition: TetVolMesh.h:4130
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::PersistentTypeID::type
std::string type
Definition: Persistent.h:72
SCIRun::TetVolMesh::Node::array_type
StackVector< index_type, 4 > array_type
Definition: TetVolMesh.h:158
SCIRun::TetVolMesh::node_type_description
static const TypeDescription * node_type_description()
Definition: TetVolMesh.h:4435
SCIRun::TetVolMesh::to_index
void to_index(typename Cell::index_type &index, index_type i) const
Definition: TetVolMesh.h:446
SCIRun::TetVolMesh::compute_node_grid
void compute_node_grid()
Definition: TetVolMesh.h:3808
CompGeom.h
SCIRun::TetVolMesh::EdgeHash::sz_half_int
static const int sz_half_int
Definition: TetVolMesh.h:2515
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::TetVolMesh::PFace::PFace
PFace(typename Node::index_type n1, typename Node::index_type n2, typename Node::index_type n3)
Definition: TetVolMesh.h:2404
SCIRun::TetVolMesh::fill_neighbors
void fill_neighbors(Iter begin, Iter end, Functor fill_ftor)
SCIRun::Core::Geometry::Transform
Definition: Transform.h:53
SCIRun::TetVolMesh::Face
Definition: TetVolMesh.h:168
SCIRun::TetVolMesh::get_edges
void get_edges(typename Edge::array_type &array, typename Node::index_type idx) const
Definition: TetVolMesh.h:463
SCIRun::TetVolMesh::type_name
static const std::string type_name(int n=-1)
Core functionality for getting the name of a templated mesh class.
Definition: TetVolMesh.h:2777
SCIRun::TetVolMesh::get_elems
void get_elems(typename Elem::array_type &array, typename Face::index_type idx) const
Definition: TetVolMesh.h:509
SCIRun::TetVolMesh::basis_type
Basis basis_type
Definition: TetVolMesh.h:151
SCIRun::Core::Geometry::Vector::x
void x(double)
Definition: Vector.h:175
SCIRun::TetVolMesh::EdgeHash::operator()
bool operator()(const PEDGE &e1, const PEDGE &e2) const
This should return less than rather than equal to.
Definition: TetVolMesh.h:2529
SCIRun::TetVolMesh::get_elem_neighbors
bool get_elem_neighbors(ARRAY &array, INDEX1 elem, INDEX2 delem) const
Definition: TetVolMesh.h:1870
SCIRun::Core::Geometry::closest_point_on_tri
void closest_point_on_tri(Point &result, const Point &orig, const Point &p0, const Point &p1, const Point &p2, const double epsilon)
Definition: CompGeom.cc:136
SCIRun::TetVolMesh::get_length
double get_length(typename Edge::index_type idx) const
More specific names for get_size.
Definition: TetVolMesh.h:593
SCIRun::TetVolMesh::epsilon_
double epsilon_
Definition: TetVolMesh.h:2623
SCIRun::TetVolMesh::PFaceCell
Face information.
Definition: TetVolMesh.h:2318
SCIRun::TetVolMesh::fill_cells
void fill_cells(Iter begin, Iter end, Functor fill_ftor)
Definition: TetVolMesh.h:2749
SCIRun::TetVolMesh::synchronize
virtual bool synchronize(mask_type mask)
Definition: TetVolMesh.h:3173
SCIRun::TetVolMesh::PEdgeNode::nodes_
Node::index_type nodes_[2]
Definition: TetVolMesh.h:2412
SCIRun::TetVolMesh::tetvolmesh_typeid
static PersistentTypeID tetvolmesh_typeid
This ID is created as soon as this class will be instantiated.
Definition: TetVolMesh.h:1385
SCIRun::TetVolMesh::locate
bool locate(typename Cell::index_type &cell, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:630
SCIRun::TetVolMesh::get_center
void get_center(Core::Geometry::Point &result, typename Edge::index_type idx) const
Definition: TetVolMesh.h:554
SCIRun::Mesh::FACES_E
Definition: Mesh.h:73
SCIRun::TetVolMesh::synchronize_cond_
Core::Thread::ConditionVariable synchronize_cond_
Definition: TetVolMesh.h:2614
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::TetVolMesh::ElemData::edge2_index
index_type edge2_index() const
Definition: TetVolMesh.h:244
readAllFields.v
v
Definition: readAllFields.py:42
SCIRun::TetVolMesh::PEdgeCell::cells_
std::vector< index_type > cells_
list of all the cells this edge is in.
Definition: TetVolMesh.h:2455
SCIRun::TetVolMesh::end
void end(typename Node::iterator &) const
Definition: TetVolMesh.h:3334
SCIRun::Core::Geometry::Transform::pre_translate
void pre_translate(const Vector &)
Definition: Transform.cc:278
SCIRun::TetVolMesh::synchronized_
mask_type synchronized_
Definition: TetVolMesh.h:2617
SCIRun::Core::Geometry::Vector::normalize
double normalize()
Definition: Vector.h:437
SCIRun::TetVolMesh::elem_grid_
boost::shared_ptr< SearchGridT< index_type > > elem_grid_
Definition: TetVolMesh.h:2610
SCIRun::TetVolMesh::jacobian_metric
double jacobian_metric(INDEX idx) const
Definition: TetVolMesh.h:816
SCIRun::TetVolMesh::face_ht
std::map< PFace, typename Face::index_type, FaceHash > face_ht
Definition: TetVolMesh.h:2553
share.h
SCIRun::TetVolMesh::get_cells_from_face
void get_cells_from_face(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1825
SCIRun::TetVolMesh::face_table_
face_nt face_table_
Definition: TetVolMesh.h:2567
SCIRun::TetVolMesh::Synchronize::run
void run()
Definition: TetVolMesh.h:307
SCIRun::Core::Geometry::Vector::y
void y(double)
Definition: Vector.h:185
SCIRun::TetVolMesh::get_faces_from_node
void get_faces_from_node(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1772
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::TetVolMesh::compute_faces
void compute_faces()
Definition: TetVolMesh.h:2997
SCIRun::TetVolMesh::FaceHash
Definition: TetVolMesh.h:2471
SCIRun::EdgeIterator
Distinct type for edge Iterator.
Definition: FieldIterator.h:105
SCIRun::TetVolMesh::locate
bool locate(typename Elem::index_type &elem, std::vector< double > &coords, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:633
SCIRun::Mesh::ELEMS_E
Definition: Mesh.h:75
ASSERTFAIL
#define ASSERTFAIL(string)
Definition: Assert.h:52
Mesh.h
SCIRun::TetVolMesh::clone
virtual TetVolMesh * clone() const
Definition: TetVolMesh.h:362
SCIRun::Pio
void Pio(Piostream &stream, Array1< T > &array)
Definition: Array1.h:65
FieldIterator.h
Some convenient simple iterators for fields.
SCIRun::TetVolMesh::get_neighbors
bool get_neighbors(std::vector< typename Elem::index_type > &array, typename Elem::index_type elem, typename DElem::index_type delem) const
Definition: TetVolMesh.h:615
SCIRun::TetVolMesh::ElemData::node1_index
index_type node1_index() const
Definition: TetVolMesh.h:217
SCIRun::TetVolMesh::locate
bool locate(typename Node::index_type &node, const Core::Geometry::Point &p) const
return false if point is out of range.
Definition: TetVolMesh.h:624
SCIRun::TetVolMesh::PFaceCell::cells_
index_type cells_[2]
Definition: TetVolMesh.h:2328
Mutex.h
SCIRun::TetVolMesh::get_nodes
void get_nodes(typename Node::array_type &array, typename Node::index_type idx) const
Get the child topology elements of the given topology.
Definition: TetVolMesh.h:450
SCIRun::TetVolMesh::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: TetVolMesh.h:552
SCIRun::TetVolMesh::get_center
void get_center(Core::Geometry::Point &result, typename Face::index_type idx) const
Definition: TetVolMesh.h:556
SCIRun::TetVolMesh::Cell
Definition: TetVolMesh.h:175
SCIRun::TetVolMesh::find_closest_node
bool find_closest_node(double &pdist, Core::Geometry::Point &result, INDEX &node, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:888
SCIRun::TetVolMesh::delete_cell_faces
void delete_cell_faces(typename Cell::index_type, bool table_only=false)
Definition: TetVolMesh.h:3527
CheckSum.h
SCIRun::TetVolMesh::get_size
double get_size(typename Edge::index_type idx) const
Definition: TetVolMesh.h:565
SCIRun::Mesh::NODE_NEIGHBORS_E
Definition: Mesh.h:80
SCIRun::Piostream::end_class
virtual void end_class()
Definition: Persistent.cc:178
SCIRun::TetVolMesh::Face::array_type
std::vector< index_type > array_type
Definition: TetVolMesh.h:172
SCIRun::TetVolMesh::find_closest_nodes
bool find_closest_nodes(ARRAY &nodes, const Core::Geometry::Point &p, double maxdist) const
Definition: TetVolMesh.h:1000
SCIRun::TetVolMesh::Face::index_type
FaceIndex< under_type > index_type
Definition: TetVolMesh.h:169
SCIRun::TetVolMesh::get_cells_from_edge
void get_cells_from_edge(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1816
SCIRun::TetVolMesh::Edge::iterator
EdgeIterator< under_type > iterator
Definition: TetVolMesh.h:163
SCIRun::TetVolMesh::maker
static Persistent * maker()
This function returns a maker for Pio.
Definition: TetVolMesh.h:1402
SCIRun::TetVolMesh::ElemData::node3
const Core::Geometry::Point & node3() const
Definition: TetVolMesh.h:280
SCIRun::TetVolMesh::ElemData::node0
const Core::Geometry::Point & node0() const
Definition: TetVolMesh.h:265
SCIRun::TetVolMesh::Edge::array_type
std::vector< index_type > array_type
Definition: TetVolMesh.h:165
SCIRun::TetVolMesh::get_elems
void get_elems(typename Elem::array_type &array, typename Node::index_type idx) const
Definition: TetVolMesh.h:503
SCIRun::TetVolMesh::mesh_maker
static MeshHandle mesh_maker()
This function returns a handle for the virtual interface.
Definition: TetVolMesh.h:1404
SCIRun::TetVolMesh::PEdgeCell
Definition: TetVolMesh.h:2452
SCIRun::index_type
long long index_type
Definition: Types.h:39
ConditionVariable.h
SCIRun::TetVolMesh::Edge::index_type
EdgeIndex< under_type > index_type
Definition: TetVolMesh.h:162
SCIRun::TetVolMesh::get_edges
void get_edges(typename Edge::array_type &array, typename Face::index_type idx) const
Definition: TetVolMesh.h:469
SCIRun::Mesh::NODE_LOCATE_E
Definition: Mesh.h:83
SCIRun::TetVolMesh::DElem
Face DElem
Definition: TetVolMesh.h:186
SCIRun::StackVector::resize
void resize(size_t size, const value_type &val=value_type())
Definition: StackVector.h:61
SCIRun::Mesh::DELEMS_E
Definition: Mesh.h:76
SCIRun::TetVolMesh::PFaceNode
Definition: TetVolMesh.h:2331
SCIRun::TetVolMesh::get_nodes_from_face
void get_nodes_from_face(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1518
SearchGridT.h
SCIRun::TetVolMesh::pwl_approx_face
void pwl_approx_face(std::vector< std::vector< VECTOR > > &coords, INDEX ci, unsigned which_face, unsigned div_per_unit) const
Definition: TetVolMesh.h:543
SCIRun::TetVolMesh::delete_cell_syncinfo_special
void delete_cell_syncinfo_special(typename Cell::index_type ci)
Definition: TetVolMesh.h:3640
SCIRun::TetVolMesh::points_
std::vector< Core::Geometry::Point > points_
all the nodes.
Definition: TetVolMesh.h:2312
SCIRun::TetVolMesh::get_cells
void get_cells(typename Cell::array_type &array, typename Cell::index_type idx) const
Definition: TetVolMesh.h:499
SCIRun::TetVolMesh::get_type_description
virtual const TypeDescription * get_type_description() const
Definition: TetVolMesh.h:4428
SCIRun::TetVolMesh::Synchronize::operator()
void operator()()
Definition: TetVolMesh.h:302
SCIRun::TetVolMesh::locate_elem
bool locate_elem(INDEX &elem, ARRAY &coords, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:2149
SCIRun::TetVolMesh::create_cell_node_neighbors
void create_cell_node_neighbors(typename Cell::index_type)
Definition: TetVolMesh.h:3541
SCIRun::TetVolMesh::fill_points
void fill_points(Iter begin, Iter end, Functor fill_ftor)
Definition: TetVolMesh.h:2732
SCIRun::TetVolMesh::get_delems
void get_delems(typename DElem::array_type &array, typename Face::index_type idx) const
Definition: TetVolMesh.h:522
BBox.h
SCIRun::Core::Geometry::BBox::reset
void reset()
Definition: BBox.h:95
SCIRun::TetVolMesh::PFaceCell::shared
bool shared() const
Definition: TetVolMesh.h:2326
SCIRun::TetVolMesh::ElemData::node1
const Core::Geometry::Point & node1() const
Definition: TetVolMesh.h:270
SCIRun::TetVolMesh::edge_table_
edge_nt edge_table_
Definition: TetVolMesh.h:2571
SCIRun::TetVolMesh::bbox_
Core::Geometry::BBox bbox_
Definition: TetVolMesh.h:2621
SCIRun::PersistentTypeID
Definition: Persistent.h:64
SCIRun::TetVolMesh::get_faces_from_cell
void get_faces_from_cell(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1656
SCIRun::TetVolMesh::ElemData::edge5_index
index_type edge5_index() const
Definition: TetVolMesh.h:259
SCIRun::TetVolMesh::to_index
void to_index(typename Edge::index_type &index, index_type i) const
Definition: TetVolMesh.h:442
SCIRun::TetVolMesh::find_closest_elem
bool find_closest_elem(double &pdist, Core::Geometry::Point &result, ARRAY &coords, INDEX &elem, const Core::Geometry::Point &p) const
Find the closest element to a point.
Definition: TetVolMesh.h:1132
SCIRun::TetVolMesh::locate_node
bool locate_node(INDEX &node, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:1936
SCIRun::TetVolMesh::epsilon2_
double epsilon2_
Definition: TetVolMesh.h:2624
SCIRun::TetVolMesh::fill_data
void fill_data(Iter begin, Iter end, Functor fill_ftor)
SCIRun::TetVolMesh::getFacade
MeshFacadeHandle getFacade() const
Definition: TetVolMesh.h:370
SCIRun::EdgeIndex
Distinct type for edge index.
Definition: FieldIndex.h:81
SCIRun::TetVolMesh::PFaceNode::nodes_
Node::index_type nodes_[3]
Definition: TetVolMesh.h:2334
eab.n
int n
Definition: eab.py:9
SCIRun::TetVolMesh::get_weights
int get_weights(const Core::Geometry::Point &p, typename Node::array_type &l, double *w)
Definition: TetVolMesh.h:3705
SCIRun::TetVolMesh::synchronize_lock_
Core::Thread::Mutex synchronize_lock_
Definition: TetVolMesh.h:2613
SCIRun::TetVolMesh::locate_elem
bool locate_elem(INDEX &elem, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:2081
SCIRun::TetVolMesh::get_size
double get_size(typename Face::index_type idx) const
Definition: TetVolMesh.h:575
SCIRun::TetVolMesh::boundary_faces_
std::vector< unsigned char > boundary_faces_
Definition: TetVolMesh.h:2602
SCIRun::TetVolMesh::find_closest_elems
bool find_closest_elems(double &, Core::Geometry::Point &, ARRAY &, const Core::Geometry::Point &) const
Find the closest elements to a point.
Definition: TetVolMesh.h:1371
SCIRun::TetVolMesh::PFace
Definition: TetVolMesh.h:2401
SCIRun::TetVolMesh::compute_edges
void compute_edges()
Definition: TetVolMesh.h:3111
SCIRun::TetVolMesh::get_weights
int get_weights(const Core::Geometry::Point &, typename Face::array_type &, double *)
Definition: TetVolMesh.h:643
SCIRun::TetVolMesh::remove_elem_from_grid
void remove_elem_from_grid(typename Elem::index_type ci)
Definition: TetVolMesh.h:3742
SCIRun::TetVolMesh::Synchronize::Synchronize
Synchronize(TetVolMesh< Basis > *mesh, mask_type sync)
Definition: TetVolMesh.h:299
SCIRun::TetVolMesh::compute_elem_grid
void compute_elem_grid()
Definition: TetVolMesh.h:3772
SCIRun::TetVolMesh::EdgeHash::min_buckets
static const size_t min_buckets
Definition: TetVolMesh.h:2511
SCIRun::TetVolMesh::locate
bool locate(typename Face::index_type &face, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:628
SCIRun::TetVolMesh::get_nodes_from_elem
void get_nodes_from_elem(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1550
SCIRun::TetVolMesh::synchronizing_
mask_type synchronizing_
Definition: TetVolMesh.h:2619
DEBUG_CONSTRUCTOR
#define DEBUG_CONSTRUCTOR(type)
Definition: Debug.h:64
SCIRun::TetVolMesh::delete_cell_edges
void delete_cell_edges(typename Cell::index_type, bool table_only=false)
Definition: TetVolMesh.h:3498
SCIRun::TetVolMesh::get_neighbors
void get_neighbors(std::vector< typename Node::index_type > &array, typename Node::index_type node) const
These are more general implementations.
Definition: TetVolMesh.h:612
SCIRun::TetVolMesh::find_closest_elem
bool find_closest_elem(double &pdist, Core::Geometry::Point &result, INDEX &elem, const Core::Geometry::Point &p) const
Definition: TetVolMesh.h:1362
SCIRun::Mesh::io
void io(Piostream &stream)
Persistent I/O.
Definition: Mesh.cc:387
SCIRun::TetVolMesh::delete_cells
void delete_cells(std::set< index_type > &to_delete)
Deleting cells, make sure mesh is detached.
Definition: TetVolMesh.h:3989
SCIRun::TetVolMesh::get_elems
void get_elems(typename Elem::array_type &array, typename Edge::index_type idx) const
Definition: TetVolMesh.h:506
SCIRun::TetVolMesh::dynamic_type_name
virtual std::string dynamic_type_name() const
Definition: TetVolMesh.h:1389
SCIRun::TetVolMesh::hash_edge
void hash_edge(typename Node::index_type n1, typename Node::index_type n2, index_type combined_index, edge_ht &table)
Definition: TetVolMesh.h:3086
SCIRun::TetVolMesh::io
virtual void io(Piostream &)
Export this class using the old Pio system.
Definition: TetVolMesh.h:4380
SCIRun::TetVolMesh::has_normals
virtual bool has_normals() const
Has this mesh normals.
Definition: TetVolMesh.h:402
SCIRun::TetVolMesh::get_size
double get_size(typename Elem::index_type idx) const
Definition: TetVolMesh.h:586
SCIRun::TetVolMesh::get_nodes_from_cell
void get_nodes_from_cell(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1539
SCIRun::TetVolMesh::sz_int
static const int sz_int
Definition: TetVolMesh.h:2470
SCIRun::TetVolMesh::add_edge
void add_edge(typename Node::index_type n1, typename Node::index_type n2, index_type combined_index)
Definition: TetVolMesh.h:3152
SCIRun::TetVolMesh::FaceHash::bucket_size
static const size_t bucket_size
Definition: TetVolMesh.h:2477
SCIRun::Core::Geometry::Vector::z
void z(double)
Definition: Vector.h:195
SCIRun::Core::Datatypes::MeshTraits< VMesh >::MeshFacadeHandle
boost::shared_ptr< MeshFacade< VMesh > > MeshFacadeHandle
Definition: MeshTraits.h:61
SCIRun::VMesh
Definition: VMesh.h:53
SCIRun::TetVolMesh::add_elem
Elem::index_type add_elem(ARRAY a)
Add a new element to the mesh.
Definition: TetVolMesh.h:691
SCIRun::TetVolMesh::get_random_point
void get_random_point(Core::Geometry::Point &p, typename Elem::index_type i, FieldRNG &r) const
Definition: TetVolMesh.h:2801
SCIRun::TetVolMesh::Node
Index and Iterator types required for Mesh Concept.
Definition: TetVolMesh.h:154
SCIRun::TetVolMesh::has_face_normals
virtual bool has_face_normals() const
Has this mesh face normals.
Definition: TetVolMesh.h:405
SCIRun::TetVolMesh::node_grid_
boost::shared_ptr< SearchGridT< index_type > > node_grid_
Definition: TetVolMesh.h:2609
SCIRun::TetVolMesh::PFaceCell::PFaceCell
PFaceCell()
Definition: TetVolMesh.h:2320
SCIRun::TetVolMesh::edge_type_description
static const TypeDescription * edge_type_description()
Definition: TetVolMesh.h:4452
SCIRun::TetVolMesh::PFaceNode::operator==
bool operator==(const PFaceNode &f) const
true if both have the same nodes (order does not matter)
Definition: TetVolMesh.h:2381
SCIRun::TetVolMesh::under_type
SCIRun::index_type under_type
Definition: TetVolMesh.h:145
VirtualMeshFacade.h
MiscMath.h
SCIRun::TetVolMesh::hash_face
void hash_face(typename Node::index_type n1, typename Node::index_type n2, typename Node::index_type n3, index_type ci, face_ht &table)
Definition: TetVolMesh.h:2957
SCIRun::TetVolMesh::Cell::size_type
CellIndex< under_type > size_type
Definition: TetVolMesh.h:178
SCIRun::TetVolMesh::get_delems
void get_delems(typename DElem::array_type &array, typename Node::index_type idx) const
Definition: TetVolMesh.h:516
SCIRun::TetVolMesh::derivate
void derivate(const VECTOR1 &coords, INDEX idx, VECTOR2 &J) const
Definition: TetVolMesh.h:731
SCIRun::TetVolMesh::transform
virtual void transform(const Core::Geometry::Transform &t)
Core::Geometry::Transform a field (transform all nodes using this transformation matrix) ...
Definition: TetVolMesh.h:2869
SCIRun::TetVolMesh::get_delems
void get_delems(typename DElem::array_type &array, typename Cell::index_type idx) const
Definition: TetVolMesh.h:525
eabLatVolData.size
int size
Definition: eabLatVolData.py:2
SCIRun::TetVolMesh::basis_
Basis basis_
Definition: TetVolMesh.h:2622
SCIRun::TetVolMesh::get_edges_from_elem
void get_edges_from_elem(ARRAY &array, INDEX idx) const
Definition: TetVolMesh.h:1650
SCIRun::TetVolMesh::to_index
void to_index(typename Face::index_type &index, index_type i) const
Definition: TetVolMesh.h:444
SCIRun::TetVolMesh::add_find_point
Node::index_type add_find_point(const Core::Geometry::Point &p, double err=1.0e-6)
Definition: TetVolMesh.h:3879
SCIRun::TetVolMesh::get_normal
void get_normal(Core::Geometry::Vector &, typename Node::index_type) const
Normals for visualizations.
Definition: TetVolMesh.h:655
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::elem_locate
bool elem_locate(INDEX &elem, MESH &msh, const Core::Geometry::Point &p)
General case locate, search each elem.
Definition: Mesh.h:188
SCIRun::TetVolMesh::get_cells
void get_cells(typename Cell::array_type &array, typename Face::index_type idx) const
Definition: TetVolMesh.h:496
Generated on Wed Jul 9 2014 14:31:46 for SCIRun by   doxygen 1.8.5