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