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