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