SIM/SIM_RawField.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * Produced by:
 *      Jeff Lait
 *      Side Effects Software Inc
 *      123 Front Street West, Suite 1401
 *      Toronto, Ontario
 *      Canada   M5J 2M2
 *      416-504-9876
 *
 * NAME:        SIM_RawField.h ( SIM Library, C++)
 *
 * COMMENTS:
 */

#ifndef __SIM_RawField__
#define __SIM_RawField__

#include "SIM_API.h"

#include <UT/UT_VoxelArray.h>
#include <UT/UT_ThreadedAlgorithm.h>
#include <UT/UT_Vector.h>
#include <UT/UT_FloatArray.h>
#include <UT/UT_SparseMatrix.h>
#include <UT/UT_Wavelet.h>


// For 32 bit windows it is important to use a paged sparse matrix
// as vista's memory model doesn't let you allocate large continguous
// structures.  There is a performance penalty for this so we only
// do it on that platform.
#if defined(WIN32) && !defined(AMD64)
typedef UT_SparseMatrixT<fpreal32, true> SIM_RawSparseMatrix;
#else
typedef UT_SparseMatrixF SIM_RawSparseMatrix;
#endif

class SIM_VectorField;
class SIM_RawField;
class SIM_RawIndexField;
class SIM_PartitionedRawIndexField;
class GEO_PrimVolume;
class UT_IntArray;

// The order of this enum is important.  It is a 3 bit array
// where each bit determines if the center value or corner value
// should be used on that axis.  0 means center, 1 means node.
// The least significant bit is x axis, then y, then z.
enum SIM_FieldSample
{
    SIM_SAMPLE_CENTER = 0,
    SIM_SAMPLE_FACEX = 1,
    SIM_SAMPLE_FACEY = 2,
    SIM_SAMPLE_EDGEXY = 3,
    SIM_SAMPLE_FACEZ = 4,
    SIM_SAMPLE_EDGEXZ = 5,
    SIM_SAMPLE_EDGEYZ = 6,
    SIM_SAMPLE_CORNER = 7
};

enum SIM_FieldBoundary
{
    SIM_BOUNDARY_NONE,  // Unchanged, free to change.
    SIM_BOUNDARY_SAME,  // Same as adjacent voxels.
    SIM_BOUNDARY_ZERO,  // Set to constant value, myBoundaryValue
    SIM_BOUNDARY_NEGATE,        // Opposite sign of adjacent voxels.
    SIM_BOUNDARY_FIXED, // Unaffected by enforcement, but should be treated
                        // as a fixed boundary by projection methods.
};

enum SIM_FieldAdvection
{
    SIM_ADVECT_SINGLE = 0,      // Take a single step
    SIM_ADVECT_TRACE = 1,       // Trace according to velocity field
    SIM_ADVECT_MIDPOINT = 2     // Trace using midpoint updates
};

// Type definition for a velocity function
typedef UT_Vector3 (*sim_PointVelocity)(const UT_Vector3 &, int);

class SIM_API SIM_RawField
{
public:
                SIM_RawField();
    virtual    ~SIM_RawField();

    /// Copy constructor:
                SIM_RawField(const SIM_RawField &src);

    /// Assigment operator:
    const SIM_RawField  &operator=(const SIM_RawField &src);            

    /// Initializes the field.
    /// The resolution given is in terms of voxels, the actual dimensions
    /// of this field may be slightly different due to the sampling
    /// choice.
    void        init(SIM_FieldSample sample, 
                     const UT_Vector3 &orig, const UT_Vector3 &size,
                     int xres, int yres, int zres);

    /// Initializes the voxel iterator to only run over a subset
    /// of tiles that correspond to the appropriate range for
    /// the jobinfo.
    void        getPartialRange(UT_VoxelArrayIteratorF &vit,
                            const UT_JobInfo &info) const;

    /// Returns true if this should be multithreaded.
    bool        shouldMultiThread() const
                {
#ifdef CELLBE
                    return false;
#else
                    return field()->numTiles() > 1;
#endif
                }

    /// Initializes this to be the same dimensions, sampling pattern,
    /// etc, of the given field.
    /// The values of this may be reset to zero.
    void        match(const SIM_RawField &src);
    void        match(const SIM_RawIndexField &src);

    /// Initializes this to be the same as the given field.  The
    /// values of this will be found by resampling the other field
    /// with the new sampling pattern.
    void        resample(SIM_FieldSample sample, const SIM_RawField *src);

    /// Returns the set of samples in *this* field which correspond
    /// to the given location & sampling pattern.
    /// ix, iy, and iz should be size 8.
    /// If you want the deltas for the sampling pattern, call with
    /// x, y, z zero and clamp to false.
    void        getSamplePattern(SIM_FieldSample sample, int x, int y, int z,
                            int &numsample, int *ix, int *iy, int *iz,
                            bool clamp) const;

