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