HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SIM_RawField.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: SIM_RawField.h ( SIM Library, C++)
7  *
8  * COMMENTS:
9  */
10 
11 #ifndef __SIM_RawField__
12 #define __SIM_RawField__
13 
14 #include "SIM_API.h"
15 
16 #include <GEO/GEO_PointTree.h>
17 #include <UT/UT_VoxelArray.h>
19 #include <UT/UT_Vector.h>
20 #include <UT/UT_FloatArray.h>
21 #include <UT/UT_Wavelet.h>
22 
23 class SIM_VectorField;
24 class SIM_RawField;
25 class SIM_RawIndexField;
26 class GEO_PrimVolume;
27 class GEO_PrimVDB;
29 class GU_Detail;
30 class CE_Grid;
31 
32 class gu_sdf_qelem;
33 
34 // The order of this enum is important. It is a 3 bit array
35 // where each bit determines if the center value or corner value
36 // should be used on that axis. 0 means center, 1 means node.
37 // The least significant bit is x axis, then y, then z.
39 {
48 };
49 
51 {
52  SIM_BOUNDARY_NONE, // Unchanged, free to change.
53  SIM_BOUNDARY_SAME, // Same as adjacent voxels.
54  SIM_BOUNDARY_ZERO, // Set to constant value, myBoundaryValue
55  SIM_BOUNDARY_NEGATE, // Opposite sign of adjacent voxels.
56  SIM_BOUNDARY_FIXED, // Unaffected by enforcement, but should be treated
57  // as a fixed boundary by projection methods.
58 };
59 
61 {
62  SIM_ADVECT_SINGLE = 0, // Take a single step
63  SIM_ADVECT_TRACE = 1, // Trace according to velocity field
64  SIM_ADVECT_MIDPOINT = 2, // Trace using midpoint updates
65  SIM_ADVECT_HJWENO = 3, // Don't trace, but use HJWENO update
66  SIM_ADVECT_UPWIND = 4, // First order upwind.
67  SIM_ADVECT_RK3 = 5, // Third order TVD-RK3
68  SIM_ADVECT_RK4 = 6 // Fourth order TVD-RK4
69 };
70 
71 //
72 // The boundary of the field can be set to be "open" above a
73 // line using the following struct. The height value specifies
74 // the dividing line between the open and closed boundary with
75 // respect to the plane positioned at the origin and oriented with
76 // planeNormal. This is useful for letting liquid splash
77 // out of a closed boundary without the expense of explicitly
78 // creating a solid boundary.
79 //
80 
82 {
84  : activated(false)
85  , height(0.)
86  , planeNormal(UT_Vector3(0.))
87  {}
88  bool activated;
91 };
92 
93 
94 // Type definition for a velocity function
96 
98 {
99 public:
100  SIM_RawField();
101  virtual ~SIM_RawField();
102 
103  /// Copy constructor:
104  SIM_RawField(const SIM_RawField &src);
105 
106  /// Assigment operator:
107  const SIM_RawField &operator=(const SIM_RawField &src);
108 
109  /// Initializes the field.
110  /// The resolution given is in terms of voxels, the actual dimensions
111  /// of this field may be slightly different due to the sampling
112  /// choice.
113  void init(SIM_FieldSample sample,
114  const UT_Vector3 &orig, const UT_Vector3 &size,
115  int xres, int yres, int zres);
116 
117  /// Initializes the field.
118  /// Will gain ownership of the given voxel array.
119  void init(SIM_FieldSample sample,
120  const UT_Vector3 &orig, const UT_Vector3 &size,
121  UT_VoxelArrayF *voxels);
122 
123  /// Initializes the voxel iterator to only run over a subset
124  /// of tiles that correspond to the appropriate range for
125  /// the jobinfo.
126  void getPartialRange(UT_VoxelArrayIteratorF &vit,
127  const UT_JobInfo &info) const;
128 
129  /// Returns true if this should be multithreaded.
130  bool shouldMultiThread() const
131  {
132 #ifdef CELLBE
133  return false;
134 #else
135  return field()->numTiles() > 1;
136 #endif
137  }
138 
139  /// Initializes this to be the same dimensions, sampling pattern,
140  /// etc, of the given field.
141  /// The values of this may be reset to zero.
142  void match(const SIM_RawField &src);
143  void match(const SIM_RawIndexField &src);
144 
145  /// Initializes this field to a constant value.
146  void makeConstant(float cval);
147 
148  /// Initializes this to be the same as the given field. The
149  /// values of this will be found by resampling the other field
150  /// with the new sampling pattern.
151  void resample(SIM_FieldSample sample, const SIM_RawField *src);
152 
153  /// Returns the set of samples in *this* field which correspond
154  /// to the given location & sampling pattern.
155  /// ix, iy, and iz should be size 8.
156  /// If you want the deltas for the sampling pattern, call with
157  /// x, y, z zero and clamp to false.
158  void getSamplePattern(SIM_FieldSample sample, int x, int y, int z,
159  int &numsample, int *ix, int *iy, int *iz,
160  bool clamp) const;
161 
162  /// Returns an offset, in voxel coords, to sample *this* field giving
163  /// a voxel based address in a field with sampling pattern sample.
164  /// The offset will be 0.5, 0, or -0.5 in each axis.
165  UT_Vector3 getSampleVoxelOffset(SIM_FieldSample sample) const;
166 
167  /// Extrapolates the values of this field to all locations which are
168  /// on the wrong side of the isocontour. dir == 1 means greater than
169  /// isocontour is wrong, dir == -1 is less than.
170  /// maxdist, if not negative, represents the greatest *absolute*
171  /// distance to extrapolate to. Values outside of this are left
172  /// unaffected, unless clamp is set, in which case they are set to
173  /// clampval.
174  void extrapolate(const SIM_RawField *depths,
175  const SIM_RawField *valid,
176  fpreal isocontour,
177  fpreal dir, fpreal maxdist,
178  bool clamp, fpreal clampval);
179 
181  {
187  const GU_Detail *gdp;
193  bool clamp;
195  };
196 
198  {
200  sampleval(v), samplew(w), ptpos(p) {}
201  sim_extrapolateCacheElem() : sampleval(0), samplew(0) {}
205  };
206 
207  /// Extrapolates the values of this field to all locations which are
208  /// on the wrong side of the isocontour. dir == 1 means greater than
209  /// isocontour is wrong, dir == -1 is less than.
210  /// maxdist, if not negative, represents the greatest *absolute*
211  /// distance to extrapolate to. Values outside of this are left
212  /// unaffected, unless clamp is set, in which case they are set to
213  /// clampval. The extrpolated values are found by interpolating at the
214  /// position of the closest points in the provided GU_Detail, as stored
215  /// in the index field by point index.
216  void extrapolateFromIndex(const SIM_RawField *depths,
217  const SIM_RawIndexField *index,
218  const SIM_RawField *collision,
219  const SIM_RawField *valid,
220  const GU_Detail *gdp,
221  const UT_DMatrix4 &xform,
222  fpreal isocontour,
223  fpreal isotol,
224  fpreal dir,
225  fpreal maxdist,
226  bool clamp,
227  fpreal clampval);
228 
229  THREADED_METHOD1(SIM_RawField, shouldMultiThread(),
230  extrapolateFromIndexInternal,
232 
233  void extrapolateFromIndexInternalPartial(sim_extrapolateFromIndexParms &p,
234  const UT_JobInfo &info);
235 
236  /// Advances source field along the field's curvature.
237  /// Takes a given b to act as the movement coefficient.
238  /// Does a single euler step.
239  /// this is filled in with the resulting field. this cannot be source.
240  /// This solves:
241  /// dphi/dt = b * K * |grad(phi)|
242  THREADED_METHOD3(SIM_RawField, shouldMultiThread(), moveAlongCurvature,
243  fpreal, b_val,
244  const SIM_RawField &, source,
245  fpreal, timestep)
246  void moveAlongCurvaturePartial(fpreal b_val,
248  fpreal timestep,
249  const UT_JobInfo &jobinfo);
250 
251  /// Advances this field along the field's curvature.
252  /// Is given a total amount to move and will calculate the appropriate
253  /// timestep according to the given cfl condition & min/max iterations.
254  void moveAlongCurvature(fpreal amount,
255  fpreal cflcond,
256  int miniter, int maxiter);
257 
258  /// Advances source field along the field's normal direction.
259  /// Takes a matching field which defines the speed to move in the
260  /// normal direction on a per cell basis.
261  /// Past the bandwidth assumes a valid SDF and skips gradient estimates.
262  /// Does a single euler step.
263  /// this is filled in with the resulting field. this cannot
264  /// be source.
265  /// This solves:
266  /// dphi/dt + speed * |grad(phi)| = 0
267  THREADED_METHOD4(SIM_RawField, shouldMultiThread(), moveAlongNormal,
268  const SIM_RawField &, speed,
269  const SIM_RawField &, source,
270  fpreal, timestep,
271  fpreal, bandwidth)
272  void moveAlongNormalPartial(const SIM_RawField &speed,
273  const SIM_RawField &source,
274  fpreal timestep, fpreal bandwidth,
275  const UT_JobInfo &info);
276 
277  /// Uses the midpoint method to do a second order temporal
278  /// update of the moveAlongNormal algorithm.
279  void moveAlongNormalMidpoint(const SIM_RawField &speed,
280  const SIM_RawField &source,
281  fpreal timestep, fpreal bandwidth);
282 
283  /// Performs the reinitialization equation.
284  /// This solves for:
285  /// dphi/dt + S(phi0) * (|grad(phi)| - 1) = 0
286  /// at steady state.
287  /// S(phi0) = phi0 / sqrt(phi0^2 + dx^2), a smeared sign function.
288  /// It is held constant for the reinitialization.
289  /// Usemidpoint indicates using midpoint integration vs. Euler.
290  /// Any parts of the SDF past the bandwidth are assumed to already
291  /// be a valid SDF and not require reinitialization.
292  void reinitializeSignedDistance(int maxiter,
293  bool usemidpoint,
294  fpreal bandwidth);
295 
296  /// Negates all values in this field
297  THREADED_METHOD(SIM_RawField, shouldMultiThread(), negate)
298  void negatePartial(const UT_JobInfo &info);
299 
300  /// Scales all values in this field
301  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), scale,
302  fpreal, scale)
303  void scalePartial(fpreal scale, const UT_JobInfo &info);
304 
305  /// Makes this field the minimum or maximum
306  /// of this field and the other field. The other field
307  /// is not assumed to be matching.
308  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), maximum,
309  const SIM_RawField *, other)
310  void maximumPartial(const SIM_RawField *other, const UT_JobInfo &info);
311 
312  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), minimum,
313  const SIM_RawField *, other)
314  void minimumPartial(const SIM_RawField *other, const UT_JobInfo &jobinfo);
315 
316  /// Averages this with the other field. Result in this. Assumes
317  /// fields match.
318  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), average,
319  const SIM_RawField &, other)
320  void averagePartial(const SIM_RawField &other, const UT_JobInfo &jobinfo);
321 
322  /// Sums the field into separate summation lists, one per component
323  /// of the comp index field. Negative components ignored, the given
324  /// arrays should already be sized to fit the maxindex.
325  THREADED_METHOD3_CONST(SIM_RawField, shouldMultiThread(), sumPerComponent,
326  UT_DoubleArray &, result,
327  UT_Int64Array &, activevoxels,
328  const SIM_RawIndexField *, comp)
329  void sumPerComponentPartial(UT_DoubleArray &result,
330  UT_Int64Array &activevoxels,
331  const SIM_RawIndexField *comp,
332  const UT_JobInfo &info) const;
333 
334  THREADED_METHOD5_CONST(SIM_RawField, shouldMultiThread(), sumComponent,
335  double &, result,
336  int64 &, activevoxels,
337  int, compnum,
338  const SIM_RawIndexField *, comp,
339  const SIM_RawField *, pressureboundary)
340  void sumComponentPartial(double &result,
341  int64 &activevoxels,
342  int compnum,
343  const SIM_RawIndexField *comp,
344  const SIM_RawField *pressureboundary,
345  const UT_JobInfo &info) const;
346 
347  /// Adds to this field the value given per component. Negative components
348  /// neglected.
349  THREADED_METHOD2(SIM_RawField, shouldMultiThread(), addValuePerComponent,
350  const UT_DoubleArray &, valuelist,
351  const SIM_RawIndexField *, comp)
352  void addValuePerComponentPartial(const UT_DoubleArray &valuelist,
353  const SIM_RawIndexField *comp,
354  const UT_JobInfo &info);
355 
356  THREADED_METHOD4(SIM_RawField, shouldMultiThread(), addValueToComponent,
357  double, value,
358  int, compnum,
359  const SIM_RawIndexField *, comp,
360  const SIM_RawField *, pressureboundary)
361  void addValueToComponentPartial(double value,
362  int compnum,
363  const SIM_RawIndexField *comp,
364  const SIM_RawField *pressureboundary,
365  const UT_JobInfo &info);
366 
367  /// Sets this to A + s*B. A & B match this. This can be either
368  /// A or B, or neither.
369  THREADED_METHOD3(SIM_RawField, shouldMultiThread(), setScaleAdd,
370  const SIM_RawField &, A,
371  fpreal, scale,
372  const SIM_RawField &, B)
373  void setScaleAddPartial(const SIM_RawField &A,
374  fpreal scale,
375  const SIM_RawField &B, const UT_JobInfo &jobinfo);
376 
377  /// Sets this to the smeared sign function of the given sdf.
378  /// this = sdf / sqrt(sdf*sdf + bandwidth*bandwidth)
379  /// this must already match sdf.
380  THREADED_METHOD2(SIM_RawField, shouldMultiThread(), smearedSign,
381  const SIM_RawField &, sdf,
382  fpreal, bandwidth)
383  void smearedSignPartial(const SIM_RawField &sdf,
384  fpreal bandwidth,
385  const UT_JobInfo &info);
386 
387  /// Sets this to the heaviside function of itself.
388  /// this = clamp(-this/diam+0.5, 0, 1)
389  THREADED_METHOD(SIM_RawField, shouldMultiThread(), convertToHeaviside)
390  void convertToHeavisidePartial(const UT_JobInfo &info);
391 
392  /// Perform heaviside & inverse heaviside in a consistent fashion.
393  static fpreal toHeaviside(fpreal val, fpreal diam);
394  static fpreal fromHeaviside(fpreal val, fpreal diam);
395  /// Defines width as maxcomponent of our voxels.
396  fpreal toHeaviside(fpreal val) const;
397  fpreal fromHeaviside(fpreal val) const;
398 
399  /// Types of reductions supported by reduceOp.
401  {
411  REDUCE_MINABS
412  };
413 
414  /// Perform reduction on the field using the given method
415  /// The internal methods are helpers to deal with the threading
416  /// collating.
417  fpreal reduceOp(REDUCE_NAMES op) const;
418 
419  /// Prebuilt versions.
420  fpreal average() const { return reduceOp(REDUCE_AVERAGE); }
421 
422  /// Performs a localized reduction on the field.
423  /// Stores in this the result of the reduction.
424  /// Radius is in world coordinates.
425  void localReduceOp(REDUCE_NAMES op, const UT_Vector3 &radius);
426 
427  /// Inplace box blur.
428  /// Radius is in world coordinates.
429  void boxBlur(float radius)
430  {
431  UT_Vector3 r(radius, radius, radius);
432  localReduceOp(REDUCE_AVERAGE, r);
433  }
434 
435  /// Performs the reduction only including parts of the field
436  /// that match the given mask. If mask is null, falls through
437  /// to normal reduction.
438  fpreal reduceMaskedOp(REDUCE_NAMES op,
439  const SIM_RawField *mask, bool maskissdf) const;
440 
441  /// Sorts the valid voxels (as defined by optional mask) and returns
442  /// the given percentile voxel. 0.5 means median. 0.25 first quartile,
443  /// etc.
444  fpreal findProportionalValue(fpreal position,
445  const SIM_RawField *mask, bool maskissdf) const;
446 
447  /// Returns true if the given field and this one match in terms
448  /// of number of voxels and bounding box size.
449  /// This means the voxel cells match - not necessarily the sample
450  /// points!
451  bool isMatching(const SIM_RawField *field) const;
452  bool isMatching(const SIM_RawIndexField *field) const;
453 
454  /// Returns true if the two fields are precisely aligned. This
455  /// means that samples are matched so a given integer index
456  /// into either field would give the same result.
457  bool isAligned(const SIM_RawField *field) const;
458  bool isAligned(const SIM_RawIndexField *field) const;
459  bool isAligned(const GEO_PrimVolume *vold) const;
460 
461  /// Returns true if our field has any NANs
462  bool hasNan() const { return field()->hasNan(); }
463 
464  /// Tests for nans, outputs text and asserts if has any
465  /// Only runs if the test environment variable is set.
466  void testForNan() const;
467 
468  /// Fetches the raw field.
469  const UT_VoxelArrayF *field() const { if (myFieldOutOfDate) { updateFieldFromGrid(); } return myField; }
470  UT_VoxelArrayF *fieldNC() const { if (myGrid) { clearGrid(); } return myField; }
471 
472  /// Fetches the compute grid, returns 0 if not present.
473  CE_Grid *grid() const;
474 
475  /// Fetches the compute grid but throws cl::Error() on failure.
476  CE_Grid *requireGrid() const;
477 
478  /// Mark the field as out of date, but only if we have a valid grid.
479  void markGridAsChanged() {if (myGrid) myFieldOutOfDate = true;}
480  void clearGrid() const;
481 
482  /// Copies our field into the grid.
483  void updateGridFromField() const;
484 
485  /// Copies our grid into our field.
486  void updateFieldFromGrid() const;
487 
488  /// Steals the voxel array, leaving this pointing to a 0 constant
489  /// array
490  UT_VoxelArrayF *steal();
491 
492  const UT_Vector3 &getOrig() const { return myOrig; }
493  const UT_Vector3 &getSize() const { return mySize; }
494  const UT_Vector3 &getBBoxOrig() const { return myBBoxOrig; }
495  const UT_Vector3 &getBBoxSize() const { return myBBoxSize; }
496  const UT_Vector3 &getVoxelSize() const { return myVoxelSize; }
497  fpreal getVoxelDiameter() const { return myVoxelDiameter; }
499  { return myVoxelSize.x()
500  * myVoxelSize.y()
501  * myVoxelSize.z(); }
502 
503  SIM_FieldSample getSample() const { return mySample; }
504 
505  int64 getMemoryUsage() const;
506 
507  /// Returns the resolution of the voxel grid that we are sampling.
508  /// This is a count of voxels, so may differ for our different
509  /// sampling methods.
510  void getVoxelRes(int &xres, int &yres, int &zres) const;
511 
512  /// Returns the actual number of samples in each resolution.
513  /// Preferred over field()->getXRes() as it doesn't require a
514  /// copy of the CE fields.
515  exint getXRes() const { return myField->getXRes(); }
516  exint getYRes() const { return myField->getYRes(); }
517  exint getZRes() const { return myField->getZRes(); }
518 
519  /// Functions to resolve quantities about the field.
520 
521  /// Sample the field at the given world space location.
522  fpreal getValue(UT_Vector3 pos) const;
523 
524  /// Returns an averaged value for the center of the given voxel.
525  fpreal getCellValue(int x, int y, int z) const;
526 
527  /// Adds the value to the given voxel cell, averaging out among
528  /// adjacent samples if we aren't sampled evenly.
529  void addToCell(int x, int y, int z, fpreal v);
530 
531  /// Ensures the given voxel cell has the given value. This
532  /// will set all of the adjacent samples if we aren't sampled
533  /// evenly.
534  void setCellValue(int x, int y, int z, fpreal v);
535 
536  /// Return the gradient of the field at the given world space location.
537  /// Uses central differencing with a sample spacing of the voxelsize.
538  UT_Vector3 getGradient(UT_Vector3 pos) const;
539 
540  /// Returns the gradient at the given voxel index.
541  /// Uses central differencing.
542  UT_Vector3 getGradientAtIndex(int x, int y, int z) const;
543 
544  /// Computes the laplacian of the field at the world space
545  /// coordinate using interpolation of getLaplacianAtIndex
546  fpreal64 getLaplacian(UT_Vector3 pos) const;
547 
548  /// Computes the laplacian of the field at the specific voxel index
549  /// specified
550  fpreal64 getLaplacianAtIndex(int x, int y, int z) const;
551 
552  /// Returns the curvature of the field at a given world space location.
553  /// Uses interpolation of index based getCurvatureAtIndex.
554  fpreal64 getCurvature(UT_Vector3 pos) const;
555 
556  /// Computes the curvature at the given voxel index.
557  /// Uses central differencing.
558  /// The resulting curvature is clamped according to the resolution
559  /// of the field to avoid abnormally large values in noisy data.
560  fpreal64 getCurvatureAtIndex(int x, int y, int z) const;
561 
562  /// Computes the curvature according to a 3^3 voxel probe
563  static fpreal64 getCurvatureAtProbe(UT_VoxelProbeCubeF &probe, const UT_Vector3 &invvoxelsize);
564 
565  /// Computes K |grad(phi|), the curvature times the length of the
566  /// gradient. By folding the two operations, we can return a non-zero
567  /// value where grad goes to zero (but curvature doesn't) using
568  /// L'Hospital's rule. This also does not clamp the curvature.
569  fpreal64 getCurvatureTimesGradAtIndex(int x, int y, int z) const;
570 
571  /// Giving the relevent divided differences, compute the HJWENO
572  /// approxmiation of the derivative.
573  fpreal64 calculateHJWENO(fpreal64 v1, fpreal64 v2, fpreal64 v3, fpreal64 v4, fpreal64 v5) const;
574 
575  /// Calculate the derivitive along the specified axis using
576  /// an HJ WENO method at the given index. The boolean
577  /// determines whether the derivitive in the positive or negative
578  /// direction is computed.
579  fpreal64 calculateDerivative(int x, int y, int z,
580  int axis, bool positivegradient) const;
581 
582  /// Converts the boundary line information from world space to local index space.
583  /// This applies an inverse transform to the plane normal and regular transform
584  /// to the height value to put both into index space.
585  SIM_BoundaryLine boundaryLineToIndex(const SIM_BoundaryLine &worldbline) const;
586 
587  /// Convert indices to world coordinates and vice-versa. Note this uses
588  /// this field's indices which change depending on sampling.
589  bool indexToPos(int x, int y, int z, UT_Vector3 &pos) const;
590 
591  /// Converts a worldspace position into an integer index.
592  bool posToIndex(UT_Vector3 pos, int &x, int &y, int &z) const;
593 
594  bool posToIndex(UT_Vector3 wpos, UT_Vector3 &ipos) const;
595 
596  /// Converts worldspace position into an integer index + the lerp values
597  /// required to interpolate.
598  /// Lerp with (1-dx) * (x,y,z) + dx * (x+1,y,z)
599  bool posToIndex(UT_Vector3 pos, int &x, int &y, int &z,
600  fpreal &dx, fpreal &dy, fpreal &dz) const;
601  /// Convert voxel cell indices to world coordinates and vice-versa.
602  /// Returns values at cell centers. Is equivalent to the indexToPos
603  /// style functions only when sampling is CENTER.
604  bool cellIndexToPos(int x, int y, int z, UT_Vector3 &pos) const;
605  bool posToCellIndex(UT_Vector3 pos, int &x, int &y, int &z) const;
606 
607  /// Verbs that can be performed on these fields.
608 
609  /// Advect a point in space according to an array of velocity
610  /// fields.
611  static void advect(UT_Vector3 &pos,
612  const SIM_RawField *velx,
613  const SIM_RawField *vely,
614  const SIM_RawField *velz, float time,
615  const SIM_RawField *collision = 0,
616  float cfl = 1.0F);
617  static void advect(UT_Vector3 &pos,
618  sim_PointVelocity getVelocity,
619  float time, float voxelsize,
620  int jobnum = 0,
621  const SIM_RawField *collision = 0,
622  float cfl = 1.0F);
623 
624  /// Advect a point with the midpoint method.
625  static void advectMidpoint(UT_Vector3 &pos,
626  const SIM_RawField *velx,
627  const SIM_RawField *vely,
628  const SIM_RawField *velz, float time,
629  const SIM_RawField *collision = 0,
630  float cfl = 1.0F);
631  static void advectMidpoint(UT_Vector3 &pos,
632  sim_PointVelocity getVelocity,
633  float time, float voxelsize,
634  int jobnum = 0,
635  const SIM_RawField *collision = 0,
636  float cfl = 1.0F);
637 
638  /// Advect a point with TVD-RK3 method.
639  static void advectRK3(UT_Vector3 &pos,
640  const SIM_RawField *velx,
641  const SIM_RawField *vely,
642  const SIM_RawField *velz, float time,
643  const SIM_RawField *collision = 0,
644  float cfl = 1.0F);
645 
646  /// Advect a point with TVD-RK4 method.
647  static void advectRK4(UT_Vector3 &pos,
648  const SIM_RawField *velx,
649  const SIM_RawField *vely,
650  const SIM_RawField *velz, float time,
651  const SIM_RawField *collision = 0,
652  float cfl = 1.0F);
653 
654  /// Move a point to the given isooffset.
655  /// Returns false if fails to complete the move in the given time.
656  bool movePtToIso(UT_Vector3 &pos,
657  fpreal goaliso,
658  fpreal maxtime,
659  fpreal tol = 1e-4) const;
660 
661  /// Advect the source field by the given set of velocity fields,
662  /// setting this to the result. Source cannot be this.
663  /// this and source are assumed to match.
664  void advect(const SIM_RawField *source,
665  const SIM_RawField *velx,
666  const SIM_RawField *vely,
667  const SIM_RawField *velz,
668  float time,
669  const SIM_RawField *collision,
670  SIM_FieldAdvection advectmethod,
671  float cfl);
672 
673  void advectMinMax(const SIM_RawField *source,
674  SIM_RawField *minf, SIM_RawField *maxf,
675  const SIM_RawField *velx,
676  const SIM_RawField *vely,
677  const SIM_RawField *velz,
678  float time,
679  const SIM_RawField *collision,
680  SIM_FieldAdvection advectmethod,
681  float cfl);
682 
684  {
686  const SIM_RawField *velx, *vely, *velz;
687  float time;
689  float cfl;
692  };
693 
694  THREADED_METHOD7(SIM_RawField, shouldMultiThread(), advectSingle,
695  const SIM_RawField *, source,
696  const SIM_RawField *, velx,
697  const SIM_RawField *, vely,
698  const SIM_RawField *, velz,
699  float, time,
700  const SIM_RawField *, collision,
701  float, cfl
702  )
703  void advectSinglePartial(const SIM_RawField *source,
705  const SIM_RawField *vely,
706  const SIM_RawField *velz,
707  float time,
708  const SIM_RawField *collision,
709  float cfl,
710  const UT_JobInfo &info);
711 
712  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), advectSingleMinMax,
713  const sim_advectParms &, parms)
714  void advectSingleMinMaxPartial(const sim_advectParms &parms,
715  const UT_JobInfo &info);
716 
717  THREADED_METHOD7(SIM_RawField, shouldMultiThread(), advectTrace,
718  const SIM_RawField *, source,
719  const SIM_RawField *, velx,
720  const SIM_RawField *, vely,
721  const SIM_RawField *, velz,
722  float, time,
723  const SIM_RawField *, collision,
724  float, cfl
725  )
726  void advectTracePartial(const SIM_RawField *source,
727  const SIM_RawField *velx,
728  const SIM_RawField *vely,
729  const SIM_RawField *velz,
730  float time,
731  const SIM_RawField *collision,
732  float cfl,
733  const UT_JobInfo &info);
734 
735  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), advectTraceMinMax,
736  const sim_advectParms &, parms)
737  void advectTraceMinMaxPartial(const sim_advectParms &parms,
738  const UT_JobInfo &info);
739 
740  THREADED_METHOD8(SIM_RawField, shouldMultiThread(), advectMultistep,
741  const SIM_RawField *, source,
742  const SIM_RawField *, velx,
743  const SIM_RawField *, vely,
744  const SIM_RawField *, velz,
745  float, time,
746  const SIM_RawField *, collision,
747  float, cfl,
748  SIM_FieldAdvection, advectmethod
749  )
750  void advectMultistepPartial(const SIM_RawField *source,
751  const SIM_RawField *velx,
752  const SIM_RawField *vely,
753  const SIM_RawField *velz,
754  float time,
755  const SIM_RawField *collision,
756  float cfl,
757  SIM_FieldAdvection advectmethod,
758  const UT_JobInfo &info);
759 
760  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), advectMultistepMinMax,
761  const sim_advectParms &, parms)
762  void advectMultistepMinMaxPartial(const sim_advectParms &parms,
763  const UT_JobInfo &info);
764 
765  /// Advects the field using HJWENO
766  THREADED_METHOD5(SIM_RawField, shouldMultiThread(), advectHJWENO,
767  const SIM_RawField &, source,
768  const SIM_RawField *, velx,
769  const SIM_RawField *, vely,
770  const SIM_RawField *, velz,
771  fpreal, timestep
772  )
773  void advectHJWENOPartial(const SIM_RawField &source,
774  const SIM_RawField *velx,
775  const SIM_RawField *vely,
776  const SIM_RawField *velz,
777  fpreal timestep,
778  const UT_JobInfo &info);
779 
780  /// Advects the field using upwind
781  THREADED_METHOD5(SIM_RawField, shouldMultiThread(), advectUpwind,
782  const SIM_RawField &, source,
783  const SIM_RawField *, velx,
784  const SIM_RawField *, vely,
785  const SIM_RawField *, velz,
786  fpreal, timestep
787  )
788  void advectUpwindPartial(const SIM_RawField &source,
789  const SIM_RawField *velx,
790  const SIM_RawField *vely,
791  const SIM_RawField *velz,
792  fpreal timestep,
793  const UT_JobInfo &info);
794 
795 
796  THREADED_METHOD4(SIM_RawField, shouldMultiThread(),
797  buoyancy,
798  const SIM_RawField *, temperature,
799  fpreal, up,
800  fpreal, Tamb,
801  fpreal, buoyancy)
802  void buoyancyPartial(const SIM_RawField *temperature,
803  fpreal up, fpreal Tamb,
804  fpreal buoyancy,
805  const UT_JobInfo &info);
806 
807  enum MIX_NAMES {
816  NUM_MIX
817  };
818 
819  /// Performs the requires mixing.
820  static fpreal mixValues(MIX_NAMES mixtype,
821  fpreal d, fpreal s)
822  {
823  switch (mixtype)
824  {
825  case MIX_COPY:
826  d = s;
827  break;
828  case MIX_ADD:
829  d += s;
830  break;
831  case MIX_SUB:
832  d -= s;
833  break;
834  case MIX_MUL:
835  d *= s;
836  break;
837  case MIX_DIV:
838  d = SYSsafediv(d, s);
839  break;
840  case MIX_MAX:
841  d = SYSmax(d, s);
842  break;
843  case MIX_MIN:
844  d = SYSmin(d, s);
845  break;
846  case MIX_AVERAGE:
847  d = d + s;
848  d *= 0.5f;
849  break;
850  case NUM_MIX:
851  UT_ASSERT(!"Invalid mix value");
852  break;
853  }
854 
855  return d;
856  }
857 
859  {
860  float threshold, bandwidth;
861  bool extrapolate, usemaxextrapolate;
863  float d_preadd, d_premul;
864  float s_preadd, s_premul;
865  float postadd, postmul;
866  float uniformrad;
868 
869  bool useattrib;
870  const char *attribname;
871  int offset;
874  bool useaffine;
875  const char *affinename;
876  };
877 
879  const sim_particleToFieldParms &parms)
880  {
881  // srcval is our source value from the particle-defined field.
882  // dstval is original value from destination field.
883  // Apply our calculation method, and return.
884  srcval *= parms.s_premul;
885  srcval += parms.s_preadd;
886  dstval *= parms.d_premul;
887  dstval += parms.d_preadd;
888  dstval = mixValues(parms.calctype, dstval, srcval);
889 
890  dstval *= parms.postmul;
891  dstval += parms.postadd;
892  return dstval;
893 
894  }
895 
896  THREADED_METHOD3(SIM_RawField, shouldMultiThread(),
897  applyParticles,
898  const GU_Detail *, particles,
899  GEO_PointTreeGAOffset *, pttree,
900  sim_particleToFieldParms &, parms)
901  void applyParticlesPartial(
902  const GU_Detail *particles,
905  const UT_JobInfo &info);
906 
907  void accumulateParticles(const GU_Detail *particles,
909  const GA_PointGroup *ptgrp = NULL);
910 
911  THREADED_METHOD6(SIM_RawField, shouldMultiThread(), advect2,
912  const SIM_RawField *, source,
913  sim_PointVelocity, getVelocity,
914  float, time, float, voxelsize,
915  const SIM_RawField *, collision,
916  float, cfl)
917  void advect2Partial(const SIM_RawField *source,
918  sim_PointVelocity getVelocity,
919  float time, float voxelsize,
920  const SIM_RawField *collision,
921  float cfl,
922  const UT_JobInfo &info);
923 
924  /// Advect this field by the given velocity fields. Invokes
925  /// advect but handles the creation of the intermediate field.
926  void advectSelf(const SIM_RawField *velx,
927  const SIM_RawField *vely,
928  const SIM_RawField *velz,
929  float time,
930  const SIM_RawField *collision,
931  SIM_FieldAdvection advectmethod,
932  float cfl);
933  void advectSelf(sim_PointVelocity getVelocity,
934  float time, float voxelsize,
935  const SIM_RawField *collision,
936  float cfl);
937 
938  /// Like advectSelf, but also generates min/max fields.
939  void advectMinMaxSelf(SIM_RawField *minfield,
940  SIM_RawField *maxfield,
941  const SIM_RawField *velx,
942  const SIM_RawField *vely,
943  const SIM_RawField *velz,
944  float time,
945  const SIM_RawField *collision,
946  SIM_FieldAdvection advectmethod,
947  float cfl);
948  /// Solve the diffusion equation over this field according to
949  /// the given diffusion rate. One likely wants to roll the
950  /// desired timestep into the diffusion rate.
951  void diffuse(fpreal diffrate,
952  int numiter,
953  const SIM_RawField *collision=0);
954 
955  /// Raw diffuse algorithm, exposed only for external performance tests
956  static void diffuse(float *dstdata, const float *srcdata[3][3][3],
957  float b, float ivsx, float ivsy, float ivsz,
958  int tx, int ty,
959  int max_xtile, int max_ytile,
960  int max_xvox, int max_yvox);
961 
962  /// Determine all components connected according to < 0 semantic
963  void computeConnectedComponents(UT_VoxelArray<int64> &comp, int &numcomponent) const;
964 
965  /// Computes the divergence of the face-centered velocity field
966  /// and stores the result in this. this must match the velocity field.
967  THREADED_METHOD2(SIM_RawField, shouldMultiThread(), buildDivergenceFace,
968  const SIM_VectorField *, vel,
969  const SIM_RawField *, surface)
970  void buildDivergenceFacePartial(const SIM_VectorField *vel,
971  const SIM_RawField *surface,
972  const UT_JobInfo &info);
973 
975  {
979  PCG_MIC
980  };
981 
982  /// Solves the pressure field that would eliminate the given
983  /// divergence field.
984  void solvePressure(const SIM_RawField *divergence,
985  const SIM_RawField *collision,
986  int numiter = 20);
987  void solvePressurePCG(const SIM_RawField &divergence,
989  SIM_VectorField *vel,
990  const SIM_RawIndexField *comp = 0,
991  const UT_IntArray *expandable = 0,
992  const SIM_RawField *surface = 0,
993  bool variational = true,
994  bool ghostfluid = true,
995  PCG_METHOD pcgmethod = PCG_MIC);
996 
997  // For each voxel v in *this and corresponding voxel b in *B
998  // sv = sum of 6-way neighbours of v
999  // v' = (sv * sumweight + b) * weight
1000  // This is only applied to the voxels whose parity (x^y^z) matches
1001  // parity, so should be called twice to complete a single iteration.
1002  THREADED_METHOD4(SIM_RawField, shouldMultiThread(),
1003  gaussSeidelIteration,
1004  const SIM_RawField *, B,
1005  fpreal32, weight,
1006  fpreal32, sumweight,
1007  int, parity)
1008  void gaussSeidelIterationPartial(
1009  const SIM_RawField * B,
1010  fpreal32 weight,
1011  fpreal32 sumweight,
1012  int parity,
1013  const UT_JobInfo &info);
1014 
1015  // Note that in the flat case the destination field is passed
1016  // parameter, as opposed to the previous function where it is
1017  // the B field that is a parameter.
1018  THREADED_METHOD4_CONST(SIM_RawField, shouldMultiThread(),
1019  gaussSeidelIterationFlat,
1020  fpreal32 *, A,
1021  fpreal32, weight,
1022  fpreal32, sumweight,
1023  int, parity)
1024  void gaussSeidelIterationFlatPartial(
1025  fpreal32 * A,
1026  fpreal32 weight,
1027  fpreal32 sumweight,
1028  int parity,
1029  const UT_JobInfo &info) const;
1030 
1031  // Note that in the flat case the destination field is passed
1032  // parameter, as opposed to the previous function where it is
1033  // the B field that is a parameter.
1034  // This does a 4-way neighbours.
1035  THREADED_METHOD4_CONST(SIM_RawField, shouldMultiThread(),
1036  gaussSeidelIterationFlat2D,
1037  fpreal32 *, A,
1038  fpreal32, weight,
1039  fpreal32, sumweight,
1040  int, parity)
1041  void gaussSeidelIterationFlat2DPartial(
1042  fpreal32 * A,
1043  fpreal32 weight,
1044  fpreal32 sumweight,
1045  int parity,
1046  const UT_JobInfo &info) const;
1047 
1048  /// Enforces the boundary conditions on this field.
1049  /// Each axis can have its own conditions
1050  /// The boundary line is expected to be in index space
1051  void enforceBoundary(SIM_FieldBoundary collisionboundary = SIM_BOUNDARY_NONE,
1052  const SIM_RawField *collision=0,
1053  const SIM_RawField *cvalue = 0,
1054  const SIM_RawField *boundary = 0,
1055  const SIM_BoundaryLine &indexbline = SIM_BoundaryLine());
1056  void enforceCollisionBoundary(SIM_FieldBoundary boundary,
1057  const SIM_RawField *collision,
1058  const SIM_RawField *cvalue = 0);
1059  void enforceSideBoundary(int axis, int side,
1060  SIM_FieldBoundary bound,
1061  fpreal boundaryvalue,
1062  const SIM_BoundaryLine &indexbline = SIM_BoundaryLine(),
1063  const SIM_RawField *boundaryfield = 0);
1064 
1065  THREADED_METHOD3(SIM_RawField, shouldMultiThread(),
1066  enforceCollisionBoundaryInternal,
1067  SIM_FieldBoundary, boundary,
1068  const SIM_RawField *, collision,
1069  const SIM_RawField *, cvalue)
1070  void enforceCollisionBoundaryInternalPartial(
1071  SIM_FieldBoundary boundary,
1072  const SIM_RawField *collision,
1073  const SIM_RawField *cvalue,
1074  const UT_JobInfo &info);
1075 
1076  /// Enforces the boundary conditions onto a flat array
1077  /// using the dimensions of this.
1078  /// Does SAME boundaries for collision objects.
1079  /// Index field is where to copy from for each collision voxel.
1080  void enforceBoundaryFlat(fpreal32 *values,
1081  const SIM_RawIndexField *collision_lookup);
1082  THREADED_METHOD2(SIM_RawField, shouldMultiThread(),
1083  enforceCollisionBoundaryFlat,
1084  fpreal32 *, values,
1085  const SIM_RawIndexField *, collision_lookup)
1086 
1087  void enforceCollisionBoundaryFlatPartial(fpreal32 *values,
1088  const SIM_RawIndexField *collision_lookup,
1089  const UT_JobInfo &info);
1090  void enforceSideBoundaryFlat(fpreal32 *values,
1091  int axis, int side,
1092  SIM_FieldBoundary bound,
1093  fpreal boundval);
1094 
1095  /// These boundary conditions do *not* apply to reading outside
1096  /// of the valid field range. The native UT_VoxelArray boundary
1097  /// condition is used for that.
1098  /// Instead, they are used for the behaviour of enforceBoundary
1099  /// and by various places where we want to distinguish open
1100  /// velocity fields from closed.
1101 
1102  void setBoundary(int axis, int side, SIM_FieldBoundary bound)
1103  { myBoundary[axis][side] = bound; }
1104  SIM_FieldBoundary getBoundary(int axis, int side) const
1105  { return myBoundary[axis][side]; }
1106  void setBoundaryValue(int axis, int side, fpreal v)
1107  { myBoundaryValue[axis][side] = v; }
1108  fpreal getBoundaryValue(int axis, int side) const
1109  { return myBoundaryValue[axis][side]; }
1110 
1111  /// Mark this field as being an extrapolated field. Out of bound
1112  /// voxels will read the clamped value. The difference between
1113  /// the clamped position and the real position is then dot producted
1114  /// with the given scale factor and added to the resulting value.
1115  /// This allows you to have rest fields that extrapolate meaningfully.
1116  void setAsExtrapolatedField(UT_Vector3 scale);
1117 
1118  /// These adjust the native UT_VoxelArray border values that *are*
1119  /// used for reading outside the valid range.
1120  UT_VoxelBorderType getBorder() const;
1121  float getBorderValue() const;
1122  void setBorder(UT_VoxelBorderType border, float bval);
1123 
1124  /// These adjust the native UT_VoxelArray comrpression options.
1125  void setCompressionOptions(const UT_VoxelCompressOptions &options);
1126  const UT_VoxelCompressOptions &getCompressionOptions() const;
1127 
1128  /// These adjust the native UT_VoxelArray comrpression tolerance..
1129  void setCompressionTolerance(fpreal tol);
1130  fpreal getCompressionTolerance() const;
1131 
1132  /// Transform turns this into a packed set of wavelet coeffecients
1133  /// from the scalar data in field. Inverse unpacks and generates
1134  /// a scalar field into this.
1135  ///
1136  /// Using a smooth edge conditon we assume that
1137  /// in case of an odd sized field, the ultimate row has
1138  /// zero diff coefficients. We thus store the extra
1139  /// column in the averages side of the matrix, for
1140  /// ceil(n/2) in the averages and floor(n/2) in the diffs.
1141  /// Note this causes 3d fields that are actually 2d
1142  /// to properly decompose.
1143  void waveletTransform(UT_Wavelet::WAVELET_NAMES wavelettype,
1144  const SIM_RawField *field,
1145  int maxpasses = -1);
1146  void waveletInverseTransform(
1147  UT_Wavelet::WAVELET_NAMES wavelettype,
1148  const SIM_RawField *wavelet,
1149  int maxpasses = -1);
1150 
1151  /// Computes the sum of squares of the given level's detail vector.
1152  void waveletComputePSD(const SIM_RawField *wavelet,
1153  int level);
1154 
1155  /// Extracts the given component from a packed wavelet array.
1156  void waveletExtractComponent(const SIM_RawField *wavelet,
1157  int level,
1158  int component);
1159 
1160  /// Builds from either a GEO_PrimVolume or a GEO_PrimVDB.
1161  /// The provided volidx is the index inside of a vector volume
1162  /// to use for the building.
1163  void buildFromPrim(const GEO_Primitive *vol,
1164  int volidx,
1165  const UT_DMatrix4 &xform,
1166  fpreal scale);
1167 
1168  /// Builds from a GEO_PrimVolume
1169  void buildFromGeo(const GEO_PrimVolume *vol,
1170  const UT_DMatrix4 &xform,
1171  fpreal scale);
1172 
1173  /// Compute fractional volume weights representing the amount the voxel
1174  /// surrounding each sample point is inside the provided SDF. This method
1175  /// just subsamples the voxel space according to the samplesperaxis
1176  /// parameter.
1177 
1178 
1180  int samplesperaxis,
1181  bool invert,
1182  fpreal minweight,
1183  fpreal dilatedist = 0)
1184  {
1185  computeSDFWeightsSampledInternal(sdf,samplesperaxis,invert,
1186  minweight, dilatedist);
1187  }
1188 
1189 private:
1190  THREADED_METHOD5(SIM_RawField, shouldMultiThread(),
1191  computeSDFWeightsSampledInternal,
1192  const SIM_RawField *, sdf,
1193  int, samplesperaxis,
1194  bool, invert,
1195  fpreal, dilatedist,
1196  fpreal, minweight);
1197 
1198  void computeSDFWeightsSampledInternalPartial(const SIM_RawField *sdf,
1199  int samplesperaxis,
1200  bool invert,
1201  fpreal minweight,
1202  fpreal dilatedist,
1203  const UT_JobInfo &info);
1204 public:
1205  /// Compute fractional volume weights representing the amount the voxel
1206  /// surrounding each sample point is inside the provided SDF. This method
1207  /// uses the normalized area of the square that is in the plane given by
1208  /// the axis parameter as an approximation.
1209  THREADED_METHOD4(SIM_RawField, shouldMultiThread(),
1210  computeSDFWeightsFace,
1211  const SIM_RawField *, sdf,
1212  int, axis,
1213  bool, invert,
1214  fpreal, minweight);
1215 
1216  void computeSDFWeightsFacePartial(const SIM_RawField *sdf,
1217  int axis,
1218  bool invert,
1219  fpreal minweight,
1220  const UT_JobInfo &info);
1221 
1222  /// Compute fractional volume weights representing the amount the voxel
1223  /// surrounding each sample point is inside the provided SDF. This method
1224  /// uses accurate volume fractions.
1225  THREADED_METHOD3(SIM_RawField, shouldMultiThread(),
1226  computeSDFWeightsVolumeFraction,
1227  const SIM_RawField *, sdf,
1228  bool, invert,
1229  fpreal, minweight);
1230 
1231  void computeSDFWeightsVolumeFractionPartial(const SIM_RawField *sdf,
1232  bool invert,
1233  fpreal minweight,
1234  const UT_JobInfo &info);
1235 
1236 
1237  /// Compute fractional volume weights for the voxel specified by x, y, z.
1238  /// This method uses the normalized area of the square that is in the
1239  /// plane given by the axis parameter as an approximation.
1240  fpreal computeVoxelSDFWeightFace(int x, int y, int z,
1241  const SIM_RawField *sdf,
1242  int axis) const;
1243 
1244  /// Scale this field by s * B / C or set to zero if scaled value is lower than
1245  /// given threshold
1246  THREADED_METHOD4(SIM_RawField, shouldMultiThread(),
1247  setScaleDivideThreshold,
1248  fpreal, scale,
1249  const SIM_RawField *, B,
1250  const SIM_RawField *, C,
1251  fpreal, threshold)
1252  void setScaleDivideThresholdPartial(fpreal scale,
1254  const SIM_RawField *C,
1255  fpreal threshold,
1256  const UT_JobInfo &jobinfo);
1257 
1258 protected:
1259  void waveletRebuildFromVoxelArray(
1260  UT_VoxelArrayF *array,
1261  float scale);
1262 
1263 
1264  static fpreal applyBoundary(SIM_FieldBoundary bound, fpreal v,
1265  fpreal boundval);
1266 
1267  THREADED_METHOD2_CONST(SIM_RawField, shouldMultiThread(), reduceOpInternal,
1268  fpreal64 *, sum,
1269  REDUCE_NAMES, op)
1270  void reduceOpInternalPartial(fpreal64 *sum, REDUCE_NAMES op, const UT_JobInfo &info) const;
1271 
1272  THREADED_METHOD5_CONST(SIM_RawField, shouldMultiThread(), reduceMaskedOpInternal,
1273  fpreal64 *, sum,
1274  fpreal64 *, masktotal,
1275  REDUCE_NAMES, op,
1276  const SIM_RawField *, mask,
1277  bool, maskissdf)
1278  void reduceMaskedOpInternalPartial(fpreal64 *sum, fpreal64 *masktotal, REDUCE_NAMES op, const SIM_RawField *maskfield, bool maskissdf, const UT_JobInfo &info) const;
1279 
1280  void sortAllVoxels(UT_FloatArray &voxelvalues) const;
1281  void sortAllVoxelsMasked(UT_FloatArray &voxelvalues,
1282  const SIM_RawField *maskfield,
1283  bool maskissdf) const;
1284 
1285  /// Methods for extrapolation
1286  THREADED_METHOD8(SIM_RawField, shouldMultiThread(), buildExtrapList,
1287  UT_Array<UT_ValArray<gu_sdf_qelem *> > &, lists,
1288  const SIM_RawField *, depths,
1289  const SIM_RawField *, valid,
1290  fpreal, isocontour,
1291  fpreal, dir,
1292  fpreal, maxdist,
1293  bool, clamp,
1294  fpreal, clampval)
1295  void buildExtrapListPartial(
1296  UT_Array<UT_ValArray<gu_sdf_qelem *> > &lists,
1297  const SIM_RawField *depths,
1298  const SIM_RawField *valid,
1299  fpreal isocontour,
1300  fpreal dir, fpreal maxdist,
1301  bool clamp, fpreal clampval,
1302  const UT_JobInfo &info);
1303 
1304  /// Reduction by reducing each axis in turn.
1305  /// This *will* change field and dst != field is required.
1306  void localReduceByAxis(REDUCE_NAMES op,
1308  UT_VoxelArrayF &field,
1309  UT_Vector3 radius) const;
1310 
1311  /// Reduce along a given axis by the specified radius in voxels.
1312  template <int AXIS, REDUCE_NAMES OP>
1313  void localReduceAxis(UT_VoxelArrayF &dstfield,
1314  const UT_VoxelArrayF &field,
1315  float radius,
1316  const UT_JobInfo &info) const;
1317  template <int AXIS, REDUCE_NAMES OP>
1318  void localMinMaxAxis(UT_VoxelArrayF &dstfield,
1319  const UT_VoxelArrayF &field,
1320  float radius,
1321  const UT_JobInfo &info) const;
1322 
1323  template <int AXIS>
1324  void localReduceAxisOp(REDUCE_NAMES op,
1325  UT_VoxelArrayF &dstfield,
1326  UT_VoxelArrayF &field,
1327  float radius,
1328  const UT_JobInfo &info) const;
1329  /// Again a triple specialization to engage threading.
1330  THREADED_METHOD4_CONST(SIM_RawField,
1331  field.getTileRes(1)*field.getTileRes(2) > 1,
1332  localReduceAxisX,
1333  REDUCE_NAMES, op,
1334  UT_VoxelArrayF &, dst,
1335  UT_VoxelArrayF &, field,
1336  float, radius)
1337  void localReduceAxisXPartial(
1338  REDUCE_NAMES op,
1339  UT_VoxelArrayF &dst,
1340  UT_VoxelArrayF &field,
1341  float radius,
1342  const UT_JobInfo &info) const
1343  { localReduceAxisOp<0>(op, dst, field, radius, info); }
1344 
1346  field.getTileRes(0)*field.getTileRes(2) > 1,
1347  localReduceAxisY,
1348  REDUCE_NAMES, op,
1349  UT_VoxelArrayF &, dst,
1350  UT_VoxelArrayF &, field,
1351  float, radius)
1352  void localReduceAxisYPartial(
1353  REDUCE_NAMES op,
1355  UT_VoxelArrayF &field,
1356  float radius,
1357  const UT_JobInfo &info) const
1358  { localReduceAxisOp<1>(op, dst, field, radius, info); }
1359 
1361  field.getTileRes(0)*field.getTileRes(1) > 1,
1362  localReduceAxisZ,
1363  REDUCE_NAMES, op,
1364  UT_VoxelArrayF &, dst,
1365  UT_VoxelArrayF &, field,
1366  float, radius)
1367  void localReduceAxisZPartial(
1368  REDUCE_NAMES op,
1370  UT_VoxelArrayF &field,
1371  float radius,
1372  const UT_JobInfo &info) const
1373  { localReduceAxisOp<2>(op, dst, field, radius, info); }
1374 
1375  THREADED_METHOD3(SIM_RawField, shouldMultiThread(), buildFromGeoSampled,
1376  const GEO_PrimVolume *, vol,
1377  const UT_DMatrix4 &, xform,
1378  fpreal, scale)
1379  void buildFromGeoSampledPartial(const GEO_PrimVolume *vol,
1381  fpreal scale,
1382  const UT_JobInfo &info);
1383 
1384  THREADED_METHOD4(SIM_RawField, shouldMultiThread(), buildFromVDBSampled,
1385  const GEO_PrimVDB *, vdb,
1386  int, vdbidx,
1387  const UT_DMatrix4 &, xform,
1388  fpreal, scale)
1389  void buildFromVDBSampledPartial(const GEO_PrimVDB *vdb,
1390  int vdbidx,
1391  const UT_DMatrix4 &xform,
1392  fpreal scale,
1393  const UT_JobInfo &info);
1394 
1396 
1397  /// We always have myField at our current resolution.
1398  /// myGrid only exists if it matches myField. If myField
1399  /// changes, myGrid will be destroyed.
1400  /// If something updates the grid the field is out of date
1401  /// so is flagged.
1402  UT_VoxelArrayF *myField;
1403  mutable CE_Grid *myGrid;
1404  mutable bool myFieldOutOfDate;
1405 
1406  // This is the size and offset which converts our UT_VoxelArray into
1407  // world coordinates.
1408  UT_Vector3 myOrig, mySize;
1409 
1410  // This is our actual official size.
1411  UT_Vector3 myBBoxOrig, myBBoxSize;
1412 
1413  UT_Vector3 myVoxelSize;
1414  fpreal myVoxelDiameter;
1415 
1416  SIM_FieldBoundary myBoundary[3][2];
1417  fpreal myBoundaryValue[3][2];
1418 };
1419 
1420 
1421 ///
1422 /// Allows you to iterate over the voxel cells of a RawField. This
1423 /// is different than iterating over the voxel samples!
1424 /// The idea of getting or setting the raw value doesn't apply
1425 /// here as we are not necessarily alligned with the sample grid.
1426 ///
1428 {
1429 public:
1431  {
1432  myRes[0] = myRes[1] = myRes[2] = 0;
1433  rewind();
1434  }
1436  {
1437  }
1438 
1439  void setArray(const SIM_RawField *field)
1440  {
1441  field->getVoxelRes(myRes[0], myRes[1], myRes[2]);
1442  rewind();
1443  }
1444 
1445  void rewind()
1446  {
1447  myIdx[0] = 0;
1448  myIdx[1] = 0;
1449  myIdx[2] = 0;
1450  }
1451 
1452  bool atEnd() const
1453  {
1454  return myIdx[2] >= myRes[2];
1455  }
1456 
1457  void advance()
1458  {
1459  myIdx[0]++;
1460  if (myIdx[0] >= myRes[0])
1461  {
1462  myIdx[0] = 0;
1463  myIdx[1]++;
1464  if (myIdx[1] >= myRes[1])
1465  {
1466  myIdx[1] = 0;
1467  myIdx[2]++;
1468  if (myIdx[2] >= myRes[2])
1469  {
1470  // All done!
1471  myIdx[2] = myRes[2];
1472  }
1473  }
1474  }
1475  }
1476 
1477  int x() const { return myIdx[0]; }
1478  int y() const { return myIdx[1]; }
1479  int z() const { return myIdx[2]; }
1480  int idx(int axis) const { return myIdx[axis]; }
1481 
1482  /// Returns true if we are at the start of a new tile.
1483  /// Since we don't work on tile basis, we define a single
1484  /// x pass to be a tile.
1485  bool isStartOfTile() const
1486  {
1487  return myIdx[0] == 0;
1488  }
1489 
1490 protected:
1491  int myIdx[3];
1492  int myRes[3];
1493 };
1494 
1495 //
1496 // These macros create the CALL_VOXELPROBE macro
1497 //
1498 // CALL_VOXELPROBE(const SIM_RawField *src,
1499 // SIM_RawField *dst,
1500 // float defaultvalue,
1501 // callFunc(dst, theprobe))
1502 //
1503 // callFunc should be
1504 //
1505 // template <typename P>
1506 // void
1507 // callFunc(SIM_RawField *dst, P &probe)
1508 //
1509 // dst->isMatching(src) must be true.
1510 //
1511 // The callFunc will receive in the parameter "theprobe" a
1512 // UT_VoxelProbeAverage properly templated against the difference
1513 // between the source and dest fields.
1514 // theprobe.setIndex(x, y, z);
1515 // will cause theprobe.getValue() to return the interpolated value
1516 // in source of the destination's voxel sample (x, y, z). This takes
1517 // into account the different sampling patterns.
1518 //
1519 #define CALL_PROBEXY(XSTEP, YSTEP, src, dst, command) \
1520 { \
1521  if ((dstsample ^ srcsample) & 4) \
1522  { \
1523  if (dstsample & 4) /* ZStep -1 */ \
1524  { \
1525  UT_VoxelProbeAverage<fpreal32, XSTEP, YSTEP, -1> theprobe; \
1526  theprobe.setArray(src->field()); \
1527  command; \
1528  } \
1529  else /* ZStep 1 */ \
1530  { \
1531  UT_VoxelProbeAverage<fpreal32, XSTEP, YSTEP, 1> theprobe; \
1532  theprobe.setArray(src->field()); \
1533  command; \
1534  } \
1535  } \
1536  else /* ZStep 0 */ \
1537  { \
1538  UT_VoxelProbeAverage<fpreal32, XSTEP, YSTEP, 0> theprobe; \
1539  theprobe.setArray(src->field()); \
1540  command; \
1541  } \
1542 }
1543 
1544 #define CALL_PROBEX(XSTEP, src, dst, command) \
1545 { \
1546  if ((dstsample ^ srcsample) & 2) \
1547  { \
1548  if (dstsample & 2) /* YStep -1 */ \
1549  { \
1550  CALL_PROBEXY(XSTEP, -1, src, dst, command); \
1551  } \
1552  else /* YStep 1 */ \
1553  { \
1554  CALL_PROBEXY(XSTEP, 1, src, dst, command); \
1555  } \
1556  } \
1557  else /* YStep 0 */ \
1558  { \
1559  CALL_PROBEXY(XSTEP, 0, src, dst, command); \
1560  } \
1561 }
1562 
1563 #define CALL_VOXELPROBE(src, dst, defval, command) \
1564 { \
1565  if (!src) \
1566  { \
1567  UT_VoxelProbeConstant<fpreal32> theprobe; \
1568  theprobe.setValue(defval); \
1569  command; \
1570  } \
1571  else \
1572  { \
1573  int dstsample, srcsample; \
1574  dstsample = dst->getSample(); \
1575  srcsample = src->getSample(); \
1576  if ((dstsample ^ srcsample) & 1) \
1577  { \
1578  if (dstsample & 1) /* XStep -1 */ \
1579  { \
1580  CALL_PROBEX(-1, src, dst, command); \
1581  } \
1582  else /* XStep 1 */ \
1583  { \
1584  CALL_PROBEX(1, src, dst, command); \
1585  } \
1586  } \
1587  else /* XStep 0 */ \
1588  { \
1589  CALL_PROBEX(0, src, dst, command); \
1590  } \
1591  } \
1592 }
1593 
1594 #endif
#define SYSmax(a, b)
Definition: SYS_Math.h:1365
#define THREADED_METHOD7(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3, PARMTYPE4, PARMNAME4, PARMTYPE5, PARMNAME5, PARMTYPE6, PARMNAME6, PARMTYPE7, PARMNAME7)
#define THREADED_METHOD4_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3, PARMTYPE4, PARMNAME4)
virtual ~SIM_RawFieldCellIterator()
void getVoxelRes(int &xres, int &yres, int &zres) const
#define THREADED_METHOD8(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3, PARMTYPE4, PARMNAME4, PARMTYPE5, PARMNAME5, PARMTYPE6, PARMNAME6, PARMTYPE7, PARMNAME7, PARMTYPE8, PARMNAME8)
GLboolean invert
Definition: glcorearb.h:548
const UT_Vector3 & getVoxelSize() const
Definition: SIM_RawField.h:496
#define THREADED_METHOD1(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1)
static fpreal applyParticleToFieldParms(fpreal srcval, fpreal dstval, const sim_particleToFieldParms &parms)
Definition: SIM_RawField.h:878
const GLdouble * v
Definition: glcorearb.h:836
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
exint getXRes() const
Definition: SIM_RawField.h:515
UT_Vector3T< float > UT_Vector3
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:847
static fpreal mixValues(MIX_NAMES mixtype, fpreal d, fpreal s)
Performs the requires mixing.
Definition: SIM_RawField.h:820
GLint level
Definition: glcorearb.h:107
UT_VoxelArrayF * fieldNC() const
Definition: SIM_RawField.h:470
UT_VoxelBorderType
Definition: UT_VoxelArray.h:69
GLint GLuint mask
Definition: glcorearb.h:123
const SIM_RawField * velz
Definition: SIM_RawField.h:686
GLint y
Definition: glcorearb.h:102
GLfloat GLfloat GLfloat v2
Definition: glcorearb.h:817
GLfloat GLfloat GLfloat GLfloat v3
Definition: glcorearb.h:818
bool shouldMultiThread() const
Returns true if this should be multithreaded.
Definition: SIM_RawField.h:130
#define THREADED_METHOD5_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3, PARMTYPE4, PARMNAME4, PARMTYPE5, PARMNAME5)
int idx(int axis) const
GLsizeiptr size
Definition: glcorearb.h:663
#define THREADED_METHOD3(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3)
SIM_FieldSample
Definition: SIM_RawField.h:38
UT_Vector3(* sim_PointVelocity)(const UT_Vector3 &, int)
Definition: SIM_RawField.h:95
long long int64
Definition: SYS_Types.h:100
GA_API const UT_StringHolder scale
fpreal getVoxelDiameter() const
Definition: SIM_RawField.h:497
void markGridAsChanged()
Mark the field as out of date, but only if we have a valid grid.
Definition: SIM_RawField.h:479
SIM_FieldSample getSample() const
Definition: SIM_RawField.h:503
exint getZRes() const
Definition: SIM_RawField.h:517
const UT_Vector3 & getBBoxOrig() const
Definition: SIM_RawField.h:494
exint getYRes() const
Definition: SIM_RawField.h:516
#define UT_ASSERT(ZZ)
Definition: UT_Assert.h:102
int64 exint
Definition: SYS_Types.h:109
const UT_Vector3 & getBBoxSize() const
Definition: SIM_RawField.h:495
double fpreal64
Definition: SYS_Types.h:185
#define THREADED_METHOD2(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2)
void computeSDFWeightsSampled(const SIM_RawField *sdf, int samplesperaxis, bool invert, fpreal minweight, fpreal dilatedist=0)
bool isStartOfTile() const
GLsizei GLsizei GLchar * source
Definition: glcorearb.h:802
const SIM_RawField * collision
Definition: SIM_RawField.h:688
void setBoundaryValue(int axis, int side, fpreal v)
const UT_Vector3 & getSize() const
Definition: SIM_RawField.h:493
UT_Vector3 planeNormal
Definition: SIM_RawField.h:90
const SIM_RawField * source
Definition: SIM_RawField.h:685
GLenum GLsizei GLsizei GLint * values
Definition: glcorearb.h:1601
#define THREADED_METHOD5(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3, PARMTYPE4, PARMNAME4, PARMTYPE5, PARMNAME5)
fpreal average() const
Prebuilt versions.
Definition: SIM_RawField.h:420
#define THREADED_METHOD4(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3, PARMTYPE4, PARMNAME4)
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
void setArray(const SIM_RawField *field)
fpreal getBoundaryValue(int axis, int side) const
const UT_VoxelArrayF * field() const
Fetches the raw field.
Definition: SIM_RawField.h:469
#define THREADED_METHOD2_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2)
GLint GLsizei GLsizei height
Definition: glcorearb.h:102
GLenum GLenum dst
Definition: glcorearb.h:1792
GLsizei const GLfloat * value
Definition: glcorearb.h:823
double fpreal
Definition: SYS_Types.h:263
GLint GLint GLsizei GLint border
Definition: glcorearb.h:107
typedef int
Definition: png.h:1175
GA_API const UT_StringHolder up
GLuint index
Definition: glcorearb.h:785
SIM_FieldAdvection
Definition: SIM_RawField.h:60
void boxBlur(float radius)
Definition: SIM_RawField.h:429
#define SIM_API
Definition: SIM_API.h:10
GLint GLenum GLint x
Definition: glcorearb.h:408
#define THREADED_METHOD6(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3, PARMTYPE4, PARMNAME4, PARMTYPE5, PARMNAME5, PARMTYPE6, PARMNAME6)
GLfloat GLfloat v1
Definition: glcorearb.h:816
GLuint GLfloat * val
Definition: glcorearb.h:1607
sim_extrapolateCacheElem(fpreal v, fpreal w, const UT_Vector3 &p)
Definition: SIM_RawField.h:199
GLubyte GLubyte GLubyte GLubyte w
Definition: glcorearb.h:856
fpreal getVoxelVolume() const
Definition: SIM_RawField.h:498
bool hasNan() const
Returns true if our field has any NANs.
Definition: SIM_RawField.h:462
VectorToScalarConverter< GridType >::Type::Ptr divergence(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the divergence of the given vector-valued grid.
SIM_FieldAdvection advectmethod
Definition: SIM_RawField.h:691
GLboolean r
Definition: glcorearb.h:1221
#define const
Definition: zconf.h:214
REDUCE_NAMES
Types of reductions supported by reduceOp.
Definition: SIM_RawField.h:400
const UT_Vector3 & getOrig() const
Definition: SIM_RawField.h:492
#define THREADED_METHOD3_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3)
#define SYSmin(a, b)
Definition: SYS_Math.h:1366
This class holds a three dimensional vector field.
float fpreal32
Definition: SYS_Types.h:184
#define THREADED_METHOD(CLASSNAME, DOMULTI, METHOD)
SIM_FieldBoundary getBoundary(int axis, int side) const
GLenum clamp
Definition: glcorearb.h:1233
GLenum src
Definition: glcorearb.h:1792
SIM_FieldBoundary
Definition: SIM_RawField.h:50