    /// Extrapolates the values of this field to all locations which are
    /// on the wrong side of the isocontour.  dir == 1 means greater than
    /// isocontour is wrong, dir == -1 is less than.
    /// maxdist, if not negative, represents the greatest *absolute*
    /// distance to extrapolate to.  Values outside of this are left
    /// unaffected, unless clamp is set, in which case they are set to
    /// clampval.
    void        extrapolate(const SIM_RawField *depths, fpreal isocontour, 
                            fpreal dir, fpreal maxdist,
                            bool clamp, fpreal clampval);

    /// Advances source field along the field's curvature.
    /// Takes a given b to act as the movement coefficient.
    /// Does a single euler step.
    /// this is filled in with the resulting field.  this cannot be source.
    /// This solves:
    ///  dphi/dt = b * K * |grad(phi)|
    THREADED_METHOD3(SIM_RawField, shouldMultiThread(), moveAlongCurvature,
                    fpreal, b_val,
                    const SIM_RawField &, source,
                    fpreal, timestep)
    void         moveAlongCurvaturePartial(fpreal b_val, 
                            const SIM_RawField &source,
                            fpreal timestep,
                            const UT_JobInfo &jobinfo);

    /// Advances this field along the field's curvature.
    /// Is given a total amount to move and will calculate the appropriate
    /// timestep according to the given cfl condition & min/max iterations.
    void         moveAlongCurvature(fpreal amount,
                            fpreal cflcond,
                            int miniter, int maxiter);

    /// Advances source field along the field's normal direction.
    /// Takes a matching field which defines the speed to move in the
    /// normal direction on a per cell basis.
    /// Does a single euler step.
    /// this is filled in with the resulting field.  this cannot
    /// be source.
    /// This solves:
    ///  dphi/dt + speed * |grad(phi)| = 0
    THREADED_METHOD3(SIM_RawField, shouldMultiThread(), moveAlongNormal,
                     const SIM_RawField &, speed,
                     const SIM_RawField &, source,
                     fpreal, timestep)
    void        moveAlongNormalPartial(const SIM_RawField &speed, 
                                const SIM_RawField &source, 
                                fpreal timestep, const UT_JobInfo &info);

    /// Uses the midpoint method to do a second order temporal
    /// update of the moveAlongNormal algorithm.
    void        moveAlongNormalMidpoint(const SIM_RawField &speed,
                                const SIM_RawField &source,
                                fpreal timestep);

    /// Performs the reinitialization equation.
    /// This solves for:
    ///  dphi/dt + S(phi0) * (|grad(phi)| - 1) = 0
    /// at steady state.
    ///  S(phi0) = phi0 / sqrt(phi0^2 + dx^2), a smeared sign function.
    /// It is held constant for the reinitialization.
    void        reinitializeSignedDistance(int maxiter);

    /// Negates all values in this field
    THREADED_METHOD(SIM_RawField, shouldMultiThread(), negate)
    void        negatePartial(const UT_JobInfo &info);

    /// Makes this field the minimum or maximum 
    /// of this field and the other field.  The other field
    /// is not assumed to be matching.
    THREADED_METHOD1(SIM_RawField, shouldMultiThread(), maximum,
                        const SIM_RawField *, other)
    void        maximumPartial(const SIM_RawField *other, const UT_JobInfo &info);

    THREADED_METHOD1(SIM_RawField, shouldMultiThread(), minimum, 
                    const SIM_RawField *, other)
    void        minimumPartial(const SIM_RawField *other, const UT_JobInfo &jobinfo);

    /// Averages this with the other field.  Result in this.  Assumes
    /// fields match.
    THREADED_METHOD1(SIM_RawField, shouldMultiThread(), average,
                    const SIM_RawField &, other)
    void        averagePartial(const SIM_RawField &other, const UT_JobInfo &jobinfo);   

    /// Sets this to A + s*B.  A & B match this.  This can be either
    /// A or B, or neither.
    THREADED_METHOD3(SIM_RawField, shouldMultiThread(), setScaleAdd, 
                    const SIM_RawField &, A,
                    fpreal, scale,
                    const SIM_RawField &, B)
    void        setScaleAddPartial(const SIM_RawField &A,
                            fpreal scale,
                            const SIM_RawField &B, const UT_JobInfo &jobinfo);

    /// Sets this to the smeared sign function of the given sdf.
    ///  this = sdf / sqrt(sdf*sdf + bandwidth*bandwidth)
    /// this must already match sdf.
    THREADED_METHOD2(SIM_RawField, shouldMultiThread(), smearedSign,
                        const SIM_RawField &, sdf,
                        fpreal, bandwidth)
    void        smearedSignPartial(const SIM_RawField &sdf,
                                    fpreal bandwidth,
                                    const UT_JobInfo &info);

    /// Sets this to the heaviside function of itself.
    /// this = clamp(-this/diam+0.5, 0, 1)
    THREADED_METHOD(SIM_RawField, shouldMultiThread(), convertToHeaviside)
    void        convertToHeavisidePartial(const UT_JobInfo &info);

