SCIRun  5.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
StructQuadSurfMesh.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 ///@details
30 /// A structured grid is a dataset with regular topology but with irregular
31 /// geometry. The grid may have any shape but can not be overlapping or
32 /// self-intersecting. The topology of a structured grid is represented
33 /// using a 2D, or 3D vector with the points being stored in an index
34 /// based array. The faces (quadrilaterals) and cells (Hexahedron) are
35 /// implicitly define based based upon their indexing.
36 ///
37 /// Structured grids are typically used in finite difference analysis.
38 ///
39 /// For more information on datatypes see Schroeder, Martin, and Lorensen,
40 /// "The Visualization Toolkit", Prentice Hall, 1998.
41 ///
42 
43 
44 #ifndef CORE_DATATYPES_STRUCTQUADSURFMESH_H
45 #define CORE_DATATYPES_STRUCTQUADSURFMESH_H 1
46 
47 /// Include what kind of support we want to have
48 /// Need to fix this and couple it to sci-defs
49 #include <Core/Datatypes/Legacy/Field/MeshSupport.h>
50 
51 #include <Core/GeometryPrimitives/Point.h>
52 #include <Core/GeometryPrimitives/Plane.h>
53 #include <Core/GeometryPrimitives/CompGeom.h>
54 
55 #include <Core/Containers/Array2.h>
56 #include <Core/GeometryPrimitives/SearchGridT.h>
57 #include <Core/Thread/Mutex.h>
58 
59 #include <Core/Datatypes/Legacy/Field/ImageMesh.h>
60 
61 #include <Core/Datatypes/Legacy/Field/share.h>
62 
63 namespace SCIRun {
64 
65 /// Declarations for virtual interface
66 
67 /// Functions for creating the virtual interface
68 /// Declare the functions that instantiate the virtual interface
69 template <class Basis>
70 class StructQuadSurfMesh;
71 
72 /// make sure any other mesh other than the preinstantiate ones
73 /// returns no virtual interface. Altering this behavior will allow
74 /// for dynamically compiling the interface if needed.
75 template<class MESH>
76 VMesh* CreateVStructQuadSurfMesh(MESH* mesh) { return (0); }
77 
78 /// These declarations are needed for a combined dynamic compilation as
79 /// as well as virtual functions solution.
80 /// Declare that these can be found in a library that is already
81 /// precompiled. So dynamic compilation will not instantiate them again.
82 
83 #if (SCIRUN_STRUCTQUADSURF_SUPPORT > 0)
84 
85 SCISHARE VMesh* CreateVStructQuadSurfMesh(StructQuadSurfMesh<Core::Basis::QuadBilinearLgn<Core::Geometry::Point> >* mesh);
86 
87 #endif
88 
89 
90 template <class Basis>
91 class StructQuadSurfMesh : public ImageMesh<Basis>
92 {
93 
94 /// Make sure the virtual interface has access
95 template <class MESH> friend class VStructQuadSurfMesh;
96 
97 public:
98  /// Types that change depending on 32 or 64bits
99  typedef SCIRun::index_type under_type;
100  typedef SCIRun::index_type index_type;
101  typedef SCIRun::size_type size_type;
102  typedef SCIRun::mask_type mask_type;
103 
104  typedef boost::shared_ptr<StructQuadSurfMesh<Basis> > handle_type;
105 
106  StructQuadSurfMesh();
107  StructQuadSurfMesh(size_type x, size_type y);
108  StructQuadSurfMesh(const StructQuadSurfMesh<Basis> &copy);
109  virtual StructQuadSurfMesh *clone() const
110  { return new StructQuadSurfMesh<Basis>(*this); }
111  virtual ~StructQuadSurfMesh()
112  {
113  DEBUG_DESTRUCTOR("StructQuadSurfMesh")
114  }
115 
116  /// get the mesh statistics
117  virtual Core::Geometry::BBox get_bounding_box() const;
118  virtual void transform(const Core::Geometry::Transform &t);
119 
120  virtual bool get_dim(std::vector<size_type>&) const;
121  virtual void set_dim(std::vector<size_type> dims) {
122  ImageMesh<Basis>::ni_ = dims[0];
123  ImageMesh<Basis>::nj_ = dims[1];
124 
125  points_.resize(dims[1], dims[0]);
126  normals_.resize(dims[1], dims[0]);
127 
128  /// Create a new virtual interface for this copy
129  /// all pointers have changed hence create a new
130  /// virtual interface class
131  ImageMesh<Basis>::vmesh_.reset(CreateVStructQuadSurfMesh(this));
132  }
133 
134  virtual int topology_geometry() const
135  {
136  return (Mesh::STRUCTURED | Mesh::IRREGULAR);
137  }
138 
139  /// Get the size of an element (length, area, volume)
140  double get_size(const typename ImageMesh<Basis>::Node::index_type &) const
141  { return 0.0; }
142  double get_size(typename ImageMesh<Basis>::Edge::index_type idx) const
143  {
144  typename ImageMesh<Basis>::Node::array_type arr;
145  this->get_nodes(arr, idx);
146  Core::Geometry::Point p0, p1;
147  get_center(p0, arr[0]);
148  get_center(p1, arr[1]);
149  return (p1 - p0).length();
150  }
151  double get_size(const typename ImageMesh<Basis>::Face::index_type &idx) const
152  {
153  /// Sum the sizes of the two triangles.
154  const Core::Geometry::Point &p0 = points_(idx.j_ + 0, idx.i_ + 0);
155  const Core::Geometry::Point &p1 = points_(idx.j_ + 0, idx.i_ + 1);
156  const Core::Geometry::Point &p2 = points_(idx.j_ + 1, idx.i_ + 1);
157  const Core::Geometry::Point &p3 = points_(idx.j_ + 1, idx.i_ + 0);
158  const double a0 = Cross(p0 - p1, p2 - p0).length();
159  const double a1 = Cross(p2 - p3, p0 - p2).length();
160  return (a0 + a1) * 0.5;
161  }
162 
163  double get_size(typename ImageMesh<Basis>::Cell::index_type) const
164  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
165  double get_length(typename ImageMesh<Basis>::Edge::index_type idx) const
166  { return get_size(idx); }
167  double get_area(const typename ImageMesh<Basis>::Face::index_type &idx) const
168  { return get_size(idx); }
169  double get_volume(typename ImageMesh<Basis>::Cell::index_type /*idx*/) const
170  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
171 
172  void get_normal(Core::Geometry::Vector &,
173  const typename ImageMesh<Basis>::Node::index_type &) const;
174 
175  void get_normal(Core::Geometry::Vector &result, std::vector<double> &coords,
176  typename ImageMesh<Basis>::Elem::index_type eidx,
177  unsigned int)
178  {
179  ElemData ed(*this, eidx);
180  std::vector<Core::Geometry::Point> Jv;
181  this->basis_.derivate(coords, ed, Jv);
182  result = Cross(Core::Geometry::Vector(Jv[0]), Core::Geometry::Vector(Jv[1]));
183  result.normalize();
184  }
185  /// get the center point (in object space) of an element
186  void get_center(Core::Geometry::Point &,
187  const typename ImageMesh<Basis>::Node::index_type &) const;
188  void get_center(Core::Geometry::Point &,
189  typename ImageMesh<Basis>::Edge::index_type) const;
190  void get_center(Core::Geometry::Point &,
191  const typename ImageMesh<Basis>::Face::index_type &) const;
192  void get_center(Core::Geometry::Point &,
193  typename ImageMesh<Basis>::Cell::index_type) const
194  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
195 
196  bool locate(typename ImageMesh<Basis>::Node::index_type &node,
197  const Core::Geometry::Point &p) const
198  { return (locate_node(node,p)); }
199  bool locate(typename ImageMesh<Basis>::Edge::index_type &edge,
200  const Core::Geometry::Point &p) const
201  { return (locate_edge(edge,p)); }
202  bool locate(typename ImageMesh<Basis>::Face::index_type &face,
203  const Core::Geometry::Point &p ) const
204  { return (locate_elem(face,p)); }
205  bool locate(typename ImageMesh<Basis>::Cell::index_type &,
206  const Core::Geometry::Point &) const
207  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
208 
209  template <class ARRAY>
210  bool locate(typename ImageMesh<Basis>::Elem::index_type &elem,
211  ARRAY& coords,
212  const Core::Geometry::Point &p ) const
213  { return (locate_elem(elem,coords,p)); }
214 
215 
216  bool find_closest_node(double& pdist, Core::Geometry::Point &result,
217  typename ImageMesh<Basis>::Node::index_type &node,
218  const Core::Geometry::Point &p) const
219  { return (find_closest_node(pdist,result,node,p,-1.0)); }
220 
221  bool find_closest_node(double& pdist, Core::Geometry::Point &result,
222  typename ImageMesh<Basis>::Node::index_type &node,
223  const Core::Geometry::Point &p, double maxdist) const;
224 
225  bool find_closest_elem(double& pdist, Core::Geometry::Point &result,
226  typename ImageMesh<Basis>::Elem::index_type &elem,
227  const Core::Geometry::Point &p) const;
228 
229  bool find_closest_elems(double& pdist, Core::Geometry::Point &result,
230  std::vector<typename ImageMesh<Basis>::Elem::index_type> &elems,
231  const Core::Geometry::Point &p) const;
232 
233  int get_weights(const Core::Geometry::Point &p,
234  typename ImageMesh<Basis>::Node::array_type &l, double *w);
235  int get_weights(const Core::Geometry::Point & ,
236  typename ImageMesh<Basis>::Edge::array_type & , double * )
237  { ASSERTFAIL("StructQuadSurfMesh::get_weights for edges isn't supported"); }
238  int get_weights(const Core::Geometry::Point &p,
239  typename ImageMesh<Basis>::Face::array_type &l, double *w);
240  int get_weights(const Core::Geometry::Point & ,
241  typename ImageMesh<Basis>::Cell::array_type & , double * )
242  { ASSERTFAIL("This mesh type does not have cells use \"elem\"."); }
243 
244 
245  bool inside3_p(typename ImageMesh<Basis>::Face::index_type i,
246  const Core::Geometry::Point &p) const;
247 
248 
249  void get_point(Core::Geometry::Point &point,
250  const typename ImageMesh<Basis>::Node::index_type &index) const
251  { get_center(point, index); }
252  void set_point(const Core::Geometry::Point &point,
253  const typename ImageMesh<Basis>::Node::index_type &index);
254 
255  void get_random_point(Core::Geometry::Point &,
256  const typename ImageMesh<Basis>::Elem::index_type &,
257  FieldRNG &rng) const;
258 
259  class ElemData
260  {
261  public:
262 
263  ElemData(const StructQuadSurfMesh<Basis>& msh,
264  const typename ImageMesh<Basis>::Elem::index_type ind) :
265  mesh_(msh),
266  index_(ind)
267  {}
268 
269  /// the following designed to coordinate with ::get_nodes
270  inline
271  index_type node0_index() const
272  {
273  return (index_.i_ + mesh_.get_ni()*index_.j_);
274  }
275  inline
276  index_type node1_index() const
277  {
278  return (index_.i_+ 1 + mesh_.get_ni()*index_.j_);
279  }
280  inline
281  index_type node2_index() const
282  {
283  return (index_.i_ + 1 + mesh_.get_ni()*(index_.j_ + 1));
284 
285  }
286  inline
287  index_type node3_index() const
288  {
289  return (index_.i_ + mesh_.get_ni()*(index_.j_ + 1));
290  }
291 
292  /// the following designed to coordinate with ::get_edges
293  inline
294  index_type edge0_index() const
295  {
296  return index_.i_ + index_.j_ * (mesh_.ni_- 1);
297  }
298  inline
299  index_type edge1_index() const
300  {
301  return index_.i_ + (index_.j_ + 1) * (mesh_.ni_ - 1);
302  }
303  inline
304  index_type edge2_index() const
305  {
306  return index_.i_ *(mesh_.nj_ - 1) + index_.j_ +
307  ((mesh_.ni_ - 1) * mesh_.nj_);
308  }
309  inline
310  index_type edge3_index() const
311  {
312  return (index_.i_ + 1) * (mesh_.nj_ - 1) + index_.j_ +
313  ((mesh_.ni_ - 1) * mesh_.nj_);
314  }
315 
316 
317  inline
318  const Core::Geometry::Point &node0() const
319  {
320  return mesh_.points_(index_.j_, index_.i_);
321  }
322  inline
323  const Core::Geometry::Point &node1() const
324  {
325  return mesh_.points_(index_.j_, index_.i_+1);
326  }
327  inline
328  const Core::Geometry::Point &node2() const
329  {
330  return mesh_.points_(index_.j_+1, index_.i_+1);
331  }
332  inline
333  const Core::Geometry::Point &node3() const
334  {
335  return mesh_.points_(index_.j_+1, index_.i_);
336  }
337  inline
338  const Core::Geometry::Point &node4() const
339  {
340  return mesh_.points_(index_.j_, index_.i_);
341  }
342 
343  private:
344  const StructQuadSurfMesh<Basis> &mesh_;
345  const typename ImageMesh<Basis>::Elem::index_type index_;
346  };
347 
348  friend class ElemData;
349 
350  /// Generate the list of points that make up a sufficiently accurate
351  /// piecewise linear approximation of an edge.
352  void pwl_approx_edge(std::vector<std::vector<double> > &coords,
353  typename ImageMesh<Basis>::Elem::index_type /*ci*/,
354  unsigned which_edge,
355  unsigned div_per_unit) const
356  {
357  /// Needs to match unit_edges in Basis/QuadBilinearLgn.cc compare
358  /// get_nodes order to the basis order
359  int emap[] = {0, 2, 3, 1};
360  this->basis_.approx_edge(emap[which_edge], div_per_unit, coords);
361  }
362 
363  /// Get the local coordinates for a certain point within an element
364  /// This function uses a couple of newton iterations to find the
365  /// local ! coordinate of a point
366  template<class VECTOR>
367  bool get_coords(VECTOR &coords, const Core::Geometry::Point &p,
368  typename ImageMesh<Basis>::Elem::index_type idx) const
369  {
370  ElemData ed(*this, idx);
371  return this->basis_.get_coords(coords, p, ed);
372  }
373 
374  /// Find the location in the global coordinate system for a local
375  /// coordinate ! This function is the opposite of get_coords.
376  template<class VECTOR>
377  void interpolate(Core::Geometry::Point &pt, const VECTOR &coords,
378  typename ImageMesh<Basis>::Elem::index_type idx) const
379  {
380  ElemData ed(*this, idx);
381  pt = this->basis_.interpolate(coords, ed);
382  }
383 
384  /// Interpolate the derivate of the function, This infact will
385  /// return the ! jacobian of the local to global coordinate
386  /// transformation. This function ! is mainly intended for the non
387  /// linear elements
388  template<class VECTOR1, class VECTOR2>
389  void derivate(const VECTOR1 &coords,
390  typename ImageMesh<Basis>::Elem::index_type idx,
391  VECTOR2 &J) const
392  {
393  ElemData ed(*this, idx);
394  this->basis_.derivate(coords, ed, J);
395  }
396 
397  /// Get the determinant of the jacobian, which is the local volume
398  /// of an element ! and is intended to help with the integration of
399  /// functions over an element.
400  template<class VECTOR>
401  double det_jacobian(const VECTOR& coords,
402  typename ImageMesh<Basis>::Elem::index_type idx) const
403  {
404  double J[9];
405  jacobian(coords,idx,J);
406  return (DetMatrix3x3(J));
407  }
408 
409  /// Get the jacobian of the transformation. In case one wants the
410  /// non inverted ! version of this matrix. This is currentl here for
411  /// completeness of the ! interface
412  template<class VECTOR>
413  void jacobian(const VECTOR& coords,
414  typename ImageMesh<Basis>::Elem::index_type idx,
415  double* J) const
416  {
417  StackVector<Core::Geometry::Point,2> Jv;
418  ElemData ed(*this,idx);
419  this->basis_.derivate(coords,ed,Jv);
420  Core::Geometry::Vector Jv2 = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1]));
421  Jv2.normalize();
422  J[0] = Jv[0].x();
423  J[1] = Jv[0].y();
424  J[2] = Jv[0].z();
425  J[3] = Jv[1].x();
426  J[4] = Jv[1].y();
427  J[5] = Jv[1].z();
428  J[6] = Jv2.x();
429  J[7] = Jv2.y();
430  J[8] = Jv2.z();
431  }
432 
433  /// Get the inverse jacobian of the transformation. This one is
434  /// needed to ! translate local gradients into global gradients. Hence
435  /// it is crucial for ! calculating gradients of fields, or
436  /// constructing finite elements.
437  template<class VECTOR>
438  double inverse_jacobian(const VECTOR& coords,
439  typename ImageMesh<Basis>::Elem::index_type idx,
440  double* Ji) const
441  {
442  StackVector<Core::Geometry::Point,2> Jv;
443  ElemData ed(*this,idx);
444  this->basis_.derivate(coords,ed,Jv);
445  double J[9];
446  Core::Geometry::Vector Jv2 = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1]));
447  Jv2.normalize();
448  J[0] = Jv[0].x();
449  J[1] = Jv[0].y();
450  J[2] = Jv[0].z();
451  J[3] = Jv[1].x();
452  J[4] = Jv[1].y();
453  J[5] = Jv[1].z();
454  J[6] = Jv2.x();
455  J[7] = Jv2.y();
456  J[8] = Jv2.z();
457 
458  return (InverseMatrix3x3(J,Ji));
459  }
460 
461 
462  template<class INDEX>
463  double scaled_jacobian_metric(const INDEX idx) const
464  {
465  StackVector<Core::Geometry::Point,3> Jv;
466  ElemData ed(*this,idx);
467 
468  double temp;
469 
470  this->basis_.derivate(this->basis_.unit_center,ed,Jv);
471  Jv.resize(3);
472  Core::Geometry::Vector v = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1])); v.normalize();
473  Jv[2] = v.asPoint();
474  double min_jacobian = ScaledDetMatrix3P(Jv);
475 
476  size_t num_vertices = this->basis_.number_of_vertices();
477  for (size_t j=0;j < num_vertices;j++)
478  {
479  this->basis_.derivate(this->basis_.unit_vertices[j],ed,Jv);
480  Jv.resize(3);
481  v = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1])); v.normalize();
482  Jv[2] = v.asPoint();
483  temp = ScaledDetMatrix3P(Jv);
484  if(temp < min_jacobian) min_jacobian = temp;
485  }
486 
487  return (min_jacobian);
488  }
489 
490 
491  template<class INDEX>
492  double jacobian_metric(const INDEX idx) const
493  {
494  StackVector<Core::Geometry::Point,3> Jv;
495  ElemData ed(*this,idx);
496 
497  double temp;
498 
499  this->basis_.derivate(this->basis_.unit_center,ed,Jv);
500  Jv.resize(3);
501  Core::Geometry::Vector v = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1])); v.normalize();
502  Jv[2] = v.asPoint();
503  double min_jacobian = DetMatrix3P(Jv);
504 
505  size_t num_vertices = this->basis_.number_of_vertices();
506  for (size_t j=0;j < num_vertices;j++)
507  {
508  this->basis_.derivate(this->basis_.unit_vertices[j],ed,Jv);
509  Jv.resize(3);
510  v = Cross(Core::Geometry::Vector(Jv[0]),Core::Geometry::Vector(Jv[1])); v.normalize();
511  Jv[2] = v.asPoint();
512  temp = DetMatrix3P(Jv);
513  if(temp < min_jacobian) min_jacobian = temp;
514  }
515 
516  return (min_jacobian);
517  }
518 
519 
520  template <class INDEX>
521  bool inside(INDEX idx, const Core::Geometry::Point &p) const
522  {
523  const Core::Geometry::Point &p0 = points_(idx.j_,idx.i_);
524  const Core::Geometry::Point &p1 = points_(idx.j_+1,idx.i_);
525  const Core::Geometry::Point &p2 = points_(idx.j_+1,idx.i_+1);
526  const Core::Geometry::Point &p3 = points_(idx.j_,idx.i_+1);
527  const double x0 = p0.x();
528  const double y0 = p0.y();
529  const double z0 = p0.z();
530  const double x1 = p1.x();
531  const double y1 = p1.y();
532  const double z1 = p1.z();
533  const double x2 = p2.x();
534  const double y2 = p2.y();
535  const double z2 = p2.z();
536  const double x3 = p3.x();
537  const double y3 = p3.y();
538  const double z3 = p3.z();
539 
540  const double a0 = + x1*(y2*z3-y3*z2) + x2*(y3*z1-y1*z3) + x3*(y1*z2-y2*z1);
541  const double a1 = - x2*(y3*z0-y0*z3) - x3*(y0*z2-y2*z0) - x0*(y2*z3-y3*z2);
542  const double a2 = + x3*(y0*z1-y1*z0) + x0*(y1*z3-y3*z1) + x1*(y3*z0-y0*z3);
543  const double a3 = - x0*(y1*z2-y2*z1) - x1*(y2*z0-y0*z2) - x2*(y0*z1-y1*z0);
544  const double iV6 = 1.0 / (a0+a1+a2+a3);
545 
546  const double b0 = - (y2*z3-y3*z2) - (y3*z1-y1*z3) - (y1*z2-y2*z1);
547  const double c0 = + (x2*z3-x3*z2) + (x3*z1-x1*z3) + (x1*z2-x2*z1);
548  const double d0 = - (x2*y3-x3*y2) - (x3*y1-x1*y3) - (x1*y2-x2*y1);
549  const double s0 = iV6 * (a0 + b0*p.x() + c0*p.y() + d0*p.z());
550  if (s0 < -1e-8) return (false);
551 
552  const double b1 = + (y3*z0-y0*z3) + (y0*z2-y2*z0) + (y2*z3-y3*z2);
553  const double c1 = - (x3*z0-x0*z3) - (x0*z2-x2*z0) - (x2*z3-x3*z2);
554  const double d1 = + (x3*y0-x0*y3) + (x0*y2-x2*y0) + (x2*y3-x3*y2);
555  const double s1 = iV6 * (a1 + b1*p.x() + c1*p.y() + d1*p.z());
556  if (s1 < -1e-8) return (false);
557 
558  const double b2 = - (y0*z1-y1*z0) - (y1*z3-y3*z1) - (y3*z0-y0*z3);
559  const double c2 = + (x0*z1-x1*z0) + (x1*z3-x3*z1) + (x3*z0-x0*z3);
560  const double d2 = - (x0*y1-x1*y0) - (x1*y3-x3*y1) - (x3*y0-x0*y3);
561  const double s2 = iV6 * (a2 + b2*p.x() + c2*p.y() + d2*p.z());
562  if (s2 < -1e-8) return (false);
563 
564  const double b3 = +(y1*z2-y2*z1) + (y2*z0-y0*z2) + (y0*z1-y1*z0);
565  const double c3 = -(x1*z2-x2*z1) - (x2*z0-x0*z2) - (x0*z1-x1*z0);
566  const double d3 = +(x1*y2-x2*y1) + (x2*y0-x0*y2) + (x0*y1-x1*y0);
567  const double s3 = iV6 * (a3 + b3*p.x() + c3*p.y() + d3*p.z());
568  if (s3 < -1e-8) return (false);
569 
570  return (true);
571  }
572 
573 
574  template <class INDEX>
575  inline bool locate_node(INDEX &node, const Core::Geometry::Point &p) const
576  {
577  /// If there are no nodes we cannot find a closest point
578  if (this->ni_ == 0 || this->nj_ == 0) return (false);
579 
580  /// Check first guess
581  if (node.i_ >= 0 && node.i_ < this->ni_ &&
582  node.j_ >= 0 && node.j_ < this->nj_)
583  {
584  node.mesh_ = this;
585  if ((p - points_(node.j_,node.i_)).length2() < epsilon2_) return (true);
586  }
587 
588  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
589  "QuadSurfMesh::locate_node requires synchronize(NODE_LOCATE_E).")
590 
591  // get grid sizes
592  const size_type ni = node_grid_->get_ni()-1;
593  const size_type nj = node_grid_->get_nj()-1;
594  const size_type nk = node_grid_->get_nk()-1;
595 
596  // Convert to grid coordinates.
597  index_type bi, bj, bk, ei, ej, ek;
598  node_grid_->unsafe_locate(bi, bj, bk, p);
599 
600  // Clamp to closest point on the grid.
601  if (bi > ni) bi =ni; if (bi < 0) bi = 0;
602  if (bj > nj) bj =nj; if (bj < 0) bj = 0;
603  if (bk > nk) bk =nk; if (bk < 0) bk = 0;
604 
605  ei = bi;
606  ej = bj;
607  ek = bk;
608 
609  double dmin = DBL_MAX;
610  bool found;
611  do
612  {
613  found = true;
614  /// This looks incorrect - but it is correct
615  /// We need to do a full shell without any elements that are closer
616  /// to make sure there no closer elements in neighboring searchgrid cells
617 
618  for (index_type i = bi; i <= ei; i++)
619  {
620  if (i < 0 || i > ni) continue;
621  for (index_type j = bj; j <= ej; j++)
622  {
623  if (j < 0 || j > nj) continue;
624  for (index_type k = bk; k <= ek; k++)
625  {
626  if (k < 0 || k > nk) continue;
627  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
628  {
629  if (node_grid_->min_distance_squared(p, i, j, k) < dmin)
630  {
631  found = false;
632  typename SearchGridT<typename ImageMesh<Basis>::Node::index_type >::iterator it, eit;
633  node_grid_->lookup_ijk(it,eit, i, j, k);
634 
635  while (it != eit)
636  {
637  const Core::Geometry::Point point = points_((*it).j_,(*it).i_);
638  const double dist = (p-point).length2();
639 
640  if (dist < dmin)
641  {
642  node = INDEX(*it);
643  dmin = dist;
644 
645  if (dist < epsilon2_) return (true);
646  }
647 
648  ++it;
649  }
650  }
651  }
652  }
653  }
654  }
655  bi--;ei++;
656  bj--;ej++;
657  bk--;ek++;
658  }
659  while ((!found)||(dmin == DBL_MAX)) ;
660 
661  return (true);
662  }
663 
664  /// This is currently implemented as an exhaustive search
665  template <class INDEX>
666  inline bool locate_edge(INDEX &idx, const Core::Geometry::Point &p) const
667  {
668  ASSERTMSG(synchronized_ & Mesh::EDGES_E,
669  "QuadSurfMesh::locate_node requires synchronize(EDGES_E).")
670 
671  typename ImageMesh<Basis>::Edge::iterator bi, ei;
672  typename ImageMesh<Basis>::Node::array_type nodes;
673  this->begin(bi);
674  this->end(ei);
675  idx = 0;
676 
677  bool found = false;
678  double mindist = 0.0;
679  while (bi != ei)
680  {
681  this->get_nodes(nodes,*bi);
682  const double dist = distance_to_line2(p, points_(nodes[0].j_,nodes[0].i_),
683  points_(nodes[1].j_,nodes[1].i_),epsilon_);
684  if (!found || dist < mindist)
685  {
686  idx = static_cast<INDEX>(*bi);
687  mindist = dist;
688  found = true;
689  }
690  ++bi;
691  }
692  return (found);
693  }
694 
695 
696  template <class INDEX>
697  inline bool locate_elem(INDEX &elem, const Core::Geometry::Point &p) const
698  {
699  if (this->basis_.polynomial_order() > 1) return elem_locate(elem, *this, p);
700 
701  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
702  "StructQuadSurfMesh::locate_node requires synchronize(ELEM_LOCATE_E).")
703 
704  typename ImageMesh<Basis>::Elem::size_type sz;
705  this->size(sz);
706  if ((elem > 0)&&(elem < sz))
707  {
708  elem.mesh_ = this;
709  if (inside(elem,p)) return (true);
710  }
711 
712  if (!elem_grid_) return (false);
713 
714  typename SearchGridT<typename ImageMesh<Basis>::Elem::index_type >::iterator it, eit;
715  if (elem_grid_->lookup(it,eit, p))
716  {
717  while (it != eit)
718  {
719  if (inside(typename ImageMesh<Basis>::Elem::index_type(*it), p))
720  {
721  elem = static_cast<INDEX>(*it);
722  return (true);
723  }
724  ++it;
725  }
726  }
727  return (false);
728  }
729 
730 
731  template <class INDEX, class ARRAY>
732  inline bool locate_elem(INDEX &elem, ARRAY& coords, const Core::Geometry::Point &p) const
733  {
734  if (this->basis_.polynomial_order() > 1) return elem_locate(elem, *this, p);
735 
736  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
737  "StructQuadSurfMesh::locate_node requires synchronize(ELEM_LOCATE_E).")
738 
739  typename ImageMesh<Basis>::Elem::size_type sz;
740  this->size(sz);
741  if ((elem > 0)&&(elem < sz))
742  {
743  elem.mesh_ = this;
744  if (inside(elem,p))
745  {
746  ElemData ed(*this,elem);
747  this->basis_.get_coords(coords,p,ed);
748  return (true);
749  }
750  }
751 
752  if (!elem_grid_) return (false);
753 
754  typename SearchGridT<typename ImageMesh<Basis>::Elem::index_type >::iterator it, eit;
755  if (elem_grid_->lookup(it,eit, p))
756  {
757  while (it != eit)
758  {
759  if (inside(typename ImageMesh<Basis>::Elem::index_type(*it), p))
760  {
761  elem = static_cast<INDEX>(*it);
762  ElemData ed(*this,elem);
763  this->basis_.get_coords(coords,p,ed);
764  return (true);
765  }
766  ++it;
767  }
768  }
769  return (false);
770  }
771 
772  template <class ARRAY>
773  bool find_closest_elem(double& pdist, Core::Geometry::Point &result, ARRAY& coords,
774  typename ImageMesh<Basis>::Elem::index_type &elem, const Core::Geometry::Point &p) const
775  {
776  return (find_closest_elem(pdist,result,coords,elem,p,-1.0));
777  }
778 
779  template <class ARRAY>
780  bool find_closest_elem(double& pdist, Core::Geometry::Point &result, ARRAY& coords,
781  typename ImageMesh<Basis>::Elem::index_type &elem,
782  const Core::Geometry::Point &p, double maxdist) const
783  {
784  /// If there are no nodes we cannot find the closest one
785  if (this->ni_ < 2 || this->nj_ < 2) return (false);
786 
787  if (maxdist < 0.0) maxdist = DBL_MAX; else maxdist = maxdist*maxdist;
788 
789  /// Test the one in face that is an initial guess
790  if (elem.i_ >= 0 && elem.i_ < (this->ni_-1) &&
791  elem.j_ >= 0 && elem.j_ < (this->nj_-1))
792  {
793  elem.mesh_ = this;
794  typename ImageMesh<Basis>::Node::array_type nodes;
795 
796  this->get_nodes(nodes,elem);
797  est_closest_point_on_quad(result, p,
798  points_(nodes[0].j_,nodes[0].i_),
799  points_(nodes[1].j_,nodes[1].i_),
800  points_(nodes[2].j_,nodes[2].i_),
801  points_(nodes[3].j_,nodes[3].i_));
802 
803  double dist = (p-result).length2();
804  if ( dist < epsilon2_ )
805  {
806  /// As we computed an estimate, we use the Newton's method in the basis functions
807  /// compute a more refined solution. This function may slow down computation.
808  /// This piece of code will calculate the coordinates in the local element framework
809  /// (the newton's method of finding a minimum), then it will project this back
810  /// THIS CODE SHOULD BE FURTHER OPTIMIZED
811 
812  ElemData ed(*this,elem);
813  this->basis_.get_coords(coords,result,ed);
814 
815  result = this->basis_.interpolate(coords,ed);
816  pdist = sqrt((result-p).length2());
817  return (true);
818  }
819  }
820 
821  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
822  "StructQuadSurfMesh::find_closest_elem requires synchronize(ELEM_LOCATE_E).")
823 
824  // get grid sizes
825  const size_type ni = elem_grid_->get_ni()-1;
826  const size_type nj = elem_grid_->get_nj()-1;
827  const size_type nk = elem_grid_->get_nk()-1;
828 
829  // Convert to grid coordinates.
830  index_type bi, bj, bk, ei, ej, ek;
831  elem_grid_->unsafe_locate(bi, bj, bk, p);
832 
833  // Clamp to closest point on the grid.
834  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
835  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
836  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
837 
838  ei = bi; ej = bj; ek = bk;
839 
840  double dmin = maxdist;
841  bool found;
842  bool found_one = false;
843 
844  typename ImageMesh<Basis>::Node::array_type nodes;
845  typename ImageMesh<Basis>::Elem::index_type idx;
846 
847  do
848  {
849  found = true;
850  /// This looks incorrect - but it is correct
851  /// We need to do a full shell without any elements that are closer
852  /// to make sure there no closer elements in neighboring searchgrid cells
853 
854  for (index_type i = bi; i <= ei; i++)
855  {
856  if (i < 0 || i > ni) continue;
857  for (index_type j = bj; j <= ej; j++)
858  {
859  if (j < 0 || j > nj) continue;
860  for (index_type k = bk; k <= ek; k++)
861  {
862  if (k < 0 || k > nk) continue;
863  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
864  {
865  if (elem_grid_->min_distance_squared(p, i, j, k) < dmin)
866  {
867  found = false;
868  typename SearchGridT<typename ImageMesh<Basis>::Elem::index_type>::iterator it, eit;
869  elem_grid_->lookup_ijk(it, eit, i, j, k);
870 
871  while (it != eit)
872  {
873  Core::Geometry::Point r;
874  this->get_nodes(nodes,*it);
875  est_closest_point_on_quad(r, p,
876  points_(nodes[0].j_,nodes[0].i_),
877  points_(nodes[1].j_,nodes[1].i_),
878  points_(nodes[2].j_,nodes[2].i_),
879  points_(nodes[3].j_,nodes[3].i_));
880  const double dtmp = (p - r).length2();
881  if (dtmp < dmin)
882  {
883  found_one = true;
884  result = r;
885  idx = *it;
886  dmin = dtmp;
887 
888  if (dmin < epsilon2_)
889  {
890  /// As we computed an estimate, we use the Newton's method in the basis functions
891  /// compute a more refined solution. This function may slow down computation.
892  /// This piece of code will calculate the coordinates in the local element framework
893  /// (the newton's method of finding a minimum), then it will project this back
894  /// THIS CODE SHOULD BE FURTHER OPTIMIZED
895 
896  ElemData ed(*this,idx);
897  this->basis_.get_coords(coords,result,ed);
898 
899  result = this->basis_.interpolate(coords,ed);
900  dmin = (result-p).length2();
901 
902  elem = idx;
903  pdist = sqrt(dmin);
904  return (true);
905  }
906  }
907  ++it;
908  }
909  }
910  }
911  }
912  }
913  }
914  bi--;ei++;
915  bj--;ej++;
916  bk--;ek++;
917  }
918  while (!found);
919 
920  if (!found_one) return (false);
921 
922  /// As we computed an estimate, we use the Newton's method in the basis functions
923  /// compute a more refined solution. This function may slow down computation.
924  /// This piece of code will calculate the coordinates in the local element framework
925  /// (the newton's method of finding a minimum), then it will project this back
926  /// THIS CODE SHOULD BE FURTHER OPTIMIZED
927 
928  ElemData ed(*this,idx);
929  this->basis_.get_coords(coords,result,ed);
930 
931  result = this->basis_.interpolate(coords,ed);
932  dmin = (result-p).length2();
933 
934  elem = idx;
935  pdist = dmin;
936  return (true);
937  }
938 
939 
940 
941  /// Export this class using the old Pio system
942  virtual void io(Piostream&);
943  /// These IDs are created as soon as this class will be instantiated
944  /// The first one is for Pio and the second for the virtual interface
945  /// These are currently different as they serve different needs. static PersistentTypeID type_idts;
946  static PersistentTypeID type_idsqs;
947  /// Core functionality for getting the name of a templated mesh class
948  static const std::string type_name(int n = -1);
949 
950  virtual bool synchronize(mask_type sync);
951  virtual bool unsynchronize(mask_type sync);
952  bool clear_synchronization();
953 
954  /// Type description, used for finding names of the mesh class for
955  /// dynamic compilation purposes. Soem of this should be obsolete
956  virtual const TypeDescription *get_type_description() const;
957  static const TypeDescription* node_type_description();
958  static const TypeDescription* edge_type_description();
959  static const TypeDescription* face_type_description();
960  static const TypeDescription* cell_type_description();
961  static const TypeDescription* elem_type_description()
962  { return face_type_description(); }
963 
964  /// This function returns a maker for Pio.
965  static Persistent *maker() { return new StructQuadSurfMesh<Basis>(); }
966  /// This function returns a handle for the virtual interface.
967  static MeshHandle mesh_maker() { return boost::make_shared<StructQuadSurfMesh<Basis>>(); }
968  /// This function returns a handle for the virtual interface.
969  static MeshHandle structquadsurf_maker(size_type x ,size_type y) { return boost::make_shared<StructQuadSurfMesh<Basis>>(x,y); }
970 
971  Array2<Core::Geometry::Point>& get_points() { return (points_); }
972 
973 
974 
975  double get_epsilon() const
976  { return (epsilon_); }
977 
978 protected:
979 
980  void compute_node_grid(Core::Geometry::BBox& bb);
981  void compute_elem_grid(Core::Geometry::BBox& bb);
982  void compute_epsilon(Core::Geometry::BBox& bb);
983 
984  void compute_normals();
985 
986  /// Used to recompute data for individual cells.
987  void insert_elem_into_grid(typename ImageMesh<Basis>::Elem::index_type ci);
988  void remove_elem_from_grid(typename ImageMesh<Basis>::Elem::index_type ci);
989 
990  void insert_node_into_grid(typename ImageMesh<Basis>::Node::index_type ci);
991  void remove_node_from_grid(typename ImageMesh<Basis>::Node::index_type ci);
992 
993  /// compute_edge_neighbors is not defined anywhere... do don't
994  /// declare it... void compute_edge_neighbors(double err = 1.0e-8);
995 
996  const Core::Geometry::Point &point(const typename ImageMesh<Basis>::Node::index_type &idx)
997  { return points_(idx.j_, idx.i_); }
998 
999  index_type next(index_type i) { return ((i%4)==3) ? (i-3) : (i+1); }
1000  index_type prev(index_type i) { return ((i%4)==0) ? (i+3) : (i-1); }
1001 
1002  Array2<Core::Geometry::Point> points_;
1003  Array2<Core::Geometry::Vector> normals_; /// normalized per node
1004 
1005  boost::shared_ptr<SearchGridT<typename ImageMesh<Basis>::Node::index_type > > node_grid_;
1006  boost::shared_ptr<SearchGridT<typename ImageMesh<Basis>::Elem::index_type > > elem_grid_;
1007 
1008  mutable Core::Thread::Mutex synchronize_lock_;
1009  mask_type synchronized_;
1010  double epsilon_;
1011  double epsilon2_;
1012 
1013 
1014 };
1015 
1016 
1017 template <class Basis>
1018 PersistentTypeID
1019 StructQuadSurfMesh<Basis>::type_idsqs(StructQuadSurfMesh<Basis>::type_name(-1),
1020  "Mesh", maker);
1021 
1022 
1023 template <class Basis>
1024 StructQuadSurfMesh<Basis>::StructQuadSurfMesh()
1025  : synchronize_lock_("StructQuadSurfMesh Normals Lock"),
1026  synchronized_(Mesh::ALL_ELEMENTS_E),
1027  epsilon_(0.0),
1028  epsilon2_(0.0)
1029 {
1030  DEBUG_CONSTRUCTOR("StructQuadSurfMesh")
1031 
1032  /// Create a new virtual interface for this copy ! all pointers have
1033  /// changed hence create a new ! virtual interface class
1034  this->vmesh_.reset(CreateVStructQuadSurfMesh(this));
1035 }
1036 
1037 
1038 template <class Basis>
1039 StructQuadSurfMesh<Basis>::StructQuadSurfMesh(size_type x, size_type y)
1040  : ImageMesh<Basis>(x, y, Core::Geometry::Point(0.0, 0.0, 0.0), Core::Geometry::Point(1.0, 1.0, 1.0)),
1041  points_( y,x),
1042  normals_( y,x),
1043  synchronize_lock_("StructQuadSurfMesh Normals Lock"),
1044  synchronized_(Mesh::ALL_ELEMENTS_E),
1045  epsilon_(0.0),
1046  epsilon2_(0.0)
1047 {
1048  DEBUG_CONSTRUCTOR("StructQuadSurfMesh")
1049 
1050  /// Create a new virtual interface for this copy ! all pointers have
1051  /// changed hence create a new ! virtual interface class
1052  this->vmesh_.reset(CreateVStructQuadSurfMesh(this));
1053 }
1054 
1055 
1056 template <class Basis>
1057 StructQuadSurfMesh<Basis>::StructQuadSurfMesh(const StructQuadSurfMesh &copy)
1058  : ImageMesh<Basis>(copy),
1059  synchronize_lock_("StructQuadSurfMesh Normals Lock"),
1060  synchronized_(Mesh::NODES_E | Mesh::FACES_E | Mesh::CELLS_E)
1061 {
1062  DEBUG_CONSTRUCTOR("StructQuadSurfMesh")
1063 
1064  copy.synchronize_lock_.lock();
1065 
1066  points_ = copy.points_;
1067 
1068  epsilon_ = copy.epsilon_;
1069  epsilon2_ = copy.epsilon2_;
1070 
1071  synchronized_ |= copy.synchronized_ & Mesh::EPSILON_E;
1072 
1073  copy.synchronize_lock_.unlock();
1074 
1075  /// Create a new virtual interface for this copy ! all pointers have
1076  /// changed hence create a new ! virtual interface class
1077  this->vmesh_.reset(CreateVStructQuadSurfMesh(this));
1078 }
1079 
1080 
1081 template <class Basis>
1082 bool
1083 StructQuadSurfMesh<Basis>::get_dim(std::vector<size_type> &array) const
1084 {
1085  array.resize(2);
1086  array.clear();
1087 
1088  array.push_back(this->ni_);
1089  array.push_back(this->nj_);
1090 
1091  return true;
1092 }
1093 
1094 
1095 template <class Basis>
1096 Core::Geometry::BBox
1097 StructQuadSurfMesh<Basis>::get_bounding_box() const
1098 {
1099  Core::Geometry::BBox result;
1100 
1101  typename ImageMesh<Basis>::Node::iterator ni, nie;
1102  this->begin(ni);
1103  this->end(nie);
1104  while (ni != nie)
1105  {
1106  Core::Geometry::Point p;
1107  get_center(p, *ni);
1108  result.extend(p);
1109  ++ni;
1110  }
1111  return result;
1112 }
1113 
1114 
1115 template <class Basis>
1116 void
1117 StructQuadSurfMesh<Basis>::transform(const Core::Geometry::Transform &t)
1118 {
1119  synchronize_lock_.lock();
1120 
1121  typename ImageMesh<Basis>::Node::iterator i, ie;
1122  this->begin(i);
1123  this->end(ie);
1124 
1125  while (i != ie)
1126  {
1127  points_((*i).j_,(*i).i_) = t.project(points_((*i).j_,(*i).i_));
1128  ++i;
1129  }
1130 
1131  if (node_grid_) { node_grid_->transform(t); }
1132  if (elem_grid_) { elem_grid_->transform(t); }
1133  synchronize_lock_.unlock();
1134 }
1135 
1136 
1137 template <class Basis>
1138 void
1139 StructQuadSurfMesh<Basis>::get_normal(Core::Geometry::Vector &result,
1140  const typename ImageMesh<Basis>::Node::index_type &idx ) const
1141 {
1142  result = normals_(idx.j_, idx.i_);
1143 }
1144 
1145 
1146 template <class Basis>
1147 void
1148 StructQuadSurfMesh<Basis>::get_center(Core::Geometry::Point &result,
1149  const typename ImageMesh<Basis>::Node::index_type &idx) const
1150 {
1151  result = points_(idx.j_, idx.i_);
1152 }
1153 
1154 
1155 template <class Basis>
1156 void
1157 StructQuadSurfMesh<Basis>::get_center(Core::Geometry::Point &result,
1158  typename ImageMesh<Basis>::Edge::index_type idx) const
1159 {
1160  typename ImageMesh<Basis>::Node::array_type arr;
1161  this->get_nodes(arr, idx);
1162  Core::Geometry::Point p1;
1163  get_center(result, arr[0]);
1164  get_center(p1, arr[1]);
1165 
1166  result += p1;
1167  result *= 0.5;
1168 }
1169 
1170 
1171 template <class Basis>
1172 void
1173 StructQuadSurfMesh<Basis>::get_center(Core::Geometry::Point &p,
1174  const typename ImageMesh<Basis>::Face::index_type &idx) const
1175 {
1176  typename ImageMesh<Basis>::Node::array_type nodes;
1177  this->get_nodes(nodes, idx);
1178  ASSERT(nodes.size() == 4);
1179  typename ImageMesh<Basis>::Node::array_type::iterator nai = nodes.begin();
1180  get_point(p, *nai);
1181  ++nai;
1182  Core::Geometry::Point pp;
1183  while (nai != nodes.end())
1184  {
1185  get_point(pp, *nai);
1186  p += pp;
1187  ++nai;
1188  }
1189  p *= (1.0 / 4.0);
1190 }
1191 
1192 template <class Basis>
1193 int
1194 StructQuadSurfMesh<Basis>::get_weights(const Core::Geometry::Point &p,
1195  typename ImageMesh<Basis>::Face::array_type &l,
1196  double *w)
1197 {
1198  typename ImageMesh<Basis>::Face::index_type idx;
1199  if (locate_elem(idx, p))
1200  {
1201  l.resize(1);
1202  l[0] = idx;
1203  w[0] = 1.0;
1204  return 1;
1205  }
1206  return 0;
1207 }
1208 
1209 
1210 template <class Basis>
1211 int
1212 StructQuadSurfMesh<Basis>::get_weights(const Core::Geometry::Point &p,
1213  typename ImageMesh<Basis>::Node::array_type &l,
1214  double *w)
1215 {
1216  typename ImageMesh<Basis>::Face::index_type idx;
1217  if (locate_elem(idx, p))
1218  {
1219  this->get_nodes(l,idx);
1220  std::vector<double> coords(2);
1221  if (get_coords(coords, p, idx)) {
1222  this->basis_.get_weights(coords, w);
1223  return this->basis_.dofs();
1224  }
1225  }
1226  return 0;
1227 }
1228 
1229 
1230 template <class Basis>
1231 bool
1232 StructQuadSurfMesh<Basis>::inside3_p(typename ImageMesh<Basis>::Face::index_type i,
1233  const Core::Geometry::Point &p) const
1234 {
1235  typename ImageMesh<Basis>::Node::array_type nodes;
1236  this->get_nodes(nodes, i);
1237 
1238  size_type n = nodes.size();
1239 
1240  std::vector<Core::Geometry::Point> pts(n);
1241 
1242  for (index_type i = 0; i < n; i++)
1243  {
1244  get_center(pts[i], nodes[i]);
1245  }
1246 
1247  for (index_type i = 0; i < n; i+=2)
1248  {
1249  Core::Geometry::Point p0 = pts[(i+0)%n];
1250  Core::Geometry::Point p1 = pts[(i+1)%n];
1251  Core::Geometry::Point p2 = pts[(i+2)%n];
1252 
1253  Core::Geometry::Vector v01(p0-p1);
1254  Core::Geometry::Vector v02(p0-p2);
1255  Core::Geometry::Vector v0(p0-p);
1256  Core::Geometry::Vector v1(p1-p);
1257  Core::Geometry::Vector v2(p2-p);
1258  const double a = Cross(v01, v02).length(); /// area of the whole triangle (2x)
1259  const double a0 = Cross(v1, v2).length(); /// area opposite p0
1260  const double a1 = Cross(v2, v0).length(); /// area opposite p1
1261  const double a2 = Cross(v0, v1).length(); /// area opposite p2
1262  const double s = a0+a1+a2;
1263 
1264  /// If the area of any of the sub triangles is very small then the point
1265  /// is on the edge of the subtriangle.
1266  /// @todo : How small is small ???
1267 /// if( a0 < MIN_ELEMENT_VAL ||
1268 /// a1 < MIN_ELEMENT_VAL ||
1269 /// a2 < MIN_ELEMENT_VAL )
1270 /// return true;
1271 
1272  /// For the point to be inside a CONVEX quad it must be inside one
1273  /// of the four triangles that can be formed by using three of the
1274  /// quad vertices and the point in question.
1275  if( std::fabs(s - a) < epsilon2_*epsilon2_ && a > epsilon2_*epsilon2_ )
1276  {
1277  return true;
1278  }
1279  }
1280  return false;
1281 }
1282 
1283 
1284 template <class Basis>
1285 void
1286 StructQuadSurfMesh<Basis>::get_random_point(Core::Geometry::Point &p,
1287  const typename ImageMesh<Basis>::Elem::index_type &ei,
1288  FieldRNG &rng) const
1289 {
1290  const Core::Geometry::Point &a0 = points_(ei.j_ + 0, ei.i_ + 0);
1291  const Core::Geometry::Point &a1 = points_(ei.j_ + 0, ei.i_ + 1);
1292  const Core::Geometry::Point &a2 = points_(ei.j_ + 1, ei.i_ + 1);
1293  const Core::Geometry::Point &a3 = points_(ei.j_ + 1, ei.i_ + 0);
1294 
1295  const double aarea = Cross(a1 - a0, a2 - a0).length();
1296  const double barea = Cross(a3 - a0, a2 - a0).length();
1297 
1298  /// Fold the quad sample into a triangle.
1299  double u = rng();
1300  double v = rng();
1301  if (u + v > 1.0) { u = 1.0 - u; v = 1.0 - v; }
1302 
1303  if (rng() * (aarea + barea) < aarea)
1304  {
1305  p = a0+((a1-a0)*u)+((a2-a0)*v);
1306  }
1307  else
1308  {
1309  p = a0+((a3-a0)*u)+((a2-a0)*v);
1310  }
1311 }
1312 
1313 
1314 template <class Basis>
1315 void
1316 StructQuadSurfMesh<Basis>::set_point(const Core::Geometry::Point &point,
1317  const typename ImageMesh<Basis>::Node::index_type &index)
1318 {
1319  points_(index.j_, index.i_) = point;
1320 }
1321 
1322 
1323 template <class Basis>
1324 bool
1325 StructQuadSurfMesh<Basis>::synchronize(mask_type sync)
1326 {
1327  synchronize_lock_.lock();
1328  if (sync & Mesh::NORMALS_E && !(synchronized_ & Mesh::NORMALS_E))
1329  {
1330  compute_normals();
1331  synchronized_ |= Mesh::NORMALS_E;
1332  }
1333 
1334  if (sync & (Mesh::NODE_LOCATE_E|Mesh::ELEM_LOCATE_E|Mesh::EPSILON_E) &&
1335  ( !(synchronized_ & Mesh::NODE_LOCATE_E) ||
1336  !(synchronized_ & Mesh::ELEM_LOCATE_E) ||
1337  !(synchronized_ & Mesh::EPSILON_E) ))
1338  {
1339  /// These computations share the evalution of the bounding box
1340  Core::Geometry::BBox bb = get_bounding_box();
1341 
1342  /// Compute the epsilon for geometrical closeness comparisons
1343  /// Mainly used by the grid lookup tables
1344  if (sync & (Mesh::EPSILON_E|Mesh::LOCATE_E|Mesh::FIND_CLOSEST_E) &&
1345  !(synchronized_ & Mesh::EPSILON_E))
1346  {
1347  compute_epsilon(bb);
1348  }
1349 
1350  /// Table for finding nodes in space
1351  if (sync & (Mesh::NODE_LOCATE_E|Mesh::FIND_CLOSEST_NODE_E) &&
1352  !(synchronized_ & Mesh::NODE_LOCATE_E))
1353  {
1354  compute_node_grid(bb);
1355  }
1356 
1357  /// Table for finding elements in space
1358  if (sync & (Mesh::ELEM_LOCATE_E|Mesh::FIND_CLOSEST_ELEM_E) &&
1359  !(synchronized_ & Mesh::ELEM_LOCATE_E))
1360  {
1361  compute_elem_grid(bb);
1362  }
1363  }
1364 
1365  synchronize_lock_.unlock();
1366  return(true);
1367 }
1368 
1369 template <class Basis>
1370 bool
1371 StructQuadSurfMesh<Basis>::unsynchronize(mask_type /*sync*/)
1372 {
1373  return (true);
1374 }
1375 
1376 
1377 template <class Basis>
1378 bool
1379 StructQuadSurfMesh<Basis>::clear_synchronization()
1380 {
1381  synchronize_lock_.lock();
1382  // Undo marking the synchronization
1383  synchronized_ = Mesh::NODES_E | Mesh::ELEMS_E | Mesh::FACES_E;
1384 
1385  // Free memory where possible
1386  node_grid_.reset();
1387  elem_grid_.reset();
1388 
1389  synchronize_lock_.unlock();
1390  return (true);
1391 }
1392 
1393 template <class Basis>
1394 void
1395 StructQuadSurfMesh<Basis>::insert_elem_into_grid(typename ImageMesh<Basis>::Elem::index_type idx)
1396 {
1397  /// @todo: This can crash if you insert a new cell outside of the grid.
1398  // Need to recompute grid at that point.
1399 
1400  Core::Geometry::BBox box;
1401  box.extend(points_(idx.j_,idx.i_));
1402  box.extend(points_(idx.j_+1,idx.i_));
1403  box.extend(points_(idx.j_,idx.i_+1));
1404  box.extend(points_(idx.j_+1,idx.i_+1));
1405  box.extend(epsilon_);
1406  elem_grid_->insert(idx, box);
1407 }
1408 
1409 
1410 template <class Basis>
1411 void
1412 StructQuadSurfMesh<Basis>::remove_elem_from_grid(typename ImageMesh<Basis>::Elem::index_type idx)
1413 {
1414  Core::Geometry::BBox box;
1415  box.extend(points_(idx.j_,idx.i_));
1416  box.extend(points_(idx.j_+1,idx.i_));
1417  box.extend(points_(idx.j_,idx.i_+1));
1418  box.extend(points_(idx.j_+1,idx.i_+1));
1419  box.extend(epsilon_);
1420  elem_grid_->remove(idx, box);
1421 }
1422 
1423 
1424 template <class Basis>
1425 void
1426 StructQuadSurfMesh<Basis>::insert_node_into_grid(typename ImageMesh<Basis>::Node::index_type ni)
1427 {
1428  /// @todo: This can crash if you insert a new cell outside of the grid.
1429  // Need to recompute grid at that point.
1430  node_grid_->insert(ni,points_[ni]);
1431 }
1432 
1433 
1434 template <class Basis>
1435 void
1436 StructQuadSurfMesh<Basis>::remove_node_from_grid(typename ImageMesh<Basis>::Node::index_type ni)
1437 {
1438  node_grid_->remove(ni,points_[ni]);
1439 }
1440 
1441 template <class Basis>
1442 void
1443 StructQuadSurfMesh<Basis>::compute_node_grid(Core::Geometry::BBox& bb)
1444 {
1445  if (bb.valid())
1446  {
1447  // Cubed root of number of elems to get a subdivision ballpark.
1448 
1449  typename ImageMesh<Basis>::Node::size_type esz;
1450  this->size(esz);
1451 
1452  const size_type s = 3*static_cast<size_type>
1453  ((ceil(pow(static_cast<double>(esz) , (1.0/3.0))))/2.0 + 1.0);
1454 
1455  Core::Geometry::Vector diag = bb.diagonal();
1456  double trace = (diag.x()+diag.y()+diag.z());
1457  size_type sx = static_cast<size_type>(ceil(diag.x()/trace*s));
1458  size_type sy = static_cast<size_type>(ceil(diag.y()/trace*s));
1459  size_type sz = static_cast<size_type>(ceil(diag.z()/trace*s));
1460 
1461  Core::Geometry::BBox b = bb; b.extend(10*epsilon_);
1462  node_grid_.reset(new SearchGridT<typename ImageMesh<Basis>::Node::index_type >(sx, sy, sz, b.min(), b.max()));
1463 
1464  typename ImageMesh<Basis>::Node::iterator ni, nie;
1465  this->begin(ni);
1466  this->end(nie);
1467  while(ni != nie)
1468  {
1469  insert_node_into_grid(*ni);
1470  ++ni;
1471  }
1472  }
1473 
1474  synchronized_ |= Mesh::NODE_LOCATE_E;
1475 }
1476 
1477 template <class Basis>
1478 void
1479 StructQuadSurfMesh<Basis>::compute_elem_grid(Core::Geometry::BBox& bb)
1480 {
1481  if (bb.valid())
1482  {
1483  // Cubed root of number of elems to get a subdivision ballpark.
1484 
1485  typename ImageMesh<Basis>::Elem::size_type esz;
1486  this->size(esz);
1487 
1488  const size_type s = 3*static_cast<size_type>
1489  ((ceil(pow(static_cast<double>(esz) , (1.0/3.0))))/2.0 + 1.0);
1490 
1491  Core::Geometry::Vector diag = bb.diagonal();
1492  double trace = (diag.x()+diag.y()+diag.z());
1493  size_type sx = static_cast<size_type>(ceil(diag.x()/trace*s));
1494  size_type sy = static_cast<size_type>(ceil(diag.y()/trace*s));
1495  size_type sz = static_cast<size_type>(ceil(diag.z()/trace*s));
1496 
1497  Core::Geometry::BBox b = bb; b.extend(10*epsilon_);
1498  elem_grid_.reset(new SearchGridT<typename ImageMesh<Basis>::Elem::index_type>(sx, sy, sz, b.min(), b.max()));
1499 
1500  typename ImageMesh<Basis>::Elem::iterator ci, cie;
1501  this->begin(ci);
1502  this->end(cie);
1503  while(ci != cie)
1504  {
1505  insert_elem_into_grid(*ci);
1506  ++ci;
1507  }
1508  }
1509 
1510  synchronized_ |= Mesh::ELEM_LOCATE_E;
1511 }
1512 
1513 template <class Basis>
1514 void
1515 StructQuadSurfMesh<Basis>::compute_epsilon(Core::Geometry::BBox& bb)
1516 {
1517  epsilon_ = bb.diagonal().length();
1518  epsilon2_ = epsilon_*epsilon_;
1519 }
1520 
1521 
1522 template <class Basis>
1523 void
1524 StructQuadSurfMesh<Basis>::compute_normals()
1525 {
1526  normals_.resize(points_.dim1(), points_.dim2()); /// 1 per node
1527 
1528  /// build table of faces that touch each node
1529  Array2< std::vector<typename ImageMesh<Basis>::Face::index_type> >
1530  node_in_faces(points_.dim1(), points_.dim2());
1531 
1532  /// face normals (not normalized) so that magnitude is also the area.
1533  Array2<Core::Geometry::Vector> face_normals((points_.dim1()-1),(points_.dim2()-1));
1534 
1535  /// Computing normal per face.
1536  typename ImageMesh<Basis>::Node::array_type nodes(4);
1537  typename ImageMesh<Basis>::Face::iterator iter, iter_end;
1538  this->begin(iter);
1539  this->end(iter_end);
1540  while (iter != iter_end)
1541  {
1542  this->get_nodes(nodes, *iter);
1543 
1544  Core::Geometry::Point p0, p1, p2, p3;
1545  get_point(p0, nodes[0]);
1546  get_point(p1, nodes[1]);
1547  get_point(p2, nodes[2]);
1548  get_point(p3, nodes[3]);
1549 
1550  /// build table of faces that touch each node
1551  node_in_faces(nodes[0].j_,nodes[0].i_).push_back(*iter);
1552  node_in_faces(nodes[1].j_,nodes[1].i_).push_back(*iter);
1553  node_in_faces(nodes[2].j_,nodes[2].i_).push_back(*iter);
1554  node_in_faces(nodes[3].j_,nodes[3].i_).push_back(*iter);
1555 
1556  Core::Geometry::Vector v0 = p1 - p0;
1557  Core::Geometry::Vector v1 = p2 - p1;
1558  Core::Geometry::Vector n = Cross(v0, v1);
1559  face_normals((*iter).j_, (*iter).i_) = n;
1560 
1561  ++iter;
1562  }
1563 
1564  /// Averaging the normals.
1565  typename ImageMesh<Basis>::Node::iterator nif_iter, nif_iter_end;
1566  this->begin( nif_iter );
1567  this->end( nif_iter_end );
1568 
1569  while (nif_iter != nif_iter_end)
1570  {
1571  std::vector<typename ImageMesh<Basis>::Face::index_type> v =
1572  node_in_faces((*nif_iter).j_, (*nif_iter).i_);
1573  typename std::vector<typename ImageMesh<Basis>::Face::index_type>::const_iterator
1574  fiter = v.begin();
1575  Core::Geometry::Vector ave(0.L,0.L,0.L);
1576  while(fiter != v.end())
1577  {
1578  ave += face_normals((*fiter).j_,(*fiter).i_);
1579  ++fiter;
1580  }
1581  ave.safe_normalize();
1582  normals_((*nif_iter).j_, (*nif_iter).i_) = ave;
1583  ++nif_iter;
1584  }
1585 }
1586 
1587 
1588 #define STRUCT_QUAD_SURF_MESH_VERSION 3
1589 
1590 template <class Basis>
1591 void
1592 StructQuadSurfMesh<Basis>::io(Piostream& stream)
1593 {
1594  int version =
1595  stream.begin_class(type_name(-1), STRUCT_QUAD_SURF_MESH_VERSION);
1596  ImageMesh<Basis>::io(stream);
1597 
1598  if (version ==2)
1599  {
1600  /// The dimensions of this array were swapped
1601  Array2<Core::Geometry::Point> tpoints;
1602  Pio(stream, tpoints);
1603 
1604  size_t dim1 = tpoints.dim1();
1605  size_t dim2 = tpoints.dim2();
1606 
1607  points_.resize(dim2,dim1);
1608  for (size_t i=0; i<dim1; i++)
1609  for (size_t j=0; j<dim2; j++)
1610  points_(j,i) = tpoints(i,j);
1611  }
1612  else
1613  {
1614  Pio(stream, points_);
1615  }
1616 
1617  stream.end_class();
1618 
1619  if (stream.reading())
1620  this->vmesh_.reset(CreateVStructQuadSurfMesh(this));
1621 }
1622 
1623 
1624 template <class Basis>
1625 const std::string
1626 StructQuadSurfMesh<Basis>::type_name(int n)
1627 {
1628  ASSERT((n >= -1) && n <= 1);
1629  if (n == -1)
1630  {
1631  static const std::string name = TypeNameGenerator::make_template_id(type_name(0), type_name(1));
1632  return name;
1633  }
1634  else if (n == 0)
1635  {
1636  static const std::string nm("StructQuadSurfMesh");
1637  return nm;
1638  }
1639  else
1640  {
1641  return find_type_name((Basis *)0);
1642  }
1643 }
1644 
1645 
1646 template <class Basis>
1647 const TypeDescription*
1648 get_type_description(StructQuadSurfMesh<Basis> *)
1649 {
1650  static TypeDescription *td = 0;
1651  if (!td)
1652  {
1653  const TypeDescription *sub = get_type_description((Basis*)0);
1654  TypeDescription::td_vec *subs = new TypeDescription::td_vec(1);
1655  (*subs)[0] = sub;
1656  td = new TypeDescription("StructQuadSurfMesh", subs,
1657  std::string(__FILE__),
1658  "SCIRun",
1659  TypeDescription::MESH_E);
1660  }
1661  return td;
1662 }
1663 
1664 
1665 template <class Basis>
1666 const TypeDescription*
1667 StructQuadSurfMesh<Basis>::get_type_description() const
1668 {
1669  return SCIRun::get_type_description((StructQuadSurfMesh<Basis> *)0);
1670 }
1671 
1672 
1673 template <class Basis>
1674 const TypeDescription*
1675 StructQuadSurfMesh<Basis>::node_type_description()
1676 {
1677  static TypeDescription *td = 0;
1678  if (!td)
1679  {
1680  const TypeDescription *me =
1681  SCIRun::get_type_description((StructQuadSurfMesh<Basis> *)0);
1682  td = new TypeDescription(me->get_name() + "::Node",
1683  std::string(__FILE__),
1684  "SCIRun",
1685  TypeDescription::MESH_E);
1686  }
1687  return td;
1688 }
1689 
1690 
1691 template <class Basis>
1692 const TypeDescription*
1693 StructQuadSurfMesh<Basis>::edge_type_description()
1694 {
1695  static TypeDescription *td = 0;
1696  if (!td)
1697  {
1698  const TypeDescription *me =
1699  SCIRun::get_type_description((StructQuadSurfMesh<Basis> *)0);
1700  td = new TypeDescription(me->get_name() + "::Edge",
1701  std::string(__FILE__),
1702  "SCIRun",
1703  TypeDescription::MESH_E);
1704  }
1705  return td;
1706 }
1707 
1708 
1709 template <class Basis>
1710 const TypeDescription*
1711 StructQuadSurfMesh<Basis>::face_type_description()
1712 {
1713  static TypeDescription *td = 0;
1714  if (!td)
1715  {
1716  const TypeDescription *me =
1717  SCIRun::get_type_description((StructQuadSurfMesh<Basis> *)0);
1718  td = new TypeDescription(me->get_name() + "::Face",
1719  std::string(__FILE__),
1720  "SCIRun",
1721  TypeDescription::MESH_E);
1722  }
1723  return td;
1724 }
1725 
1726 
1727 template <class Basis>
1728 const TypeDescription*
1729 StructQuadSurfMesh<Basis>::cell_type_description()
1730 {
1731  static TypeDescription *td = 0;
1732  if (!td)
1733  {
1734  const TypeDescription *me =
1735  SCIRun::get_type_description((StructQuadSurfMesh<Basis> *)0);
1736  td = new TypeDescription(me->get_name() + "::Cell",
1737  std::string(__FILE__),
1738  "SCIRun",
1739  TypeDescription::MESH_E);
1740  }
1741  return td;
1742 }
1743 
1744 
1745 template <class Basis>
1746 bool
1747 StructQuadSurfMesh<Basis>::find_closest_node(double& pdist, Core::Geometry::Point &result,
1748  typename ImageMesh<Basis>::Node::index_type &node, const Core::Geometry::Point &p,
1749  double maxdist) const
1750 {
1751  /// If there are no nodes we cannot find a closest point
1752  if (this->ni_ == 0 || this->nj_ == 0) return (false);
1753  if (maxdist < 0.0) maxdist = DBL_MAX; else maxdist = maxdist*maxdist;
1754 
1755  /// Check first guess
1756  if (node.i_ >= 0 && node.i_ < this->ni_ &&
1757  node.j_ >= 0 && node.j_ < this->nj_)
1758  {
1759  node.mesh_ = this;
1760  Core::Geometry::Point point = points_(node.j_,node.i_);
1761  double dist = (point-p).length2();
1762 
1763  if ( dist < epsilon2_ )
1764  {
1765  result = point;
1766  pdist = sqrt(dist);
1767  return (true);
1768  }
1769  }
1770 
1771  ASSERTMSG(synchronized_ & Mesh::NODE_LOCATE_E,
1772  "StructQuadSurfMesh::find_closest_elem requires synchronize(NODELOCATE_E).")
1773 
1774  // get grid sizes
1775  const size_type ni = node_grid_->get_ni()-1;
1776  const size_type nj = node_grid_->get_nj()-1;
1777  const size_type nk = node_grid_->get_nk()-1;
1778 
1779  // Convert to grid coordinates.
1780  index_type bi, bj, bk, ei, ej, ek;
1781  node_grid_->unsafe_locate(bi, bj, bk, p);
1782 
1783  // Clamp to closest point on the grid.
1784  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
1785  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
1786  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
1787 
1788  ei = bi; ej = bj; ek = bk;
1789 
1790  typename ImageMesh<Basis>::Node::array_type nodes;
1791 
1792  double dmin = maxdist;
1793  bool found = true;
1794  bool found_one = false;
1795 
1796  do
1797  {
1798  found = true;
1799  /// This looks incorrect - but it is correct
1800  /// We need to do a full shell without any elements that are closer
1801  /// to make sure there no closer elements in neighboring searchgrid cells
1802 
1803  for (index_type i = bi; i <= ei; i++)
1804  {
1805  if (i < 0 || i > ni) continue;
1806  for (index_type j = bj; j <= ej; j++)
1807  {
1808  if (j < 0 || j > nj) continue;
1809  for (index_type k = bk; k <= ek; k++)
1810  {
1811  if (k < 0 || k > nk) continue;
1812  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
1813  {
1814  if (node_grid_->min_distance_squared(p, i, j, k) < dmin)
1815  {
1816  typename SearchGridT<typename ImageMesh<Basis>::Node::index_type>::iterator it, eit;
1817  node_grid_->lookup_ijk(it,eit, i, j, k);
1818 
1819  while (it != eit)
1820  {
1821  const typename ImageMesh<Basis>::Node::index_type idx = *it;
1822  const Core::Geometry::Point pnt = points_(idx.j_,idx.i_);
1823  const double dist = (p-pnt).length2();
1824  if (dist < dmin)
1825  {
1826  found_one = true;
1827  result = pnt;
1828  node = idx;
1829  dmin = dist;
1830 
1831  /// If we are closer than eps^2 we found a node close enough
1832  if (dmin < epsilon2_)
1833  {
1834  pdist = sqrt(dmin);
1835  return (true);
1836  }
1837  }
1838  ++it;
1839  }
1840  }
1841  }
1842  }
1843  }
1844  }
1845  bi--;ei++;
1846  bj--;ej++;
1847  bk--;ek++;
1848  }
1849  while (!found) ;
1850 
1851  if (!found_one) return (false);
1852 
1853  pdist = sqrt(dmin);
1854  return (true);
1855 }
1856 
1857 
1858 /// This function will find the closest element and the location on that
1859 /// element that is the closest
1860 
1861 template <class Basis>
1862 bool
1863 StructQuadSurfMesh<Basis>::find_closest_elem(double& pdist,
1864  Core::Geometry::Point &result,
1865  typename ImageMesh<Basis>::Elem::index_type &elem,
1866  const Core::Geometry::Point &p) const
1867 {
1868  StackVector<double,2> coords;
1869  return(find_closest_elem(pdist,result,coords,elem,p,-1.0));
1870 }
1871 
1872 
1873 template <class Basis>
1874 bool
1875 StructQuadSurfMesh<Basis>::find_closest_elems(double& pdist,
1876  Core::Geometry::Point &result,
1877  std::vector<typename ImageMesh<Basis>::Elem::index_type> &elems,
1878  const Core::Geometry::Point &p) const
1879 {
1880  elems.clear();
1881 
1882  /// If there are no nodes we cannot find the closest one
1883  if (this->ni_ < 2 || this->nj_ < 2) return (false);
1884 
1885  ASSERTMSG(synchronized_ & Mesh::ELEM_LOCATE_E,
1886  "StructQuadSurfMesh::find_closest_elem requires synchronize(ELEM_LOCATE_E).")
1887 
1888  // get grid sizes
1889  const size_type ni = elem_grid_->get_ni()-1;
1890  const size_type nj = elem_grid_->get_nj()-1;
1891  const size_type nk = elem_grid_->get_nk()-1;
1892 
1893  // Convert to grid coordinates.
1894  index_type bi, bj, bk, ei, ej, ek;
1895  elem_grid_->unsafe_locate(bi, bj, bk, p);
1896 
1897  // Clamp to closest point on the grid.
1898  if (bi > ni) bi = ni; if (bi < 0) bi = 0;
1899  if (bj > nj) bj = nj; if (bj < 0) bj = 0;
1900  if (bk > nk) bk = nk; if (bk < 0) bk = 0;
1901 
1902  ei = bi; ej = bj; ek = bk;
1903 
1904  double dmin = DBL_MAX;
1905 
1906  typename ImageMesh<Basis>::Node::array_type nodes;
1907  typename ImageMesh<Basis>::Elem::index_type last_idx;
1908 
1909  bool found = true;
1910  do
1911  {
1912  found = true;
1913  /// This looks incorrect - but it is correct
1914  /// We need to do a full shell without any elements that are closer
1915  /// to make sure there no closer elements
1916  for (index_type i = bi; i <= ei; i++)
1917  {
1918  if (i < 0 || i > ni) continue;
1919  for (index_type j = bj; j <= ej; j++)
1920  {
1921  if (j < 0 || j > nj) continue;
1922  for (index_type k = bk; k <= ek; k++)
1923  {
1924  if (k < 0 || k > nk) continue;
1925  if (i == bi || i == ei || j == bj || j == ej || k == bk || k == ek)
1926  {
1927  if (elem_grid_->min_distance_squared(p, i, j, k) < dmin)
1928  {
1929  found = false;
1930  typename SearchGridT<typename ImageMesh<Basis>::Elem::index_type>::iterator it, eit;
1931  elem_grid_->lookup_ijk(it,eit, i, j, k);
1932 
1933  while (it != eit)
1934  {
1935  Core::Geometry::Point r;
1936  this->get_nodes(nodes,*it);
1937  est_closest_point_on_quad(r, p,
1938  points_(nodes[0].j_,nodes[0].i_),
1939  points_(nodes[1].j_,nodes[1].i_),
1940  points_(nodes[2].j_,nodes[2].i_),
1941  points_(nodes[3].j_,nodes[3].i_));
1942  const double dtmp = (p - r).length2();
1943 
1944  if (dtmp < dmin - epsilon2_)
1945  {
1946  elems.clear();
1947  result = r;
1948  last_idx = *it;
1949  elems.push_back(*it);
1950  dmin = dtmp;
1951  }
1952  else if (dtmp < dmin + epsilon2_)
1953  {
1954  elems.push_back(*it);
1955  }
1956  ++it;
1957  }
1958  }
1959  }
1960  }
1961  }
1962  }
1963  bi--;ei++;
1964  bj--;ej++;
1965  bk--;ek++;
1966  }
1967  while ((!found)||(dmin == DBL_MAX)) ;
1968 
1969 
1970  // As we computed an estimate, we use the Newton's method in the basis functions
1971  // compute a more refined solution. This function may slow down computation.
1972  // This piece of code will calculate the coordinates in the local element framework
1973  // (the newton's method of finding a minimum), then it will project this back
1974  // THIS CODE SHOULD BE FURTHER OPTIMIZED
1975 
1976  if (elems.size() == 1)
1977  {
1978  // if the number of faces is more then one the point we found is located
1979  // on the node or on the edge, which means the estimate is correct.
1980  StackVector<double,3> coords;
1981  ElemData ed(*this,last_idx);
1982  this->basis_.get_coords(coords,result,ed);
1983  result = this->basis_.interpolate(coords,ed);
1984  dmin = (result-p).length2();
1985  }
1986 
1987  pdist = sqrt(dmin);
1988  return (true);
1989 }
1990 
1991 
1992 } /// namespace SCIRun
1993 
1994 #endif /// SCI_project_StructQuadSurfMesh_h
SCIRun::StructQuadSurfMesh::jacobian_metric
double jacobian_metric(const INDEX idx) const
Definition: StructQuadSurfMesh.h:492
SCIRun::ImageMesh::Cell::array_type
std::vector< index_type > array_type
Definition: ImageMesh.h:278
SCIRun::ImageMesh::begin
void begin(typename Node::iterator &) const
Definition: ImageMesh.h:1543
MeshSupport.h
SCIRun::StructQuadSurfMesh::ElemData::node0_index
index_type node0_index() const
the following designed to coordinate with ::get_nodes
Definition: StructQuadSurfMesh.h:271
SCIRun::Piostream::reading
bool reading() const
Definition: Persistent.h:164
SCIRun::StructQuadSurfMesh::get_normal
void get_normal(Core::Geometry::Vector &result, std::vector< double > &coords, typename ImageMesh< Basis >::Elem::index_type eidx, unsigned int)
Definition: StructQuadSurfMesh.h:175
SCIRun::StructQuadSurfMesh::face_type_description
static const TypeDescription * face_type_description()
Definition: StructQuadSurfMesh.h:1711
SCIRun::StructQuadSurfMesh::locate
bool locate(typename ImageMesh< Basis >::Elem::index_type &elem, ARRAY &coords, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:210
SCIRun::StructQuadSurfMesh::get_epsilon
double get_epsilon() const
Definition: StructQuadSurfMesh.h:975
Point.h
SCIRun::FieldRNG
Definition: FieldRNG.h:37
SCIRun::TypeDescription::get_name
std::string get_name(const std::string &type_sep_start="<", const std::string &type_sep_end="> ") const
Definition: TypeDescription.cc:135
SCIRun::StructQuadSurfMesh::ElemData
Definition: StructQuadSurfMesh.h:259
SCIRun::Mesh::LOCATE_E
Definition: Mesh.h:85
SCIRun::StructQuadSurfMesh::~StructQuadSurfMesh
virtual ~StructQuadSurfMesh()
Definition: StructQuadSurfMesh.h:111
SCIRun::Core::Geometry::Transform::project
Vector project(const Vector &p) const
Definition: Transform.cc:425
SCIRun::ImageMesh::end
void end(typename Node::iterator &) const
Definition: ImageMesh.h:1551
SCIRun::Mesh::IRREGULAR
Definition: Mesh.h:64
SCIRun::StructQuadSurfMesh::locate_edge
bool locate_edge(INDEX &idx, const Core::Geometry::Point &p) const
This is currently implemented as an exhaustive search.
Definition: StructQuadSurfMesh.h:666
SCIRun::Core::Geometry::BBox::max
Point max() const
Definition: BBox.h:195
SCIRun::StructQuadSurfMesh::get_coords
bool get_coords(VECTOR &coords, const Core::Geometry::Point &p, typename ImageMesh< Basis >::Elem::index_type idx) const
Definition: StructQuadSurfMesh.h:367
ASSERTMSG
#define ASSERTMSG(condition, message)
Definition: Exception.h:113
SCIRun::StructQuadSurfMesh::get_length
double get_length(typename ImageMesh< Basis >::Edge::index_type idx) const
Definition: StructQuadSurfMesh.h:165
SCIRun::StructQuadSurfMesh::elem_type_description
static const TypeDescription * elem_type_description()
Definition: StructQuadSurfMesh.h:961
SCIRun::StructQuadSurfMesh::compute_elem_grid
void compute_elem_grid(Core::Geometry::BBox &bb)
Definition: StructQuadSurfMesh.h:1479
SCIRun::ImageMesh::ni_
size_type ni_
Definition: ImageMesh.h:942
SCIRun::StructQuadSurfMesh::compute_epsilon
void compute_epsilon(Core::Geometry::BBox &bb)
Definition: StructQuadSurfMesh.h:1515
SCIRun::StructQuadSurfMesh::insert_node_into_grid
void insert_node_into_grid(typename ImageMesh< Basis >::Node::index_type ci)
Definition: StructQuadSurfMesh.h:1426
SCIRun::StructQuadSurfMesh::next
index_type next(index_type i)
Definition: StructQuadSurfMesh.h:999
SCIRun::StructQuadSurfMesh::edge_type_description
static const TypeDescription * edge_type_description()
Definition: StructQuadSurfMesh.h:1693
SCIRun::ImageMesh
Definition: ImageMesh.h:61
SCIRun::Mesh::FIND_CLOSEST_E
Definition: Mesh.h:90
SCIRun::StructQuadSurfMesh::locate
bool locate(typename ImageMesh< Basis >::Cell::index_type &, const Core::Geometry::Point &) const
Definition: StructQuadSurfMesh.h:205
SCIRun::StructQuadSurfMesh::locate
bool locate(typename ImageMesh< Basis >::Node::index_type &node, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:196
SCIRun::Piostream
Definition: Persistent.h:89
SCIRun::StructQuadSurfMesh::index_type
SCIRun::index_type index_type
Definition: StructQuadSurfMesh.h:100
SCIRun::StructQuadSurfMesh::cell_type_description
static const TypeDescription * cell_type_description()
Definition: StructQuadSurfMesh.h:1729
SCIRun::StructQuadSurfMesh::ElemData::edge3_index
index_type edge3_index() const
Definition: StructQuadSurfMesh.h:310
SCIRun::ImageMesh::IFaceIter
Definition: ImageMesh.h:191
SCIRun::StructQuadSurfMesh::locate_elem
bool locate_elem(INDEX &elem, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:697
SCIRun::ImageMesh::get_nodes
void get_nodes(typename Node::array_type &, typename Edge::index_type) const
get the child elements of the given index
Definition: ImageMesh.h:1155
SCIRun::StructQuadSurfMesh::transform
virtual void transform(const Core::Geometry::Transform &t)
Definition: StructQuadSurfMesh.h:1117
SCIRun::Core::Geometry::Point::y
void y(const double)
Definition: Point.h:135
SCIRun::StructQuadSurfMesh::derivate
void derivate(const VECTOR1 &coords, typename ImageMesh< Basis >::Elem::index_type idx, VECTOR2 &J) const
Definition: StructQuadSurfMesh.h:389
SCIRun::ImageMesh::ImageIndex::j_
size_type j_
Definition: ImageMesh.h:119
SCIRun::Persistent
Definition: Persistent.h:187
SCIRun::SearchGridT::lookup_ijk
void lookup_ijk(iterator &begin, iterator &end, size_type i, size_type j, size_type k)
Definition: SearchGridT.h:207
SCIRun::Core::Geometry::BBox::extend
BBox & extend(const Point &p)
Expand the bounding box to include point p.
Definition: BBox.h:98
SCIRun::StructQuadSurfMesh::get_point
void get_point(Core::Geometry::Point &point, const typename ImageMesh< Basis >::Node::index_type &index) const
Definition: StructQuadSurfMesh.h:249
SCIRun::StructQuadSurfMesh::get_center
void get_center(Core::Geometry::Point &, typename ImageMesh< Basis >::Cell::index_type) const
Definition: StructQuadSurfMesh.h:192
SCIRun::TypeDescription::MESH_E
Definition: TypeDescription.h:50
TDCSInputData::mesh
FieldHandle mesh()
Definition: BuildTDCSMatrixTests.cc:56
SCIRun::Mesh::ELEM_LOCATE_E
Definition: Mesh.h:84
SCIRun::StructQuadSurfMesh::ElemData::node0
const Core::Geometry::Point & node0() const
Definition: StructQuadSurfMesh.h:318
SCIRun::Core::Geometry::Point
Definition: Point.h:49
SCIRun::TypeDescription
Definition: TypeDescription.h:45
SCISHARE
#define SCISHARE
Definition: share.h:39
SCIRun::TypeDescription::td_vec
std::vector< const TypeDescription * > td_vec
Definition: TypeDescription.h:56
SCIRun::StructQuadSurfMesh::synchronized_
mask_type synchronized_
Definition: StructQuadSurfMesh.h:1009
SCIRun::Mesh
Definition: Mesh.h:45
SCIRun::Core::Geometry::BBox
Definition: BBox.h:46
SCIRun::ImageMesh::basis_
Basis basis_
The basis class.
Definition: ImageMesh.h:951
SCIRun::MeshHandle
boost::shared_ptr< Mesh > MeshHandle
Definition: DatatypeFwd.h:67
SCIRun::Core::Thread::Mutex
Definition: Mutex.h:43
SCIRun::StructQuadSurfMesh::points_
Array2< Core::Geometry::Point > points_
Definition: StructQuadSurfMesh.h:1002
SCIRun::StructQuadSurfMesh::find_closest_elems
bool find_closest_elems(double &pdist, Core::Geometry::Point &result, std::vector< typename ImageMesh< Basis >::Elem::index_type > &elems, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:1875
SCIRun::Mesh::STRUCTURED
Definition: Mesh.h:61
SCIRun::StructQuadSurfMesh::synchronize
virtual bool synchronize(mask_type sync)
Definition: StructQuadSurfMesh.h:1325
SCIRun::StructQuadSurfMesh::VStructQuadSurfMesh
friend class VStructQuadSurfMesh
Make sure the virtual interface has access.
Definition: StructQuadSurfMesh.h:95
SCIRun::mask_type
unsigned int mask_type
Definition: Types.h:45
SCIRun::StructQuadSurfMesh::find_closest_elem
bool find_closest_elem(double &pdist, Core::Geometry::Point &result, typename ImageMesh< Basis >::Elem::index_type &elem, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:1863
SCIRun::StructQuadSurfMesh::node_grid_
boost::shared_ptr< SearchGridT< typename ImageMesh< Basis >::Node::index_type > > node_grid_
normalized per node
Definition: StructQuadSurfMesh.h:1005
SCIRun::Core::Geometry::Vector::asPoint
Point & asPoint() const
Definition: Vector.h:457
ASSERT
#define ASSERT(condition)
Definition: Assert.h:110
Plane.h
SCIRun::Piostream::begin_class
virtual int begin_class(const std::string &name, int current_version)
Definition: Persistent.cc:143
SCIRun::StructQuadSurfMesh::prev
index_type prev(index_type i)
Definition: StructQuadSurfMesh.h:1000
SCIRun::ImageMesh::vmesh_
boost::shared_ptr< VMesh > vmesh_
Virtual mesh.
Definition: ImageMesh.h:954
SCIRun::ImageMesh::IFaceIndex
Definition: ImageMesh.h:124
SCIRun::Array2::dim1
size_t dim1() const
Returns number of rows.
Definition: Array2.h:97
SCIRun::Core::Geometry::Vector
Definition: Vector.h:63
SCIRun::StructQuadSurfMesh::ElemData::edge0_index
index_type edge0_index() const
the following designed to coordinate with ::get_edges
Definition: StructQuadSurfMesh.h:294
SCIRun::DetMatrix3P
double DetMatrix3P(const VectorOfPoints &p)
Inline templated determinant of matrix.
Definition: Locate.h:120
SCIRun::StructQuadSurfMesh::get_random_point
void get_random_point(Core::Geometry::Point &, const typename ImageMesh< Basis >::Elem::index_type &, FieldRNG &rng) const
Definition: StructQuadSurfMesh.h:1286
SCIRun::find_type_name
const string find_type_name(float *)
Definition: TypeName.cc:63
SCIRun::Core::Geometry::Point::z
void z(const double)
Definition: Point.h:145
SCIRun::Mesh::EDGES_E
Definition: Mesh.h:72
SCIRun::Core::Geometry::BBox::valid
bool valid() const
Definition: BBox.h:93
SCIRun::StructQuadSurfMesh::remove_elem_from_grid
void remove_elem_from_grid(typename ImageMesh< Basis >::Elem::index_type ci)
Definition: StructQuadSurfMesh.h:1412
SCIRun::StructQuadSurfMesh::scaled_jacobian_metric
double scaled_jacobian_metric(const INDEX idx) const
Definition: StructQuadSurfMesh.h:463
SCIRun::DetMatrix3x3
T DetMatrix3x3(const T *p)
Definition: Locate.h:95
SCIRun::Mesh::EPSILON_E
Definition: Mesh.h:86
SCIRun::StructQuadSurfMesh::find_closest_node
bool find_closest_node(double &pdist, Core::Geometry::Point &result, typename ImageMesh< Basis >::Node::index_type &node, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:216
SCIRun::StructQuadSurfMesh::find_closest_elem
bool find_closest_elem(double &pdist, Core::Geometry::Point &result, ARRAY &coords, typename ImageMesh< Basis >::Elem::index_type &elem, const Core::Geometry::Point &p, double maxdist) const
Definition: StructQuadSurfMesh.h:780
SCIRun::SearchGridT
Definition: SearchGridT.h:47
SCIRun::StructQuadSurfMesh::ElemData::node1
const Core::Geometry::Point & node1() const
Definition: StructQuadSurfMesh.h:323
SCIRun::StructQuadSurfMesh::ElemData::node2_index
index_type node2_index() const
Definition: StructQuadSurfMesh.h:281
SCIRun::StructQuadSurfMesh::get_volume
double get_volume(typename ImageMesh< Basis >::Cell::index_type) const
Definition: StructQuadSurfMesh.h:169
name
const char * name[]
Definition: BoostGraphExampleTests.cc:87
SCIRun::Core::Geometry::Point::x
void x(const double)
Definition: Point.h:125
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::CreateVStructQuadSurfMesh
VMesh * CreateVStructQuadSurfMesh(MESH *mesh)
Definition: StructQuadSurfMesh.h:76
SCIRun::StructQuadSurfMesh::mask_type
SCIRun::mask_type mask_type
Definition: StructQuadSurfMesh.h:102
ImageMesh.h
SCIRun::size_type
long long size_type
Definition: Types.h:40
SCIRun::StructQuadSurfMesh::clear_synchronization
bool clear_synchronization()
Definition: StructQuadSurfMesh.h:1379
SCIRun::StructQuadSurfMesh::ElemData::ElemData
ElemData(const StructQuadSurfMesh< Basis > &msh, const typename ImageMesh< Basis >::Elem::index_type ind)
Definition: StructQuadSurfMesh.h:263
SCIRun::StackVector
Definition: StackVector.h:50
SCIRun::Core::Geometry::Cross
Vector Cross(const Vector &v1, const Vector &v2)
Definition: Vector.h:378
SCIRun::ScaledDetMatrix3P
double ScaledDetMatrix3P(const VectorOfPoints &p)
Inline templated determinant of matrix.
Definition: Locate.h:132
SCIRun::StructQuadSurfMesh::compute_normals
void compute_normals()
Definition: StructQuadSurfMesh.h:1524
SCIRun::ImageMesh::Face::array_type
std::vector< index_type > array_type
Definition: ImageMesh.h:271
SCIRun::Mesh::NODES_E
Definition: Mesh.h:71
SCIRun::StructQuadSurfMesh::ElemData::node3_index
index_type node3_index() const
Definition: StructQuadSurfMesh.h:287
SCIRun::Core::Geometry::BBox::diagonal
Vector diagonal() const
Definition: BBox.h:198
SCIRun::StructQuadSurfMesh::get_size
double get_size(typename ImageMesh< Basis >::Edge::index_type idx) const
Definition: StructQuadSurfMesh.h:142
SCIRun::StructQuadSurfMesh::type_idsqs
static PersistentTypeID type_idsqs
Definition: StructQuadSurfMesh.h:946
SCIRun::StructQuadSurfMesh::synchronize_lock_
Core::Thread::Mutex synchronize_lock_
Definition: StructQuadSurfMesh.h:1008
SCIRun::StructQuadSurfMesh::locate
bool locate(typename ImageMesh< Basis >::Face::index_type &face, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:202
SCIRun::StructQuadSurfMesh::normals_
Array2< Core::Geometry::Vector > normals_
Definition: StructQuadSurfMesh.h:1003
SCIRun::ImageMesh::nj_
size_type nj_
Definition: ImageMesh.h:943
SCIRun::Mesh::NORMALS_E
Definition: Mesh.h:79
SCIRun::StructQuadSurfMesh::pwl_approx_edge
void pwl_approx_edge(std::vector< std::vector< double > > &coords, typename ImageMesh< Basis >::Elem::index_type, unsigned which_edge, unsigned div_per_unit) const
Definition: StructQuadSurfMesh.h:352
SCIRun::ImageMesh::io
virtual void io(Piostream &)
Export this class using the old Pio system.
Definition: ImageMesh.h:1479
SCIRun::StructQuadSurfMesh::type_name
static const std::string type_name(int n=-1)
Core functionality for getting the name of a templated mesh class.
Definition: StructQuadSurfMesh.h:1626
CompGeom.h
SCIRun::Core::Geometry::BBox::min
Point min() const
Definition: BBox.h:192
SCIRun::Core::Geometry::Transform
Definition: Transform.h:53
SCIRun::StructQuadSurfMesh::structquadsurf_maker
static MeshHandle structquadsurf_maker(size_type x, size_type y)
This function returns a handle for the virtual interface.
Definition: StructQuadSurfMesh.h:969
SCIRun::Core::Geometry::Vector::x
void x(double)
Definition: Vector.h:175
STRUCT_QUAD_SURF_MESH_VERSION
#define STRUCT_QUAD_SURF_MESH_VERSION
Definition: StructQuadSurfMesh.h:1588
SCIRun::StructQuadSurfMesh::get_size
double get_size(typename ImageMesh< Basis >::Cell::index_type) const
Definition: StructQuadSurfMesh.h:163
SCIRun::Mesh::FACES_E
Definition: Mesh.h:73
SCIRun::TypeNameGenerator::make_template_id
static const std::string make_template_id(const std::string &templateName, const std::string &templateParam)
Definition: TypeName.h:62
SCIRun::StructQuadSurfMesh::inside
bool inside(INDEX idx, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:521
SCIRun::StructQuadSurfMesh::get_bounding_box
virtual Core::Geometry::BBox get_bounding_box() const
get the mesh statistics
Definition: StructQuadSurfMesh.h:1097
SCIRun::ImageMesh::Edge::array_type
std::vector< index_type > array_type
Definition: ImageMesh.h:264
SCIRun::StructQuadSurfMesh
Declarations for virtual interface.
Definition: StructQuadSurfMesh.h:70
readAllFields.v
v
Definition: readAllFields.py:42
SCIRun::Core::Geometry::Vector::normalize
double normalize()
Definition: Vector.h:437
SCIRun::StructQuadSurfMesh::get_size
double get_size(const typename ImageMesh< Basis >::Face::index_type &idx) const
Definition: StructQuadSurfMesh.h:151
SCIRun::StructQuadSurfMesh::get_center
void get_center(Core::Geometry::Point &, const typename ImageMesh< Basis >::Node::index_type &) const
get the center point (in object space) of an element
Definition: StructQuadSurfMesh.h:1148
SCIRun::ImageMesh::INodeIndex
Definition: ImageMesh.h:136
SCIRun::StructQuadSurfMesh::elem_grid_
boost::shared_ptr< SearchGridT< typename ImageMesh< Basis >::Elem::index_type > > elem_grid_
Definition: StructQuadSurfMesh.h:1006
SCIRun::ImageMesh::IFaceSize
Definition: ImageMesh.h:245
share.h
SCIRun::StructQuadSurfMesh::compute_node_grid
void compute_node_grid(Core::Geometry::BBox &bb)
Definition: StructQuadSurfMesh.h:1443
SCIRun::Core::Geometry::Vector::y
void y(double)
Definition: Vector.h:185
SCIRun::StructQuadSurfMesh::ElemData::node3
const Core::Geometry::Point & node3() const
Definition: StructQuadSurfMesh.h:333
SCIRun::StructQuadSurfMesh::jacobian
void jacobian(const VECTOR &coords, typename ImageMesh< Basis >::Elem::index_type idx, double *J) const
Definition: StructQuadSurfMesh.h:413
SCIRun::StructQuadSurfMesh::ElemData::edge1_index
index_type edge1_index() const
Definition: StructQuadSurfMesh.h:299
SCIRun::Core::Geometry::Vector::safe_normalize
double safe_normalize()
Definition: Vector.h:236
SCIRun::Array2
Definition: Array2.h:51
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::Mesh::ELEMS_E
Definition: Mesh.h:75
ASSERTFAIL
#define ASSERTFAIL(string)
Definition: Assert.h:52
SCIRun::StructQuadSurfMesh::locate
bool locate(typename ImageMesh< Basis >::Edge::index_type &edge, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:199
SCIRun::StructQuadSurfMesh::get_weights
int get_weights(const Core::Geometry::Point &, typename ImageMesh< Basis >::Cell::array_type &, double *)
Definition: StructQuadSurfMesh.h:240
SCIRun::Pio
void Pio(Piostream &stream, Array1< T > &array)
Definition: Array1.h:65
SCIRun::StructQuadSurfMesh::det_jacobian
double det_jacobian(const VECTOR &coords, typename ImageMesh< Basis >::Elem::index_type idx) const
Definition: StructQuadSurfMesh.h:401
SCIRun::ImageMesh::size
void size(typename Node::size_type &) const
Definition: ImageMesh.h:1559
Mutex.h
SCIRun::StructQuadSurfMesh::locate_elem
bool locate_elem(INDEX &elem, ARRAY &coords, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:732
SCIRun::StructQuadSurfMesh::topology_geometry
virtual int topology_geometry() const
Definition: StructQuadSurfMesh.h:134
SCIRun::Core::Geometry::Vector::length
double length() const
Definition: Vector.h:365
SCIRun::Piostream::end_class
virtual void end_class()
Definition: Persistent.cc:178
SCIRun::StructQuadSurfMesh::StructQuadSurfMesh
StructQuadSurfMesh()
Definition: StructQuadSurfMesh.h:1024
SCIRun::StructQuadSurfMesh::remove_node_from_grid
void remove_node_from_grid(typename ImageMesh< Basis >::Node::index_type ci)
Definition: StructQuadSurfMesh.h:1436
SCIRun::index_type
long long index_type
Definition: Types.h:39
SCIRun::StructQuadSurfMesh::mesh_maker
static MeshHandle mesh_maker()
This function returns a handle for the virtual interface.
Definition: StructQuadSurfMesh.h:967
SCIRun::ImageMesh::ImageIndex::i_
size_type i_
Definition: ImageMesh.h:119
SCIRun::StructQuadSurfMesh::epsilon_
double epsilon_
Definition: StructQuadSurfMesh.h:1010
SCIRun::StructQuadSurfMesh::get_type_description
virtual const TypeDescription * get_type_description() const
Definition: StructQuadSurfMesh.h:1667
SCIRun::StructQuadSurfMesh::insert_elem_into_grid
void insert_elem_into_grid(typename ImageMesh< Basis >::Elem::index_type ci)
Used to recompute data for individual cells.
Definition: StructQuadSurfMesh.h:1395
SCIRun::Mesh::NODE_LOCATE_E
Definition: Mesh.h:83
SCIRun::StackVector::resize
void resize(size_t size, const value_type &val=value_type())
Definition: StackVector.h:61
SCIRun::StructQuadSurfMesh::point
const Core::Geometry::Point & point(const typename ImageMesh< Basis >::Node::index_type &idx)
Definition: StructQuadSurfMesh.h:996
SearchGridT.h
SCIRun::ImageMesh::ImageIndex::mesh_
const ImageMesh * mesh_
Definition: ImageMesh.h:121
SCIRun::StructQuadSurfMesh::set_point
void set_point(const Core::Geometry::Point &point, const typename ImageMesh< Basis >::Node::index_type &index)
Definition: StructQuadSurfMesh.h:1316
SCIRun::StructQuadSurfMesh::get_area
double get_area(const typename ImageMesh< Basis >::Face::index_type &idx) const
Definition: StructQuadSurfMesh.h:167
SCIRun::Core::Geometry::BBox::reset
void reset()
Definition: BBox.h:95
SCIRun::StructQuadSurfMesh::ElemData::edge2_index
index_type edge2_index() const
Definition: StructQuadSurfMesh.h:304
SCIRun::StructQuadSurfMesh::handle_type
boost::shared_ptr< StructQuadSurfMesh< Basis > > handle_type
Definition: StructQuadSurfMesh.h:104
SCIRun::PersistentTypeID
Definition: Persistent.h:64
SCIRun::ImageMesh::size_type
SCIRun::size_type size_type
Definition: ImageMesh.h:93
SCIRun::ImageMesh::INodeSize
Definition: ImageMesh.h:238
SCIRun::ImageMesh::mask_type
SCIRun::mask_type mask_type
Definition: ImageMesh.h:94
SCIRun::StructQuadSurfMesh::get_dim
virtual bool get_dim(std::vector< size_type > &) const
Definition: StructQuadSurfMesh.h:1083
SCIRun::StructQuadSurfMesh::clone
virtual StructQuadSurfMesh * clone() const
Definition: StructQuadSurfMesh.h:109
SCIRun::EdgeIndex
Distinct type for edge index.
Definition: FieldIndex.h:81
SCIRun::StructQuadSurfMesh::get_weights
int get_weights(const Core::Geometry::Point &, typename ImageMesh< Basis >::Edge::array_type &, double *)
Definition: StructQuadSurfMesh.h:235
eab.n
int n
Definition: eab.py:9
SCIRun::Array2::dim2
size_t dim2() const
Returns number of cols.
Definition: Array2.h:101
SCIRun::InverseMatrix3x3
T InverseMatrix3x3(const T *p, T *q)
Definition: Locate.h:47
SCIRun::StructQuadSurfMesh::inside3_p
bool inside3_p(typename ImageMesh< Basis >::Face::index_type i, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:1232
SCIRun::StructQuadSurfMesh::ElemData::node4
const Core::Geometry::Point & node4() const
Definition: StructQuadSurfMesh.h:338
SCIRun::ImageMesh::INodeIter
Definition: ImageMesh.h:162
SCIRun::StructQuadSurfMesh::under_type
SCIRun::index_type under_type
Types that change depending on 32 or 64bits.
Definition: StructQuadSurfMesh.h:99
DEBUG_CONSTRUCTOR
#define DEBUG_CONSTRUCTOR(type)
Definition: Debug.h:64
SCIRun::StructQuadSurfMesh::io
virtual void io(Piostream &)
Export this class using the old Pio system.
Definition: StructQuadSurfMesh.h:1592
SCIRun::StructQuadSurfMesh::get_size
double get_size(const typename ImageMesh< Basis >::Node::index_type &) const
Get the size of an element (length, area, volume)
Definition: StructQuadSurfMesh.h:140
SCIRun::StructQuadSurfMesh::get_points
Array2< Core::Geometry::Point > & get_points()
Definition: StructQuadSurfMesh.h:971
SCIRun::StructQuadSurfMesh::find_closest_elem
bool find_closest_elem(double &pdist, Core::Geometry::Point &result, ARRAY &coords, typename ImageMesh< Basis >::Elem::index_type &elem, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:773
SCIRun::StructQuadSurfMesh::inverse_jacobian
double inverse_jacobian(const VECTOR &coords, typename ImageMesh< Basis >::Elem::index_type idx, double *Ji) const
Definition: StructQuadSurfMesh.h:438
SCIRun::StructQuadSurfMesh::interpolate
void interpolate(Core::Geometry::Point &pt, const VECTOR &coords, typename ImageMesh< Basis >::Elem::index_type idx) const
Definition: StructQuadSurfMesh.h:377
SCIRun::Core::Geometry::Vector::z
void z(double)
Definition: Vector.h:195
SCIRun::StructQuadSurfMesh::node_type_description
static const TypeDescription * node_type_description()
Definition: StructQuadSurfMesh.h:1675
SCIRun::StructQuadSurfMesh::maker
static Persistent * maker()
This function returns a maker for Pio.
Definition: StructQuadSurfMesh.h:965
SCIRun::VMesh
Definition: VMesh.h:53
SCIRun::StructQuadSurfMesh::get_normal
void get_normal(Core::Geometry::Vector &, const typename ImageMesh< Basis >::Node::index_type &) const
Definition: StructQuadSurfMesh.h:1139
SCIRun::StructQuadSurfMesh::set_dim
virtual void set_dim(std::vector< size_type > dims)
Definition: StructQuadSurfMesh.h:121
SCIRun::StructQuadSurfMesh::epsilon2_
double epsilon2_
Definition: StructQuadSurfMesh.h:1011
eabLatVolData.size
int size
Definition: eabLatVolData.py:2
SCIRun::StructQuadSurfMesh::locate_node
bool locate_node(INDEX &node, const Core::Geometry::Point &p) const
Definition: StructQuadSurfMesh.h:575
SCIRun::StructQuadSurfMesh::size_type
SCIRun::size_type size_type
Definition: StructQuadSurfMesh.h:101
DEBUG_DESTRUCTOR
#define DEBUG_DESTRUCTOR(type)
Definition: Debug.h:65
SCIRun::get_type_description
const TypeDescription * get_type_description(Core::Basis::ConstantBasis< T > *)
Definition: Constant.h:209
SCIRun::StructQuadSurfMesh::ElemData::node1_index
index_type node1_index() const
Definition: StructQuadSurfMesh.h:276
SCIRun::StructQuadSurfMesh::unsynchronize
virtual bool unsynchronize(mask_type sync)
Definition: StructQuadSurfMesh.h:1371
Array2.h
Interface to dynamic 2D array class.
SCIRun::StructQuadSurfMesh::get_weights
int get_weights(const Core::Geometry::Point &p, typename ImageMesh< Basis >::Node::array_type &l, double *w)
Definition: StructQuadSurfMesh.h:1212
SCIRun::elem_locate
bool elem_locate(INDEX &elem, MESH &msh, const Core::Geometry::Point &p)
General case locate, search each elem.
Definition: Mesh.h:188
SCIRun::StructQuadSurfMesh::ElemData::node2
const Core::Geometry::Point & node2() const
Definition: StructQuadSurfMesh.h:328
Generated on Wed Jul 9 2014 14:31:45 for SCIRun by   doxygen 1.8.5