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 two fields have voxels aligned.
487  /// They do not have to be of the same resolution or share origin.
488  /// However, using the `offset` allows to access this field with index
489  /// of the other field.
490  /// i.e. field->indexToPos(index) == this->indexToPos(index + offset)
491  /// and this->indexToPos(index) == field->indexToPos(index - offset)
492  bool isColocated(const SIM_RawField *field, UT_Vector3I &offset) const;
493  bool isColocated(const SIM_RawIndexField *field, UT_Vector3I &offset) const;
494  bool isColocated(const GEO_Primitive *prim, UT_Vector3I &offset) const;
495 
496  /// Returns true if our field has any NANs
497  bool hasNan() const { return field()->hasNan(); }
498 
499  /// Tests for nans, outputs text and asserts if has any
500  /// Only runs if the test environment variable is set.
501  void testForNan() const;
502 
503  /// Fetches the raw field.
504  const UT_VoxelArrayF *field() const { if (myFieldOutOfDate) { updateFieldFromGrid(); } return myField; }
505  UT_VoxelArrayF *fieldNC() const { if (myGrid) { clearGrid(); } return myField; }
506 
507  /// Fetches the compute grid, returns 0 if not present.
508  CE_Grid *grid() const;
509 
510  /// Fetches the compute grid but throws cl::Error() on failure.
511  CE_Grid *requireGrid() const;
512 
513  /// Mark the field as out of date, but only if we have a valid grid.
514  void markGridAsChanged() {if (myGrid) myFieldOutOfDate = true;}
515  void clearGrid() const;
516 
517  /// Returns true if the OpenCL grid has more up-to-date data than the voxel array.
518  bool gridIsAuthoritative() const { return myGrid && myFieldOutOfDate; }
519 
520  /// Copies our field into the grid.
521  void updateGridFromField() const;
522 
523  /// Copies our grid into our field.
524  void updateFieldFromGrid() const;
525 
526  /// Steals the voxel array, leaving this pointing to a 0 constant
527  /// array
528  UT_VoxelArrayF *steal();
529  /// Adopts the heap-allocated voxel array. If the new voxel array is not
530  /// of the right size, this field is left untouched, and the new voxel array
531  /// is deleted.
532  void adopt(UT_VoxelArrayF* newf);
533 
534  const UT_Vector3 &getOrig() const { return myOrig; }
535  const UT_Vector3 &getSize() const { return mySize; }
536  const UT_Vector3 &getBBoxOrig() const { return myBBoxOrig; }
537  const UT_Vector3 &getBBoxSize() const { return myBBoxSize; }
538  const UT_Vector3 &getVoxelSize() const { return myVoxelSize; }
539  void setVoxelSize(const UT_Vector3 &voxelsize)
540  { myVoxelSize = voxelsize; myVoxelDiameter = voxelsize.length(); }
541  fpreal getVoxelDiameter() const { return myVoxelDiameter; }
543  { return myVoxelSize.x()
544  * myVoxelSize.y()
545  * myVoxelSize.z(); }
546 
547  SIM_FieldSample getSample() const { return mySample; }
548 
549  int64 getMemoryUsage() const;
550 
551  /// Returns the resolution of the voxel grid that we are sampling.
552  /// This is a count of voxels, so may differ for our different
553  /// sampling methods.
554  void getVoxelRes(int &xres, int &yres, int &zres) const;
555  UT_Vector3I getVoxelRes() const;
556 
557  /// Returns the actual number of samples in each resolution.
558  /// Preferred over field()->getXRes() as it doesn't require a
559  /// copy of the CE fields.
560  exint getXRes() const { return myField->getXRes(); }
561  exint getYRes() const { return myField->getYRes(); }
562  exint getZRes() const { return myField->getZRes(); }
563 
564  /// Functions to resolve quantities about the field.
565 
566  /// Sample the field at the given world space location.
567  fpreal getValue(UT_Vector3 pos) const;
568 
569  /// Returns the field value at the given voxel.
570  fpreal getValueIndex(UT_Vector3I index) const;
571 
572  /// Returns an averaged value for the center of the given voxel.
573  fpreal getCellValue(int x, int y, int z) const;
574 
575  /// Adds the value to the given voxel cell, averaging out among
576  /// adjacent samples if we aren't sampled evenly.
577  void addToCell(int x, int y, int z, fpreal v);
578 
579  /// Ensures the given voxel cell has the given value. This
580  /// will set all of the adjacent samples if we aren't sampled
581  /// evenly.
582  void setCellValue(int x, int y, int z, fpreal v);
583 
584  /// Return the gradient of the field at the given world space location.
585  /// Uses central differencing with a sample spacing of the voxelsize.
586  UT_Vector3 getGradient(UT_Vector3 pos) const;
587 
588  /// Returns the gradient at the given voxel index.
589  /// Uses central differencing.
590  UT_Vector3 getGradientAtIndex(int x, int y, int z) const;
591 
592  /// Computes the laplacian of the field at the world space
593  /// coordinate using interpolation of getLaplacianAtIndex
594  fpreal64 getLaplacian(UT_Vector3 pos) const;
595 
596  /// Computes the laplacian of the field at the specific voxel index
597  /// specified
598  fpreal64 getLaplacianAtIndex(int x, int y, int z) const;
599 
600  /// Returns the curvature of the field at a given world space location.
601  /// Uses interpolation of index based getCurvatureAtIndex.
602  fpreal64 getCurvature(UT_Vector3 pos) const;
603 
604  /// Computes the curvature at the given voxel index.
605  /// Uses central differencing.
606  /// The resulting curvature is clamped according to the resolution
607  /// of the field to avoid abnormally large values in noisy data.
608  fpreal64 getCurvatureAtIndex(int x, int y, int z) const;
609 
610  /// Computes the curvature according to a 3^3 voxel probe
611  static fpreal64 getCurvatureAtProbe(UT_VoxelROProbeCubeF &probe, const UT_Vector3 &invvoxelsize);
612 
613  /// Computes K |grad(phi|), the curvature times the length of the
614  /// gradient. By folding the two operations, we can return a non-zero
615  /// value where grad goes to zero (but curvature doesn't) using
616  /// L'Hospital's rule. This also does not clamp the curvature.
617  fpreal64 getCurvatureTimesGradAtIndex(int x, int y, int z) const;
618 
619  /// Giving the relevent divided differences, compute the HJWENO
620  /// approxmiation of the derivative.
621  fpreal64 calculateHJWENO(fpreal64 v1, fpreal64 v2, fpreal64 v3, fpreal64 v4, fpreal64 v5) const;
622 
623  /// Calculate the derivitive along the specified axis using
624  /// an HJ WENO method at the given index. The boolean
625  /// determines whether the derivitive in the positive or negative
626  /// direction is computed.
627  fpreal64 calculateDerivative(int x, int y, int z,
628  int axis, bool positivegradient) const;
629 
630  /// Converts the boundary line information from world space to local index space.
631  /// This applies an inverse transform to the plane normal and regular transform
632  /// to the height value to put both into index space.
633  SIM_BoundaryLine boundaryLineToIndex(const SIM_BoundaryLine &worldbline) const;
634 
635  /// Convert indices to world coordinates and vice-versa. Note this uses
636  /// this field's indices which change depending on sampling.
637  bool indexToPos(int x, int y, int z, UT_Vector3 &pos) const;
638  bool indexToPos(exint x, exint y, exint z, UT_Vector3D &pos) const;
639  bool indexToPos(UT_Vector3I index, UT_Vector3 &pos) const;
640  UT_Vector3 indexToPos(UT_Vector3I index) const;
641 
642 
643  /// Converts a worldspace position into an integer index.
644  bool posToIndex(UT_Vector3 pos, int &x, int &y, int &z) const;
645 
646  bool posToIndex(UT_Vector3 wpos, UT_Vector3 &ipos) const;
647  UT_Vector3I posToIndex(UT_Vector3 pos) const;
648 
649  /// Converts worldspace position into an integer index + the lerp values
650  /// required to interpolate.
651  /// Lerp with (1-dx) * (x,y,z) + dx * (x+1,y,z)
652  bool posToIndex(UT_Vector3 pos, int &x, int &y, int &z,
653  fpreal &dx, fpreal &dy, fpreal &dz) const;
654  /// Convert voxel cell indices to world coordinates and vice-versa.
655  /// Returns values at cell centers. Is equivalent to the indexToPos
656  /// style functions only when sampling is CENTER.
657  bool cellIndexToPos(int x, int y, int z, UT_Vector3 &pos) const;
658  bool posToCellIndex(UT_Vector3 pos, int &x, int &y, int &z) const;
659 
660  /// Verbs that can be performed on these fields.
661 
662  /// Advect a point in space according to an array of velocity
663  /// fields.
664  static void advect(UT_Vector3 &pos,
665  const SIM_RawField *velx,
666  const SIM_RawField *vely,
667  const SIM_RawField *velz, float time,
668  const SIM_RawField *collision = 0,
669  float cfl = 1.0F);
670  static void advect(UT_Vector3 &pos,
671  sim_PointVelocity getVelocity,
672  float time, float voxelsize,
673  int jobnum = 0,
674  const SIM_RawField *collision = 0,
675  float cfl = 1.0F);
676 
677  /// Advect a point with the midpoint method.
678  static void advectMidpoint(UT_Vector3 &pos,
679  const SIM_RawField *velx,
680  const SIM_RawField *vely,
681  const SIM_RawField *velz, float time,
682  const SIM_RawField *collision = 0,
683  float cfl = 1.0F);
684  static void advectMidpoint(UT_Vector3 &pos,
685  sim_PointVelocity getVelocity,
686  float time, float voxelsize,
687  int jobnum = 0,
688  const SIM_RawField *collision = 0,
689  float cfl = 1.0F);
690 
691  /// Advect a point with TVD-RK3 method.
692  static void advectRK3(UT_Vector3 &pos,
693  const SIM_RawField *velx,
694  const SIM_RawField *vely,
695  const SIM_RawField *velz, float time,
696  const SIM_RawField *collision = 0,
697  float cfl = 1.0F);
698 
699  /// Advect a point with TVD-RK4 method.
700  static void advectRK4(UT_Vector3 &pos,
701  const SIM_RawField *velx,
702  const SIM_RawField *vely,
703  const SIM_RawField *velz, float time,
704  const SIM_RawField *collision = 0,
705  float cfl = 1.0F);
706 
707  /// Move a point to the given isooffset.
708  /// Returns false if fails to complete the move in the given time.
709  bool movePtToIso(UT_Vector3 &pos,
710  fpreal goaliso,
711  fpreal maxtime,
712  fpreal tol = 1e-4) const;
713 
714  /// Advect the source field by the given set of velocity fields,
715  /// setting this to the result. Source cannot be this.
716  /// this and source are assumed to match.
717  void advect(const SIM_RawField *source,
718  const SIM_RawField *velx,
719  const SIM_RawField *vely,
720  const SIM_RawField *velz,
721  float time,
722  const SIM_RawField *collision,
723  SIM_FieldAdvection advectmethod,
724  float cfl);
725 
726  void advectMinMax(const SIM_RawField *source,
727  SIM_RawField *minf, SIM_RawField *maxf,
728  const SIM_RawField *velx,
729  const SIM_RawField *vely,
730  const SIM_RawField *velz,
731  float time,
732  const SIM_RawField *collision,
733  SIM_FieldAdvection advectmethod,
734  float cfl);
735 
737  {
739  const SIM_RawField *velx, *vely, *velz;
740  float time;
742  float cfl;
745  };
746 
747  THREADED_METHOD7(SIM_RawField, shouldMultiThread(), advectSingle,
748  const SIM_RawField *, source,
749  const SIM_RawField *, velx,
750  const SIM_RawField *, vely,
751  const SIM_RawField *, velz,
752  float, time,
753  const SIM_RawField *, collision,
754  float, cfl
755  )
756  void advectSinglePartial(const SIM_RawField *source,
758  const SIM_RawField *vely,
759  const SIM_RawField *velz,
760  float time,
761  const SIM_RawField *collision,
762  float cfl,
763  const UT_JobInfo &info);
764 
765  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), advectSingleMinMax,
766  const sim_advectParms &, parms)
767  void advectSingleMinMaxPartial(const sim_advectParms &parms,
768  const UT_JobInfo &info);
769 
770  THREADED_METHOD7(SIM_RawField, shouldMultiThread(), advectTrace,
771  const SIM_RawField *, source,
772  const SIM_RawField *, velx,
773  const SIM_RawField *, vely,
774  const SIM_RawField *, velz,
775  float, time,
776  const SIM_RawField *, collision,
777  float, cfl
778  )
779  void advectTracePartial(const SIM_RawField *source,
780  const SIM_RawField *velx,
781  const SIM_RawField *vely,
782  const SIM_RawField *velz,
783  float time,
784  const SIM_RawField *collision,
785  float cfl,
786  const UT_JobInfo &info);
787 
788  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), advectTraceMinMax,
789  const sim_advectParms &, parms)
790  void advectTraceMinMaxPartial(const sim_advectParms &parms,
791  const UT_JobInfo &info);
792 
793  THREADED_METHOD8(SIM_RawField, shouldMultiThread(), advectMultistep,
794  const SIM_RawField *, source,
795  const SIM_RawField *, velx,
796  const SIM_RawField *, vely,
797  const SIM_RawField *, velz,
798  float, time,
799  const SIM_RawField *, collision,
800  float, cfl,
801  SIM_FieldAdvection, advectmethod
802  )
803  void advectMultistepPartial(const SIM_RawField *source,
804  const SIM_RawField *velx,
805  const SIM_RawField *vely,
806  const SIM_RawField *velz,
807  float time,
808  const SIM_RawField *collision,
809  float cfl,
810  SIM_FieldAdvection advectmethod,
811  const UT_JobInfo &info);
812 
813  THREADED_METHOD1(SIM_RawField, shouldMultiThread(), advectMultistepMinMax,
814  const sim_advectParms &, parms)
815  void advectMultistepMinMaxPartial(const sim_advectParms &parms,
816  const UT_JobInfo &info);
817 
818  /// Advects the field using HJWENO
819  THREADED_METHOD5(SIM_RawField, shouldMultiThread(), advectHJWENO,
820  const SIM_RawField &, source,
821  const SIM_RawField *, velx,
822  const SIM_RawField *, vely,
823  const SIM_RawField *, velz,
824  fpreal, timestep
825  )
826  void advectHJWENOPartial(const SIM_RawField &source,
827  const SIM_RawField *velx,
828  const SIM_RawField *vely,
829  const SIM_RawField *velz,
830  fpreal timestep,
831  const UT_JobInfo &info);
832 
833  /// Advects the field using upwind
834  THREADED_METHOD5(SIM_RawField, shouldMultiThread(), advectUpwind,
835  const SIM_RawField &, source,
836  const SIM_RawField *, velx,
837  const SIM_RawField *, vely,
838  const SIM_RawField *, velz,
839  fpreal, timestep
840  )
841  void advectUpwindPartial(const SIM_RawField &source,
842  const SIM_RawField *velx,
843  const SIM_RawField *vely,
844  const SIM_RawField *velz,
845  fpreal timestep,
846  const UT_JobInfo &info);
847 
848 
849  THREADED_METHOD5(SIM_RawField, shouldMultiThread(),
850  buoyancy,
851  const SIM_RawField *, stencil,
852  const SIM_RawField *, temperature,
853  fpreal, up,
854  fpreal, Tamb,
855  fpreal, buoyancy)
856  void buoyancyPartial(const SIM_RawField *stencil,
857  const SIM_RawField *temperature,
858  fpreal up, fpreal Tamb,
859  fpreal buoyancy,
860  const UT_JobInfo &info);
861 
862  enum MIX_NAMES {
871  NUM_MIX
872  };
873 
874  /// Performs the requires mixing.
875  static fpreal mixValues(MIX_NAMES mixtype,
876  fpreal d, fpreal s)
877  {
878  switch (mixtype)
879  {
880  case MIX_COPY:
881  d = s;
882  break;
883  case MIX_ADD:
884  d += s;
885  break;
886  case MIX_SUB:
887  d -= s;
888  break;
889  case MIX_MUL:
890  d *= s;
891  break;
892  case MIX_DIV:
893  d = SYSsafediv(d, s);
894  break;
895  case MIX_MAX:
896  d = SYSmax(d, s);
897  break;
898  case MIX_MIN:
899  d = SYSmin(d, s);
900  break;
901  case MIX_AVERAGE:
902  d = d + s;
903  d *= 0.5f;
904  break;
905  case NUM_MIX:
906  UT_ASSERT(!"Invalid mix value");
907  break;
908  }
909 
910  return d;
911  }
912 
914  {
915  float threshold, bandwidth;
916  bool extrapolate, usemaxextrapolate;
918  float d_preadd, d_premul;
919  float s_preadd, s_premul;
920  float postadd, postmul;
921  float uniformrad;
923 
924  bool useattrib;
925  const char *attribname;
926  int offset;
929  bool useaffine;
930  const char *affinename;
931  };
932 
934  const sim_particleToFieldParms &parms)
935  {
936  // srcval is our source value from the particle-defined field.
937  // dstval is original value from destination field.
938  // Apply our calculation method, and return.
939  srcval *= parms.s_premul;
940  srcval += parms.s_preadd;
941  dstval *= parms.d_premul;
942  dstval += parms.d_preadd;
943  dstval = mixValues(parms.calctype, dstval, srcval);
944 
945  dstval *= parms.postmul;
946  dstval += parms.postadd;
947  return dstval;
948 
949  }
950 
951  THREADED_METHOD3(SIM_RawField, shouldMultiThread(),
952  applyParticles,
953  const GU_Detail *, particles,
954  GEO_PointTreeGAOffset *, pttree,
955  sim_particleToFieldParms &, parms)
956  void applyParticlesPartial(
957  const GU_Detail *particles,
960  const UT_JobInfo &info);
961 
962  void accumulateParticles(const GU_Detail *particles,
964  const GA_PointGroup *ptgrp = NULL);
965 
966  THREADED_METHOD6(SIM_RawField, shouldMultiThread(), advect2,
967  const SIM_RawField *, source,
968  sim_PointVelocity, getVelocity,
969  float, time, float, voxelsize,
970  const SIM_RawField *, collision,
971  float, cfl)
972  void advect2Partial(const SIM_RawField *source,
973  sim_PointVelocity getVelocity,
974  float time, float voxelsize,
975  const SIM_RawField *collision,
976  float cfl,
977  const UT_JobInfo &info);
978 
979  /// Advect this field by the given velocity fields. Invokes
980  /// advect but handles the creation of the intermediate field.
981  void advectSelf(const SIM_RawField *velx,
982  const SIM_RawField *vely,
983  const SIM_RawField *velz,
984  float time,
985  const SIM_RawField *collision,
986  SIM_FieldAdvection advectmethod,
987  float cfl);
988  void advectSelf(sim_PointVelocity getVelocity,
989  float time, float voxelsize,
990  const SIM_RawField *collision,
991  float cfl);
992 
993  /// Like advectSelf, but also generates min/max fields.
994  void advectMinMaxSelf(SIM_RawField *minfield,
995  SIM_RawField *maxfield,
996  const SIM_RawField *velx,
997  const SIM_RawField *vely,
998  const SIM_RawField *velz,
999  float time,
1000  const SIM_RawField *collision,
1001  SIM_FieldAdvection advectmethod,
1002  float cfl);
1003  /// Solve the diffusion equation over this field according to
1004  /// the given diffusion rate. One likely wants to roll the
1005  /// desired timestep into the diffusion rate.
1006  void diffuse(fpreal diffrate,
1007  int numiter,
1008  const SIM_RawField *collision=0);
1009 
1010  /// Raw diffuse algorithm, exposed only for external performance tests
1011  static void diffuse(float *dstdata, const float *srcdata[3][3][3],
1012  float b, float ivsx, float ivsy, float ivsz,
1013  int tx, int ty,
1014  int max_xtile, int max_ytile,
1015  int max_xvox, int max_yvox);
1016 
1017  /// Determine all components connected according to < 0 semantic
1018  void computeConnectedComponents(UT_VoxelArray<int64> &comp, int &numcomponent) const;
1019 
1020  /// Computes the divergence of the face-centered velocity field
1021  /// and stores the result in this. this must match the velocity field.
1022  /// If stencil is provided (and surface is not), the divergence is only
1023  /// calculated in the areas where stencil value is greater than 0.5 (with
1024  /// the remaining regions left untouched). stencil field must be aligned
1025  /// with this.
1026  THREADED_METHOD3(SIM_RawField, shouldMultiThread(), buildDivergenceFace,
1027  const SIM_VectorField *, vel,
1028  const SIM_RawField *, surface,
1029  const SIM_RawField *, stencil)
1030  void buildDivergenceFacePartial(const SIM_VectorField *vel,
1031  const SIM_RawField *surface,
1032  const SIM_RawField *stencil,
1033  const UT_JobInfo &info);
1034 
1036  {
1040  PCG_MIC
1041  };
1042 
1043  /// Solves the pressure field that would eliminate the given
1044  /// divergence field.
1045  void solvePressure(const SIM_RawField *divergence,
1046  const SIM_RawField *collision,
1047  int numiter = 20);
1048  void solvePressurePCG(const SIM_RawField &divergence,
1050  SIM_VectorField *vel,
1051  const SIM_RawIndexField *comp = 0,
1052  const UT_IntArray *expandable = 0,
1053  const SIM_RawField *surface = 0,
1054  bool variational = true,
1055  bool ghostfluid = true,
1056  PCG_METHOD pcgmethod = PCG_MIC);
1057 
1058  // For each voxel v in *this and corresponding voxel b in *B
1059  // sv = sum of 6-way neighbours of v
1060  // v' = (sv * sumweight + b) * weight
1061  // This is only applied to the voxels whose parity (x^y^z) matches
1062  // parity, so should be called twice to complete a single iteration.
1063  THREADED_METHOD4(SIM_RawField, shouldMultiThread(),
1064  gaussSeidelIteration,
1065  const SIM_RawField *, B,
1066  fpreal32, weight,
1067  fpreal32, sumweight,
1068  int, parity)
1069  void gaussSeidelIterationPartial(
1070  const SIM_RawField * B,
1072  fpreal32 sumweight,
1073  int parity,
1074  const UT_JobInfo &info);
1075 
1076  // Note that in the flat case the destination field is passed
1077  // parameter, as opposed to the previous function where it is
1078  // the B field that is a parameter.
1079  THREADED_METHOD4_CONST(SIM_RawField, shouldMultiThread(),
1080  gaussSeidelIterationFlat,
1081  fpreal32 *, A,
1082  fpreal32, weight,
1083  fpreal32, sumweight,
1084  int, parity)
1085  void gaussSeidelIterationFlatPartial(
1086  fpreal32 * A,
1087  fpreal32 weight,
1088  fpreal32 sumweight,
1089  int parity,
1090  const UT_JobInfo &info) const;
1091 
1092  // Note that in the flat case the destination field is passed
1093  // parameter, as opposed to the previous function where it is
1094  // the B field that is a parameter.
1095  // This does a 4-way neighbours.
1096  THREADED_METHOD4_CONST(SIM_RawField, shouldMultiThread(),
1097  gaussSeidelIterationFlat2D,
1098  fpreal32 *, A,
1099  fpreal32, weight,
1100  fpreal32, sumweight,
1101  int, parity)
1102  void gaussSeidelIterationFlat2DPartial(
1103  fpreal32 * A,
1104  fpreal32 weight,
1105  fpreal32 sumweight,
1106  int parity,
1107  const UT_JobInfo &info) const;
1108 
1109  /// Enforces the boundary conditions on this field.
1110  /// Each axis can have its own conditions
1111  /// The boundary line is expected to be in index space
1112  void enforceBoundary(SIM_FieldBoundary collisionboundary = SIM_BOUNDARY_NONE,
1113  const SIM_RawField *collision=0,
1114  const SIM_RawField *cvalue = 0,
1115  const SIM_RawField *boundary = 0,
1116  const SIM_BoundaryLine &indexbline = SIM_BoundaryLine());
1117  void enforceCollisionBoundary(SIM_FieldBoundary boundary,
1118  const SIM_RawField *collision,
1119  const SIM_RawField *cvalue = 0);
1120  void enforceSideBoundary(int axis, int side,
1121  SIM_FieldBoundary bound,
1122  fpreal boundaryvalue,
1123  const SIM_BoundaryLine &indexbline = SIM_BoundaryLine(),
1124  const SIM_RawField *boundaryfield = 0);
1125 
1126  THREADED_METHOD3(SIM_RawField, shouldMultiThread(),
1127  enforceCollisionBoundaryInternal,
1128  SIM_FieldBoundary, boundary,
1129  const SIM_RawField *, collision,
1130  const SIM_RawField *, cvalue)
1131  void enforceCollisionBoundaryInternalPartial(
1132  SIM_FieldBoundary boundary,
1133  const SIM_RawField *collision,
1134  const SIM_RawField *cvalue,
1135  const UT_JobInfo &info);
1136 
1137  /// Enforces the boundary conditions onto a flat array
1138  /// using the dimensions of this.
1139  /// Does SAME boundaries for collision objects.
1140  /// Index field is where to copy from for each collision voxel.
1141  void enforceBoundaryFlat(fpreal32 *values,
1142  const SIM_RawIndexField *collision_lookup);
1143  THREADED_METHOD2(SIM_RawField, shouldMultiThread(),
1144  enforceCollisionBoundaryFlat,
1145  fpreal32 *, values,
1146  const SIM_RawIndexField *, collision_lookup)
1147 
1148  void enforceCollisionBoundaryFlatPartial(fpreal32 *values,
1149  const SIM_RawIndexField *collision_lookup,
1150  const UT_JobInfo &info);
1151  void enforceSideBoundaryFlat(fpreal32 *values,
1152  int axis, int side,
1153  SIM_FieldBoundary bound,
1154  fpreal boundval);
1155 
1156  /// These boundary conditions do *not* apply to reading outside
1157  /// of the valid field range. The native UT_VoxelArray boundary
1158  /// condition is used for that.
1159  /// Instead, they are used for the behaviour of enforceBoundary
1160  /// and by various places where we want to distinguish open
1161  /// velocity fields from closed.
1162 
1163  void setBoundary(int axis, int side, SIM_FieldBoundary bound)
1164  { myBoundary[axis][side] = bound; }
1165  SIM_FieldBoundary getBoundary(int axis, int side) const
1166  { return myBoundary[axis][side]; }
1167  void setBoundaryValue(int axis, int side, fpreal v)
1168  { myBoundaryValue[axis][side] = v; }
1169  fpreal getBoundaryValue(int axis, int side) const
1170  { return myBoundaryValue[axis][side]; }
1171 
1172  /// Mark this field as being an extrapolated field. Out of bound
1173  /// voxels will read the clamped value. The difference between
1174  /// the clamped position and the real position is then dot producted
1175  /// with the given scale factor and added to the resulting value.
1176  /// This allows you to have rest fields that extrapolate meaningfully.
1177  void setAsExtrapolatedField(UT_Vector3 scale);
1178 
1179  /// These adjust the native UT_VoxelArray border values that *are*
1180  /// used for reading outside the valid range.
1181  UT_VoxelBorderType getBorder() const;
1182  float getBorderValue() const;
1183  void setBorder(UT_VoxelBorderType border, float bval);
1184 
1185  /// These adjust the native UT_VoxelArray comrpression options.
1186  void setCompressionOptions(const UT_VoxelCompressOptions &options);
1187  const UT_VoxelCompressOptions &getCompressionOptions() const;
1188 
1189  /// These adjust the native UT_VoxelArray comrpression tolerance..
1190  void setCompressionTolerance(fpreal tol);
1191  fpreal getCompressionTolerance() const;
1192 
1193  /// Transform turns this into a packed set of wavelet coeffecients
1194  /// from the scalar data in field. Inverse unpacks and generates
1195  /// a scalar field into this.
1196  ///
1197  /// Using a smooth edge conditon we assume that
1198  /// in case of an odd sized field, the ultimate row has
1199  /// zero diff coefficients. We thus store the extra
1200  /// column in the averages side of the matrix, for
1201  /// ceil(n/2) in the averages and floor(n/2) in the diffs.
1202  /// Note this causes 3d fields that are actually 2d
1203  /// to properly decompose.
1204  void waveletTransform(UT_Wavelet::WAVELET_NAMES wavelettype,
1205  const SIM_RawField *field,
1206  int maxpasses = -1);
1207  void waveletInverseTransform(
1208  UT_Wavelet::WAVELET_NAMES wavelettype,
1209  const SIM_RawField *wavelet,
1210  int maxpasses = -1);
1211 
1212  /// Computes the sum of squares of the given level's detail vector.
1213  void waveletComputePSD(const SIM_RawField *wavelet,
1214  int level);
1215 
1216  /// Extracts the given component from a packed wavelet array.
1217  void waveletExtractComponent(const SIM_RawField *wavelet,
1218  int level,
1219  int component);
1220 
1221  /// Builds from either a GEO_PrimVolume or a GEO_PrimVDB.
1222  /// The provided volidx is the index inside of a vector volume
1223  /// to use for the building.
1224  void buildFromPrim(const GEO_Primitive *vol,
1225  int volidx,
1226  const UT_DMatrix4 &xform,
1227  fpreal scale);
1228 
1229  /// Builds from a GEO_PrimVolume
1230  void buildFromGeo(const GEO_PrimVolume *vol,
1231  const UT_DMatrix4 &xform,
1232  fpreal scale);
1233 
1234  /// Compute fractional volume weights representing the amount the voxel
1235  /// surrounding each sample point is inside the provided SDF. This method
1236  /// just subsamples the voxel space according to the samplesperaxis
1237  /// parameter.
1238 
1239 
1241  int samplesperaxis,
1242  bool invert,
1243  fpreal minweight,
1244  fpreal dilatedist = 0)
1245  {
1246  computeSDFWeightsSampledInternal(sdf,samplesperaxis,invert,
1247  minweight, dilatedist);
1248  }
1249 
1250 private:
1251  THREADED_METHOD5(SIM_RawField, shouldMultiThread(),
1252  computeSDFWeightsSampledInternal,
1253  const SIM_RawField *, sdf,
1254  int, samplesperaxis,
1255  bool, invert,
1256  fpreal, dilatedist,
1257  fpreal, minweight);
1258 
1259  void computeSDFWeightsSampledInternalPartial(const SIM_RawField *sdf,
1260  int samplesperaxis,
1261  bool invert,
1262  fpreal minweight,
1263  fpreal dilatedist,
1264  const UT_JobInfo &info);
1265 public:
1266  /// Compute fractional volume weights representing the amount the voxel
1267  /// surrounding each sample point is inside the provided SDF. This method
1268  /// uses the normalized area of the square that is in the plane given by
1269  /// the axis parameter as an approximation.
1270  THREADED_METHOD4(SIM_RawField, shouldMultiThread(),
1271  computeSDFWeightsFace,
1272  const SIM_RawField *, sdf,
1273  int, axis,
1274  bool, invert,
1275  fpreal, minweight);
1276 
1277  void computeSDFWeightsFacePartial(const SIM_RawField *sdf,
1278  int axis,
1279  bool invert,
1280  fpreal minweight,
1281  const UT_JobInfo &info);
1282 
1283  /// Compute fractional volume weights representing the amount the voxel
1284  /// surrounding each sample point is inside the provided SDF. This method
1285  /// uses accurate volume fractions.
1286  THREADED_METHOD3(SIM_RawField, shouldMultiThread(),
1287  computeSDFWeightsVolumeFraction,
1288  const SIM_RawField *, sdf,
1289  bool, invert,
1290  fpreal, minweight);
1291 
1292  void computeSDFWeightsVolumeFractionPartial(const SIM_RawField *sdf,
1293  bool invert,
1294  fpreal minweight,
1295  const UT_JobInfo &info);
1296 
1297 
1298  /// Compute fractional volume weights for the voxel specified by x, y, z.
1299  /// This method uses the normalized area of the square that is in the
1300  /// plane given by the axis parameter as an approximation.
1301  fpreal computeVoxelSDFWeightFace(int x, int y, int z,
1302  const SIM_RawField *sdf,
1303  int axis) const;
1304 
1305  /// Scale this field by s * B / C or set to zero if scaled value is lower than
1306  /// given threshold
1307  THREADED_METHOD4(SIM_RawField, shouldMultiThread(),
1308  setScaleDivideThreshold,
1309  fpreal, scale,
1310  const SIM_RawField *, B,
1311  const SIM_RawField *, C,
1312  fpreal, threshold)
1313  void setScaleDivideThresholdPartial(fpreal scale,
1315  const SIM_RawField *C,
1316  fpreal threshold,
1317  const UT_JobInfo &jobinfo);
1318 
1319 protected:
1320  void waveletRebuildFromVoxelArray(
1322  float scale);
1323 
1324 
1325  static fpreal applyBoundary(SIM_FieldBoundary bound, fpreal v,
1326  fpreal boundval);
1327 
1328  THREADED_METHOD2_CONST(SIM_RawField, shouldMultiThread(), reduceOpInternal,
1329  fpreal64 *, sum,
1330  REDUCE_NAMES, op)
1331  void reduceOpInternalPartial(fpreal64 *sum, REDUCE_NAMES op, const UT_JobInfo &info) const;
1332 
1333  THREADED_METHOD5_CONST(SIM_RawField, shouldMultiThread(), reduceMaskedOpInternal,
1334  fpreal64 *, sum,
1335  fpreal64 *, masktotal,
1336  REDUCE_NAMES, op,
1337  const SIM_RawField *, mask,
1338  bool, maskissdf)
1339  void reduceMaskedOpInternalPartial(fpreal64 *sum, fpreal64 *masktotal, REDUCE_NAMES op, const SIM_RawField *maskfield, bool maskissdf, const UT_JobInfo &info) const;
1340 
1341  void sortAllVoxels(UT_FloatArray &voxelvalues) const;
1342  void sortAllVoxelsMasked(UT_FloatArray &voxelvalues,
1343  const SIM_RawField *maskfield,
1344  bool maskissdf) const;
1345 
1346  /// Methods for extrapolation
1347  THREADED_METHOD8(SIM_RawField, shouldMultiThread(), buildExtrapList,
1349  const SIM_RawField *, depths,
1350  const SIM_RawField *, valid,
1351  fpreal, isocontour,
1352  fpreal, dir,
1353  fpreal, maxdist,
1354  bool, clamp,
1355  fpreal, clampval)
1356  void buildExtrapListPartial(
1357  UT_Array<UT_ValArray<gu_sdf_qelem *> > &lists,
1358  const SIM_RawField *depths,
1359  const SIM_RawField *valid,
1360  fpreal isocontour,
1361  fpreal dir, fpreal maxdist,
1362  bool clamp, fpreal clampval,
1363  const UT_JobInfo &info);
1364 
1365  /// Reduction by reducing each axis in turn.
1366  /// This *will* change field and dst != field is required.
1367  void localReduceByAxis(REDUCE_NAMES op,
1369  UT_VoxelArrayF &field,
1370  UT_Vector3 radius) const;
1371 
1372  /// Reduce along a given axis by the specified radius in voxels.
1373  template <int AXIS, REDUCE_NAMES OP>
1374  void localReduceAxis(UT_VoxelArrayF &dstfield,
1375  const UT_VoxelArrayF &field,
1376  float radius,
1377  const UT_JobInfo &info) const;
1378  template <int AXIS, REDUCE_NAMES OP>
1379  void localMinMaxAxis(UT_VoxelArrayF &dstfield,
1380  const UT_VoxelArrayF &field,
1381  float radius,
1382  const UT_JobInfo &info) const;
1383 
1384  template <int AXIS>
1385  void localReduceAxisOp(REDUCE_NAMES op,
1386  UT_VoxelArrayF &dstfield,
1387  UT_VoxelArrayF &field,
1388  float radius,
1389  const UT_JobInfo &info) const;
1390  /// Again a triple specialization to engage threading.
1391  THREADED_METHOD4_CONST(SIM_RawField,
1392  field.getTileRes(1)*field.getTileRes(2) > 1,
1393  localReduceAxisX,
1394  REDUCE_NAMES, op,
1395  UT_VoxelArrayF &, dst,
1396  UT_VoxelArrayF &, field,
1397  float, radius)
1398  void localReduceAxisXPartial(
1399  REDUCE_NAMES op,
1400  UT_VoxelArrayF &dst,
1401  UT_VoxelArrayF &field,
1402  float radius,
1403  const UT_JobInfo &info) const
1404  { localReduceAxisOp<0>(op, dst, field, radius, info); }
1405 
1407  field.getTileRes(0)*field.getTileRes(2) > 1,
1408  localReduceAxisY,
1409  REDUCE_NAMES, op,
1410  UT_VoxelArrayF &, dst,
1411  UT_VoxelArrayF &, field,
1412  float, radius)
1413  void localReduceAxisYPartial(
1414  REDUCE_NAMES op,
1416  UT_VoxelArrayF &field,
1417  float radius,
1418  const UT_JobInfo &info) const
1419  { localReduceAxisOp<1>(op, dst, field, radius, info); }
1420 
1422  field.getTileRes(0)*field.getTileRes(1) > 1,
1423  localReduceAxisZ,
1424  REDUCE_NAMES, op,
1425  UT_VoxelArrayF &, dst,
1426  UT_VoxelArrayF &, field,
1427  float, radius)
1428  void localReduceAxisZPartial(
1429  REDUCE_NAMES op,
1431  UT_VoxelArrayF &field,
1432  float radius,
1433  const UT_JobInfo &info) const
1434  { localReduceAxisOp<2>(op, dst, field, radius, info); }
1435 
1436  THREADED_METHOD3(SIM_RawField, shouldMultiThread(), buildFromGeoSampled,
1437  const GEO_PrimVolume *, vol,
1438  const UT_DMatrix4 &, xform,
1439  fpreal, scale)
1440  void buildFromGeoSampledPartial(const GEO_PrimVolume *vol,
1442  fpreal scale,
1443  const UT_JobInfo &info);
1444 
1445  THREADED_METHOD4(SIM_RawField, shouldMultiThread(), buildFromVDBSampled,
1446  const GEO_PrimVDB *, vdb,
1447  int, vdbidx,
1448  const UT_DMatrix4 &, xform,
1449  fpreal, scale)
1450  void buildFromVDBSampledPartial(const GEO_PrimVDB *vdb,
1451  int vdbidx,
1452  const UT_DMatrix4 &xform,
1453  fpreal scale,
1454  const UT_JobInfo &info);
1455 
1457 
1458  /// We always have myField at our current resolution.
1459  /// myGrid only exists if it matches myField. If myField
1460  /// changes, myGrid will be destroyed.
1461  /// If something updates the grid the field is out of date
1462  /// so is flagged.
1463  UT_VoxelArrayF *myField;
1464  mutable CE_Grid *myGrid;
1465  mutable bool myFieldOutOfDate;
1466 
1467  // This is the size and offset which converts our UT_VoxelArray into
1468  // world coordinates.
1469  UT_Vector3 myOrig, mySize;
1470 
1471  // This is our actual official size.
1472  UT_Vector3 myBBoxOrig, myBBoxSize;
1473 
1474  UT_Vector3 myVoxelSize;
1475  fpreal myVoxelDiameter;
1476 
1477  SIM_FieldBoundary myBoundary[3][2];
1478  fpreal myBoundaryValue[3][2];
1479 };
1480 
1481 
1482 ///
1483 /// Allows you to iterate over the voxel cells of a RawField. This
1484 /// is different than iterating over the voxel samples!
1485 /// The idea of getting or setting the raw value doesn't apply
1486 /// here as we are not necessarily alligned with the sample grid.
1487 ///
1489 {
1490 public:
1492  {
1493  myRes[0] = myRes[1] = myRes[2] = 0;
1494  rewind();
1495  }
1497  {
1498  }
1499 
1500  void setArray(const SIM_RawField *field)
1501  {
1502  field->getVoxelRes(myRes[0], myRes[1], myRes[2]);
1503  rewind();
1504  }
1505 
1506  void rewind()
1507  {
1508  myIdx[0] = 0;
1509  myIdx[1] = 0;
1510  myIdx[2] = 0;
1511  }
1512 
1513  bool atEnd() const
1514  {
1515  return myIdx[2] >= myRes[2];
1516  }
1517 
1518  void advance()
1519  {
1520  myIdx[0]++;
1521  if (myIdx[0] >= myRes[0])
1522  {
1523  myIdx[0] = 0;
1524  myIdx[1]++;
1525  if (myIdx[1] >= myRes[1])
1526  {
1527  myIdx[1] = 0;
1528  myIdx[2]++;
1529  if (myIdx[2] >= myRes[2])
1530  {
1531  // All done!
1532  myIdx[2] = myRes[2];
1533  }
1534  }
1535  }
1536  }
1537 
1538  int x() const { return myIdx[0]; }
1539  int y() const { return myIdx[1]; }
1540  int z() const { return myIdx[2]; }
1541  int idx(int axis) const { return myIdx[axis]; }
1542 
1543  /// Returns true if we are at the start of a new tile.
1544  /// Since we don't work on tile basis, we define a single
1545  /// x pass to be a tile.
1546  bool isStartOfTile() const
1547  {
1548  return myIdx[0] == 0;
1549  }
1550 
1551 protected:
1552  int myIdx[3];
1553  int myRes[3];
1554 };
1555 
1556 //
1557 // These macros create the CALL_VOXELPROBE macro
1558 //
1559 // CALL_VOXELPROBE(const SIM_RawField *src,
1560 // SIM_RawField *dst,
1561 // float defaultvalue,
1562 // callFunc(dst, theprobe))
1563 //
1564 // callFunc should be
1565 //
1566 // template <typename P>
1567 // void
1568 // callFunc(SIM_RawField *dst, P &probe)
1569 //
1570 // dst->isMatching(src) must be true.
1571 //
1572 // The callFunc will receive in the parameter "theprobe" a
1573 // UT_VoxelProbeAverage properly templated against the difference
1574 // between the source and dest fields.
1575 // theprobe.setIndex(x, y, z);
1576 // will cause theprobe.getValue() to return the interpolated value
1577 // in source of the destination's voxel sample (x, y, z). This takes
1578 // into account the different sampling patterns.
1579 //
1580 #define CALL_PROBEXY(XSTEP, YSTEP, src, dst, command) \
1581 { \
1582  if ((dstsample ^ srcsample) & 4) \
1583  { \
1584  if (dstsample & 4) /* ZStep -1 */ \
1585  { \
1586  UT_VoxelProbeAverage<fpreal32, XSTEP, YSTEP, -1> theprobe; \
1587  theprobe.setArray(src->field()); \
1588  command; \
1589  } \
1590  else /* ZStep 1 */ \
1591  { \
1592  UT_VoxelProbeAverage<fpreal32, XSTEP, YSTEP, 1> theprobe; \
1593  theprobe.setArray(src->field()); \
1594  command; \
1595  } \
1596  } \
1597  else /* ZStep 0 */ \
1598  { \
1599  UT_VoxelProbeAverage<fpreal32, XSTEP, YSTEP, 0> theprobe; \
1600  theprobe.setArray(src->field()); \
1601  command; \
1602  } \
1603 }
1604 
1605 #define CALL_PROBEX(XSTEP, src, dst, command) \
1606 { \
1607  if ((dstsample ^ srcsample) & 2) \
1608  { \
1609  if (dstsample & 2) /* YStep -1 */ \
1610  { \
1611  CALL_PROBEXY(XSTEP, -1, src, dst, command); \
1612  } \
1613  else /* YStep 1 */ \
1614  { \
1615  CALL_PROBEXY(XSTEP, 1, src, dst, command); \
1616  } \
1617  } \
1618  else /* YStep 0 */ \
1619  { \
1620  CALL_PROBEXY(XSTEP, 0, src, dst, command); \
1621  } \
1622 }
1623 
1624 #define CALL_VOXELPROBE(src, dst, defval, command) \
1625 { \
1626  if (!src) \
1627  { \
1628  UT_VoxelProbeConstant<fpreal32> theprobe; \
1629  theprobe.setValue(defval); \
1630  command; \
1631  } \
1632  else \
1633  { \
1634  int dstsample, srcsample; \
1635  dstsample = dst->getSample(); \
1636  srcsample = src->getSample(); \
1637  if ((dstsample ^ srcsample) & 1) \
1638  { \
1639  if (dstsample & 1) /* XStep -1 */ \
1640  { \
1641  CALL_PROBEX(-1, src, dst, command); \
1642  } \
1643  else /* XStep 1 */ \
1644  { \
1645  CALL_PROBEX(1, src, dst, command); \
1646  } \
1647  } \
1648  else /* XStep 0 */ \
1649  { \
1650  CALL_PROBEX(0, src, dst, command); \
1651  } \
1652  } \
1653 }
1654 
1655 #endif
#define SYSmax(a, b)
Definition: SYS_Math.h:1513
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
#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()
GLint GLint GLsizei GLint border
Definition: glcorearb.h:108
GLenum GLenum GLenum GLenum GLenum scale
Definition: glew.h:14163
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 UT_Vector3 & getVoxelSize() const
Definition: SIM_RawField.h:538
#define THREADED_METHOD1(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1)
static fpreal applyParticleToFieldParms(fpreal srcval, fpreal dstval, const sim_particleToFieldParms &parms)
Definition: SIM_RawField.h:933
exint getXRes() const
Definition: SIM_RawField.h:560
UT_Vector3T< float > UT_Vector3
static fpreal mixValues(MIX_NAMES mixtype, fpreal d, fpreal s)
Performs the requires mixing.
Definition: SIM_RawField.h:875
int64 exint
Definition: SYS_Types.h:125
GLint level
Definition: glcorearb.h:108
UT_VoxelArrayF * fieldNC() const
Definition: SIM_RawField.h:505
UT_VoxelBorderType
Definition: UT_VoxelArray.h:67
const SIM_RawField * velz
Definition: SIM_RawField.h:739
constexpr SYS_FORCE_INLINE T length() const noexcept
Definition: UT_Vector3.h:352
GLenum src
Definition: glcorearb.h:1793
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
GLuint GLuint GLfloat weight
Definition: glew.h:13892
#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:518
SIM_FieldSample
Definition: SIM_RawField.h:38
GLenum GLsizei GLsizei GLint * values
Definition: glcorearb.h:1602
GLint GLenum GLint x
Definition: glcorearb.h:409
GLsizeiptr size
Definition: glcorearb.h:664
GLubyte GLubyte GLubyte GLubyte w
Definition: glcorearb.h:857
double fpreal64
Definition: SYS_Types.h:201
sim_PointVelocityCB cb
Definition: SIM_RawField.h:104
GLenum array
Definition: glew.h:9108
fpreal getVoxelDiameter() const
Definition: SIM_RawField.h:541
void markGridAsChanged()
Mark the field as out of date, but only if we have a valid grid.
Definition: SIM_RawField.h:514
SIM_FieldSample getSample() const
Definition: SIM_RawField.h:547
GLuint * lists
Definition: glew.h:12689
sim_PointVelocity(sim_PointVelocityCB _cb, void *_data)
Definition: SIM_RawField.h:99
exint getZRes() const
Definition: SIM_RawField.h:562
GLuint64EXT * result
Definition: glew.h:14311
const UT_Vector3 & getBBoxOrig() const
Definition: SIM_RawField.h:536
exint getYRes() const
Definition: SIM_RawField.h:561
GLfloat GLfloat GLfloat v2
Definition: glcorearb.h:818
GLint GLuint mask
Definition: glcorearb.h:124
const UT_Vector3 & getBBoxSize() const
Definition: SIM_RawField.h:537
const GLdouble * v
Definition: glcorearb.h:837
GLsizei GLsizei GLchar * source
Definition: glcorearb.h:803
#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
const SIM_RawField * collision
Definition: SIM_RawField.h:741
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:848
void setBoundaryValue(int axis, int side, fpreal v)
long long int64
Definition: SYS_Types.h:116
GLfloat GLfloat p
Definition: glew.h:16656
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
const UT_Vector3 & getSize() const
Definition: SIM_RawField.h:535
UT_Vector3 planeNormal
Definition: SIM_RawField.h:90
const SIM_RawField * source
Definition: SIM_RawField.h:738
#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
GLfloat GLfloat GLfloat GLfloat v3
Definition: glcorearb.h:819
#define THREADED_METHOD4(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3, PARMTYPE4, PARMNAME4)
void setArray(const SIM_RawField *field)
GLint GLsizei GLsizei height
Definition: glcorearb.h:103
fpreal getBoundaryValue(int axis, int side) const
const UT_VoxelArrayF * field() const
Fetches the raw field.
Definition: SIM_RawField.h:504
#define THREADED_METHOD2_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2)
GLuint GLfloat * val
Definition: glcorearb.h:1608
GA_API const UT_StringHolder up
fpreal64 fpreal
Definition: SYS_Types.h:277
GLuint index
Definition: glcorearb.h:786
SIM_FieldAdvection
Definition: SIM_RawField.h:60
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
#define UT_ASSERT(ZZ)
Definition: UT_Assert.h:153
Definition: core.h:1131
fpreal getVoxelVolume() const
Definition: SIM_RawField.h:542
bool hasNan() const
Returns true if our field has any NANs.
Definition: SIM_RawField.h:497
GLintptr offset
Definition: glcorearb.h:665
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:744
#define const
Definition: zconf.h:214
GLfloat GLfloat v1
Definition: glcorearb.h:817
GLint GLfloat GLint stencil
Definition: glcorearb.h:1278
REDUCE_NAMES
Types of reductions supported by reduceOp.
Definition: SIM_RawField.h:425
const UT_Vector3 & getOrig() const
Definition: SIM_RawField.h:534
GLenum GLenum dst
Definition: glcorearb.h:1793
GLboolean r
Definition: glcorearb.h:1222
#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:1514
This class holds a three dimensional vector field.
GLdouble s
Definition: glew.h:1395
GLint y
Definition: glcorearb.h:103
#define THREADED_METHOD(CLASSNAME, DOMULTI, METHOD)
SIM_FieldBoundary getBoundary(int axis, int side) const
void setVoxelSize(const UT_Vector3 &voxelsize)
Definition: SIM_RawField.h:539
Definition: format.h:895
SIM_FieldBoundary
Definition: SIM_RawField.h:50