    /// Perform heaviside & inverse heaviside in a consistent fashion.
    static fpreal       toHeaviside(fpreal val, fpreal diam);
    static fpreal       fromHeaviside(fpreal val, fpreal diam);
    /// Defines width as maxcomponent of our voxels.
    fpreal              toHeaviside(fpreal val) const;
    fpreal              fromHeaviside(fpreal val) const;

    /// Types of reductions supported by reduceOp.
    /// *not* necessarily matching GAS_Reduce!
    enum REDUCE_NAMES
    {
        REDUCE_MAX,
        REDUCE_MIN,
        REDUCE_AVERAGE,
        REDUCE_SUM,
        REDUCE_SUMABS,
        REDUCE_SUMSQUARE,
        REDUCE_RMS
    };
    
    /// Perform reduction on the field using the given method
    /// The internal methods are helpers to deal with the threading
    /// collating.
    fpreal      reduceOp(REDUCE_NAMES op) const;

    /// Prebuilt versions.
    fpreal      average() const { return reduceOp(REDUCE_AVERAGE); }

    /// Performs the reduction only including parts of the field
    /// that match the given mask.  If mask is null, falls through
    /// to normal reduction.
    fpreal      reduceMaskedOp(REDUCE_NAMES op, 
                        const SIM_RawField *mask, bool maskissdf) const;

    /// Sorts the valid voxels (as defined by optional mask) and returns
    /// the given percentile voxel.  0.5 means median. 0.25 first quartile,
    /// etc.
    fpreal      findProportionalValue(fpreal position, 
                        const SIM_RawField *mask, bool maskissdf) const;

    /// Returns true if the given field and this one match in terms
    /// of number of voxels and bounding box size.
    /// This means the voxel cells match - not necessarily the sample
    /// points!
    bool        isMatching(const SIM_RawField *field) const;
    bool        isMatching(const SIM_RawIndexField *field) const;

    /// Returns true if the two fields are precisely aligned.  This
    /// means that samples are matched so a given integer index
    /// into either field would give the same result.
    bool        isAligned(const SIM_RawField *field) const;
    bool        isAligned(const SIM_RawIndexField *field) const;

    /// Fetches the raw field.
    UT_VoxelArrayF              *field() const { return myField; }

    const UT_Vector3            &getOrig() const { return myOrig; }
    const UT_Vector3            &getSize() const { return mySize; }
    const UT_Vector3            &getBBoxOrig() const { return myBBoxOrig; }
    const UT_Vector3            &getBBoxSize() const { return myBBoxSize; }
    const UT_Vector3            &getVoxelSize() const { return myVoxelSize; }
    fpreal                       getVoxelDiameter() const { return myVoxelDiameter; }
    fpreal                       getVoxelVolume() const
                                 { return myVoxelSize.x() 
                                        * myVoxelSize.y() 
                                        * myVoxelSize.z(); }

    SIM_FieldSample              getSample() const { return mySample; }

    int64                        getMemoryUsage() const;

    /// Returns the resolution of the voxel grid that we are sampling.
    /// This is a count of voxels, so may differ for our different
    /// sampling methods.
    void                 getVoxelRes(int &xres, int &yres, int &zres) const;

    /// Functions to resolve quantities about the field.

    /// Sample the field at the given world space location.
    fpreal               getValue(UT_Vector3 pos) const;

    /// Returns an averaged value for the center of the given voxel.
    fpreal               getCellValue(int x, int y, int z) const;

    /// Adds the value to the given voxel cell, averaging out among
    /// adjacent samples if we aren't sampled evenly.
    void                 addToCell(int x, int y, int z, fpreal v);

    /// Ensures the given voxel cell has the given value.  This
    /// will set all of the adjacent samples if we aren't sampled
    /// evenly.
    void                 setCellValue(int x, int y, int z, fpreal v);

    /// Return the gradient of the field at the given world space location.
    /// Uses central differencing with a sample spacing of the voxelsize.
    UT_Vector3           getGradient(UT_Vector3 pos) const;

    /// Returns the gradient at the given voxel index.
    /// Uses central differencing.
    UT_Vector3           getGradientAtIndex(int x, int y, int z) const;

    /// Computes the laplacian of the field at the world space
    /// coordinate using interpolation of getLaplacianAtIndex
    fpreal64             getLaplacian(UT_Vector3 pos) const;

    /// Computes the laplacian of the field at the specific voxel index
    /// specified
    fpreal64             getLaplacianAtIndex(int x, int y, int z) const;

    /// Returns the curvature of the field at a given world space location.
    /// Uses interpolation of index based getCurvatureAtIndex.
    fpreal64             getCurvature(UT_Vector3 pos) const;

