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