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