    /// Computes the curvature at the given voxel index.
    /// Uses central differencing.
    /// The resulting curvature is clamped according to the resolution
    /// of the field to avoid abnormally large values in noisy data.
    fpreal64             getCurvatureAtIndex(int x, int y, int z) const;

    /// Computes the curvature according to a 3^3 voxel probe
    static fpreal64      getCurvatureAtProbe(UT_VoxelProbeCubeF &probe, const UT_Vector3 &invvoxelsize);

    /// Computes K |grad(phi|), the curvature times the length of the
    /// gradient.  By folding the two operations, we can return a non-zero
    /// value where grad goes to zero (but curvature doesn't) using
    /// L'Hospital's rule.  This also does not clamp the curvature.
    fpreal64             getCurvatureTimesGradAtIndex(int x, int y, int z) const;

    /// Giving the relevent divided differences, compute the HJWENO
    /// approxmiation of the derivative.
    fpreal64             calculateHJWENO(fpreal64 v1, fpreal64 v2, fpreal64 v3, fpreal64 v4, fpreal64 v5) const;
    
    /// Calculate the derivitive along the specified axis using
    /// an HJ WENO method at the given index.  The boolean
    /// determines whether the derivitive in the positive or negative
    /// direction is computed.
    fpreal64             calculateDerivative(int x, int y, int z,
                                    int axis, bool positivegradient) const;

    /// Convert indices to world coordinates and vice-versa.  Note this uses
    /// this field's indices which change depending on sampling.
    bool                 indexToPos(int x, int y, int z, UT_Vector3 &pos) const;

    /// Converts a worldspace position into an integer index.
    bool                 posToIndex(UT_Vector3 pos, int &x, int &y, int &z) const;

    /// Converts worldspace position into an integer index + the lerp values
    /// required to interpolate.
    /// Lerp with (1-dx) * (x,y,z) + dx * (x+1,y,z)
    bool                 posToIndex(UT_Vector3 pos, int &x, int &y, int &z,
                                    fpreal &dx, fpreal &dy, fpreal &dz) const;
    /// Convert voxel cell indices to world coordinates and vice-versa.
    /// Returns values at cell centers.  Is equivalent to the indexToPos
    /// style functions only when sampling is CENTER.
    bool                 cellIndexToPos(int x, int y, int z, UT_Vector3 &pos) const;
    bool                 posToCellIndex(UT_Vector3 pos, int &x, int &y, int &z) const;

    /// Verbs that can be performed on these fields.

    /// Advect a point in space according to an array of velocity
    /// fields.
    static void          advect(UT_Vector3 &pos, 
                                const SIM_RawField *velx, 
                                const SIM_RawField *vely, 
                                const SIM_RawField *velz, fpreal time,
                                const SIM_RawField *collision = 0);
    static void          advect(UT_Vector3 &pos,
                                sim_PointVelocity getVelocity,
                                fpreal time, fpreal voxelsize,
                                int jobnum = 0,
                                const SIM_RawField *collision = 0);

    /// Advect a point with the midpoint method.
    static void         advectMidpoint(UT_Vector3 &pos,
                                const SIM_RawField *velx, 
                                const SIM_RawField *vely, 
                                const SIM_RawField *velz, fpreal time,
                                const SIM_RawField *collision = 0);
    static void         advectMidpoint(UT_Vector3 &pos,
                                sim_PointVelocity getVelocity,
                                fpreal time, fpreal voxelsize,
                                int jobnum = 0,
                                const SIM_RawField *collision = 0);

    /// Move a point to the given isooffset.
    /// Returns false if fails to complete the move in the given time.
    bool                movePtToIso(UT_Vector3 &pos,
                                    fpreal goaliso,
                                    fpreal maxtime,
                                    fpreal tol = 1e-4) const;

    /// Advect the source field by the given set of velocity fields,
    /// setting this to the result.  Source cannot be this.
    /// this and source are assumed to match.
    THREADED_METHOD7(SIM_RawField, shouldMultiThread(), advect,
                    const SIM_RawField *, source,
                    const SIM_RawField *, velx,
                    const SIM_RawField *, vely,
                    const SIM_RawField *, velz,
                    fpreal, time,
                    const SIM_RawField *, collision,
                    SIM_FieldAdvection, advectmethod
                    )
    void                 advectPartial(const SIM_RawField *source,
                                const SIM_RawField *velx, 
                                const SIM_RawField *vely, 
                                const SIM_RawField *velz, 
                                fpreal time,
                                const SIM_RawField *collision,
                                SIM_FieldAdvection advectmethod,
                                const UT_JobInfo &info);

    THREADED_METHOD4(SIM_RawField, shouldMultiThread(), 
                        buoyancy,
                        const SIM_RawField *, temperature,
                        fpreal, up, 
                        fpreal, Tamb,
                        fpreal, buoyancy)
    void                 buoyancyPartial(const SIM_RawField *temperature,
                                  fpreal up, fpreal Tamb,
                                  fpreal buoyancy,
                                  const UT_JobInfo &info);

    THREADED_METHOD5(SIM_RawField, shouldMultiThread(), advect2,
                    const SIM_RawField *, source,
                    sim_PointVelocity, getVelocity,
                    fpreal, time, fpreal, voxelsize,
                    const SIM_RawField *, collision)
    void                advect2Partial(const SIM_RawField *source,
                                sim_PointVelocity getVelocity,
                                fpreal time, fpreal voxelsize,
                                const SIM_RawField *collision,
                                const UT_JobInfo &info);

    /// Advect this field by the given velocity fields.  Invokes
    /// advect but handles the creation of the intermediate field.
    void                 advectSelf(const SIM_RawField *velx, 
                                const SIM_RawField *vely, 
                                const SIM_RawField *velz, 
                                fpreal time,
                                const SIM_RawField *collision,
                                SIM_FieldAdvection advectmethod);
    void                 advectSelf(sim_PointVelocity getVelocity,
                                fpreal time, fpreal voxelsize,
                                const SIM_RawField *collision);

    /// Solve the diffusion equation over this field according to
    /// the given diffusion rate.  One likely wants to roll the
    /// desired timestep into the diffusion rate.
    void                 diffuse(fpreal diffrate,
                                int numiter,
                                const SIM_RawField *collision=0);

    /// Raw diffuse algorithm, exposed only for external performance tests
    static void          diffuse(float *dstdata, const float *srcdata[3][3][3],
                                 float b, float ivsx, float ivsy, float ivsz,
                                 int tx, int ty,
                                 int max_xtile, int max_ytile,
                                 int max_xvox, int max_yvox);

    /// Determine all components connected according to < 0 semantic
    void                 computeConnectedComponents(UT_VoxelArray<int64> &comp, int &numcomponent) const;

    /// Computes the divergence of the face-centered velocity field
    /// and stores the result in this.  this must match the velocity field.
    THREADED_METHOD2(SIM_RawField, shouldMultiThread(), buildDivergenceFace,
                        const SIM_VectorField *, vel,
                        const SIM_RawField *, surface)
    void                 buildDivergenceFacePartial(const SIM_VectorField *vel,
                                    const SIM_RawField *surface,
                                    const UT_JobInfo &info);

    enum PCG_METHOD
    {
        PCG_NONE,
        PCG_JACOBI,
        PCG_CHOLESKY,
        PCG_MIC
    };

    /// Solves the pressure field that would eliminate the given
    /// divergence field.
    void                 solvePressure(const SIM_RawField *divergence,
                                        const SIM_RawField *collision,
                                        int numiter = 20);
    void                 solvePressurePCG(const SIM_RawField &divergence,
                                        SIM_PartitionedRawIndexField &index,
                                        SIM_VectorField *vel,
                                        const SIM_RawIndexField *comp = 0,
                                        const UT_IntArray *expandable = 0,
                                        const SIM_RawField *surface = 0,
                                        bool variational = true,
                                        bool ghostfluid = true,
                                        PCG_METHOD pcgmethod = PCG_MIC);

    // For each voxel v in *this and corresponding voxel b in *B
    // sv = sum of 6-way neighbours of v
    // v' = (sv * sumweight + b) * weight
    // This is only applied to the voxels whose parity (x^y^z) matches
    // parity, so should be called twice to complete a single iteration.
    THREADED_METHOD4(SIM_RawField, shouldMultiThread(),
            gaussSeidelIteration,
            const SIM_RawField *, B,
            fpreal32,  weight,
            fpreal32,  sumweight,
            int,  parity)
    void gaussSeidelIterationPartial(
            const SIM_RawField * B,
            fpreal32  weight,
            fpreal32  sumweight,
            int  parity,
            const UT_JobInfo &info);

    // Note that in the flat case the destination field is passed
    // parameter, as opposed to the previous function where it is
    // the B field that is a parameter.
    THREADED_METHOD4_CONST(SIM_RawField, shouldMultiThread(),
            gaussSeidelIterationFlat,
            fpreal32 *, A,
            fpreal32,  weight,
            fpreal32,  sumweight,
            int,  parity)
    void gaussSeidelIterationFlatPartial(
            fpreal32 * A,
            fpreal32  weight,
            fpreal32  sumweight,
            int  parity,
            const UT_JobInfo &info) const;

    // Note that in the flat case the destination field is passed
    // parameter, as opposed to the previous function where it is
    // the B field that is a parameter.
    // This does a 4-way neighbours.
    THREADED_METHOD4_CONST(SIM_RawField, shouldMultiThread(),
            gaussSeidelIterationFlat2D,
            fpreal32 *, A,
            fpreal32,  weight,
            fpreal32,  sumweight,
            int,  parity)
    void gaussSeidelIterationFlat2DPartial(
            fpreal32 * A,
            fpreal32  weight,
            fpreal32  sumweight,
            int  parity,
            const UT_JobInfo &info) const;

    // Solve independent systems of equations
    //
    //     ld[s] ys[s] = bs[s],
    //
    // where each of the matrices ld[s] is lower-triangular.
    // This function is used by the threaded implementation of the
    // incomplete cholesky preconditioner.
    THREADED_METHOD3_CONST(SIM_RawField, shouldMultiThread(),
        solveLowerTriangularSystems,
        const UT_PtrArray<SIM_RawSparseMatrix*>&, ld,
        UT_RefArray<UT_VectorF>&, ys,
        const UT_RefArray<UT_VectorF>&, bs
    )
    void solveLowerTriangularSystemsPartial(
        const UT_PtrArray<SIM_RawSparseMatrix*>& ld,
        UT_RefArray<UT_VectorF>& ys,
        const UT_RefArray<UT_VectorF>& bs,
        const UT_JobInfo& info
    ) const;

    // Solve independent systems of equations
    //
    //     ud[s] xs[s] + us[s] xc = ys[s],
    //
    // where each ud[s] is upper triangular.
    // This function is used by the threaded implementation of the
    // incomplete cholesky preconditioner.
    // The vector proxied by the ys's is modified by this function!
    THREADED_METHOD5_CONST(SIM_RawField, shouldMultiThread(),
        solveUpperTriangularSystems,
        const UT_PtrArray<SIM_RawSparseMatrix*>&, ud,
        UT_RefArray<UT_VectorF>&, xs,
        const UT_PtrArray<SIM_RawSparseMatrix*>&, us,
        const UT_VectorF&, xc,
        UT_RefArray<UT_VectorF>&, ys
    )
    void solveUpperTriangularSystemsPartial(
        const UT_PtrArray<SIM_RawSparseMatrix*>& ud,
        UT_RefArray<UT_VectorF>& xs,
        const UT_PtrArray<SIM_RawSparseMatrix*>& us,
        const UT_VectorF& xc,
        UT_RefArray<UT_VectorF>& ys,
        const UT_JobInfo& info
    ) const;

    /// Applies the incomplete cholesky preconditioner.
    THREADED_METHOD2_CONST(SIM_RawField, shouldMultiThread(),
                    incompleteCholeskyPreconditioner,
                    const UT_VectorF &, b,
                    UT_VectorF &, x)
    void                 incompleteCholeskyPreconditionerPartial(
                            const UT_VectorF &b,
                            UT_VectorF &x,
                            const UT_JobInfo &info) const;

    /// Calculates the incomplete cholesky factorization.
    THREADED_METHOD2_CONST(SIM_RawField, shouldMultiThread(),
                    incompleteCholeskyFactorization,
                    const SIM_RawSparseMatrix &, A,
                    PCG_METHOD, pcgmethod)
    void                incompleteCholeskyFactorizationPartial(
                            const SIM_RawSparseMatrix &A,
                            PCG_METHOD pcgmethod,
                            const UT_JobInfo &info) const;
        
    /// Enforces the boundary conditions on this field.
    /// Each axis can have its own conditions
    void                 enforceBoundary(SIM_FieldBoundary collisionboundary = SIM_BOUNDARY_NONE, 
                                        const SIM_RawField *collision=0,
                                        const SIM_RawField *cvalue = 0);
    void                 enforceCollisionBoundary(SIM_FieldBoundary boundary,
                                        const SIM_RawField *collision,
                                        const SIM_RawField *cvalue = 0);
    void                 enforceSideBoundary(int axis, int side,
                                        SIM_FieldBoundary bound,
                                         fpreal boundaryvalue);

    THREADED_METHOD3(SIM_RawField, shouldMultiThread(),
                        enforceCollisionBoundaryInternal,
                        SIM_FieldBoundary, boundary,
                        const SIM_RawField *, collision,
                        const SIM_RawField *, cvalue)
    void                 enforceCollisionBoundaryInternalPartial(
                                SIM_FieldBoundary boundary,
                                const SIM_RawField *collision,
                                const SIM_RawField *cvalue,
                                const UT_JobInfo &info);

    /// Enforces the boundary conditions onto a flat array
    /// using the dimensions of this.
    /// Does SAME boundaries for collision objects.
    /// Index field is where to copy from for each collision voxel.
    void                enforceBoundaryFlat(fpreal32 *values,
                                const SIM_RawIndexField *collision_lookup);
    THREADED_METHOD2(SIM_RawField, shouldMultiThread(), 
                        enforceCollisionBoundaryFlat,
                        fpreal32 *, values,
                        const SIM_RawIndexField *, collision_lookup)

    void                enforceCollisionBoundaryFlatPartial(fpreal32 *values,
                                const SIM_RawIndexField *collision_lookup,
                                const UT_JobInfo &info);
    void                enforceSideBoundaryFlat(fpreal32 *values,
                                int axis, int side,
                                SIM_FieldBoundary bound,
                                fpreal boundval);

    /// These boundary conditions do *not* apply to reading outside
    /// of the valid field range.  The native UT_VoxelArray boundary
    /// condition is used for that.
    /// Instead, they are used for the behaviour of enforceBoundary
    /// and by various places where we want to distinguish open 
    /// velocity fields from closed.

    void                 setBoundary(int axis, int side, SIM_FieldBoundary bound)
                                    { myBoundary[axis][side] = bound; }
    SIM_FieldBoundary    getBoundary(int axis, int side) const 
                                    { return myBoundary[axis][side]; }
    void                 setBoundaryValue(int axis, int side, fpreal v)
                                    { myBoundaryValue[axis][side] = v; }
    fpreal               getBoundaryValue(int axis, int side) const
                                    { return myBoundaryValue[axis][side]; }

    /// Transform turns this into a packed set of wavelet coeffecients
    /// from the scalar data in field.  Inverse unpacks and generates
    /// a scalar field into this.
    ///
    /// Using a smooth edge conditon we assume that
    /// in case of an odd sized field, the ultimate row has
    /// zero diff coefficients.   We thus store the extra
    /// column in the averages side of the matrix, for
    /// ceil(n/2) in the averages and floor(n/2) in the diffs.
    /// Note this causes 3d fields that are actually 2d
    /// to properly decompose.
    void                 waveletTransform(UT_Wavelet::WAVELET_NAMES wavelettype,
                                        const SIM_RawField *field,
                                        int maxpasses = -1);
    void                 waveletInverseTransform(
                                        UT_Wavelet::WAVELET_NAMES wavelettype,
                                        const SIM_RawField *wavelet,
                                        int maxpasses = -1);

    /// Computes the sum of squares of the given level's detail vector.
    void                 waveletComputePSD(const SIM_RawField *wavelet,
                                            int level);

    /// Extracts the given component from a packed wavelet array.
    void                 waveletExtractComponent(const SIM_RawField *wavelet,
                                            int level,
                                            int component);

    THREADED_METHOD2(SIM_RawField, shouldMultiThread(), buildFromGeo,
                    const GEO_PrimVolume *, vol,
                    const UT_DMatrix4 &, xform)
    void                 buildFromGeoPartial(const GEO_PrimVolume *vol,
                                    const UT_DMatrix4 &xform,
                                    const UT_JobInfo &info);

protected:
    void                 waveletRebuildFromVoxelArray(
                                UT_VoxelArrayF *array,
                                float scale);


    static fpreal        applyBoundary(SIM_FieldBoundary bound, fpreal v,
                                        fpreal boundval);

    /// Reduce a into b, b is assumed to have our cumaltive total.
    /// If final is true, it is assumed as is also an accumulation.
    static fpreal64     applyReduction(REDUCE_NAMES op, 
                                    fpreal64 a, fpreal64 b,
                                    bool final);
    /// Similar but also interpets the mask.
    static fpreal64     applyMaskedReduction(REDUCE_NAMES op, 
                                    fpreal64 a, fpreal64 mask, fpreal64 b, 
                                    bool final);
    static fpreal64     getReductionDefault(REDUCE_NAMES op);

    THREADED_METHOD2_CONST(SIM_RawField, shouldMultiThread(), reduceOpInternal,
                    fpreal64 *, sum,
                    REDUCE_NAMES, op)
    void        reduceOpInternalPartial(fpreal64 *sum, REDUCE_NAMES op, const UT_JobInfo &info) const;

    THREADED_METHOD5_CONST(SIM_RawField, shouldMultiThread(), reduceMaskedOpInternal,
                    fpreal64 *, sum,
                    fpreal64 *, masktotal,
                    REDUCE_NAMES, op,
                    const SIM_RawField *, mask,
                    bool, maskissdf)
    void        reduceMaskedOpInternalPartial(fpreal64 *sum, fpreal64 *masktotal, REDUCE_NAMES op, const SIM_RawField *maskfield, bool maskissdf, const UT_JobInfo &info) const;

    void                sortAllVoxels(UT_FloatArray &voxelvalues) const;
    void                sortAllVoxelsMasked(UT_FloatArray &voxelvalues, 
                                            const SIM_RawField *maskfield, 
                                            bool maskissdf) const;


    SIM_FieldSample              mySample;
    UT_VoxelArrayF              *myField;

    // This is the size and offset which converts our UT_VoxelArray into
    // world coordinates.
    UT_Vector3                   myOrig, mySize;
    
    // This is our actual official size.
    UT_Vector3                   myBBoxOrig, myBBoxSize;

    UT_Vector3                   myVoxelSize;
    fpreal                       myVoxelDiameter;

    SIM_FieldBoundary            myBoundary[3][2];
    fpreal                       myBoundaryValue[3][2];
};


///
/// Allows you to iterate over the voxel cells of a RawField.  This
/// is different than iterating over the voxel samples!
/// The idea of getting or setting the raw value doesn't apply
/// here as we are not necessarily alligned with the sample grid.
///
class SIM_API SIM_RawFieldCellIterator
{
public:
             SIM_RawFieldCellIterator() 
             { 
                 myRes[0] = myRes[1] = myRes[2] = 0;
                 rewind(); 
             }
    virtual ~SIM_RawFieldCellIterator()
             {
             }
    
    void        setArray(const SIM_RawField *field)
                {
                    field->getVoxelRes(myRes[0], myRes[1], myRes[2]);
                    rewind();
                }

    void        rewind() 
                { 
                    myIdx[0] = 0; 
                    myIdx[1] = 0; 
                    myIdx[2] = 0; 
                }

    bool        atEnd() const 
                { 
                    return myIdx[2] >= myRes[2];
                }

    void        advance() 
                { 
                    myIdx[0]++;
                    if (myIdx[0] >= myRes[0])
                    {
                        myIdx[0] = 0;
                        myIdx[1]++;
                        if (myIdx[1] >= myRes[1])
                        {
                            myIdx[1] = 0;
                            myIdx[2]++;
                            if (myIdx[2] >= myRes[2])
                            {
                                // All done!
                                myIdx[2] = myRes[2];
                            }
                        }
                    }
                }

    int         x() const { return myIdx[0]; }
    int         y() const { return myIdx[1]; }
    int         z() const { return myIdx[2]; }
    int         idx(int axis) const { return myIdx[axis]; }
    
    /// Returns true if we are at the start of a new tile.
    /// Since we don't work on tile basis, we define a single
    /// x pass to be a tile.
    bool        isStartOfTile() const
                {
                    return myIdx[0] == 0;
                }

protected:
    int                                  myIdx[3];
    int                                  myRes[3];
};

//
// These macros create the CALL_VOXELPROBE macro
//
// CALL_VOXELPROBE(const SIM_RawField *src, 
//                 SIM_RawField *dst, 
//                 callFunc(dst, theprobe))
//
// callFunc should be
//
// template <typename P>
// void
// callFunc(SIM_RawField *dst, P &probe)
//
// dst->isMatching(src) must be true.
//
// The callFunc will receive in the parameter "theprobe" a
// UT_VoxelProbeAverage properly templated against the difference
// between the source and dest fields.
// theprobe.setIndex(x, y, z);
// will cause theprobe.getValue() to return the interpolated value
// in source of the destination's voxel sample (x, y, z).  This takes
// into account the different sampling patterns.
//
#define CALL_PROBEXY(XSTEP, YSTEP, src, dst, command)   \
{                                                       \
    if ((dstsample ^ srcsample) & 4)                    \
    {                                                   \
        if (dstsample & 4)      /* ZStep -1 */          \
        {                                               \
            UT_VoxelProbeAverage<fpreal32, XSTEP, YSTEP, -1> theprobe;  \
            theprobe.setArray(src->field());            \
            command;                                    \
        }                                               \
        else                    /* ZStep 1 */           \
        {                                               \
            UT_VoxelProbeAverage<fpreal32, XSTEP, YSTEP, 1> theprobe;   \
            theprobe.setArray(src->field());            \
            command;                                    \
        }                                               \
    }                                                   \
    else                        /* ZStep 0 */           \
    {                                                   \
        UT_VoxelProbeAverage<fpreal32, XSTEP, YSTEP, 0> theprobe;       \
        theprobe.setArray(src->field());                \
        command;                                        \
    }                                                   \
}

#define CALL_PROBEX(XSTEP, src, dst, command)           \
{                                                       \
    if ((dstsample ^ srcsample) & 2)                    \
    {                                                   \
        if (dstsample & 2)      /* YStep -1 */          \
        {                                               \
            CALL_PROBEXY(XSTEP, -1, src, dst, command); \
        }                                               \
        else                    /* YStep 1 */           \
        {                                               \
            CALL_PROBEXY(XSTEP, 1, src, dst, command);  \
        }                                               \
    }                                                   \
    else                        /* YStep 0 */           \
    {                                                   \
        CALL_PROBEXY(XSTEP, 0, src, dst, command);      \
    }                                                   \
}

#define CALL_VOXELPROBE(src, dst, command)              \
{                                                       \
    int         dstsample, srcsample;                   \
    dstsample = dst->getSample();                       \
    srcsample = src->getSample();                       \
    if ((dstsample ^ srcsample) & 1)                    \
    {                                                   \
        if (dstsample & 1)      /* XStep -1 */          \
        {                                               \
            CALL_PROBEX(-1, src, dst, command);         \
        }                                               \
        else                    /* XStep 1 */           \
        {                                               \
            CALL_PROBEX(1, src, dst, command);          \
        }                                               \
    }                                                   \
    else                        /* XStep 0 */           \
    {                                                   \
        CALL_PROBEX(0, src, dst, command);              \
    }                                                   \
}

#endif

---