SIM/SIM_SolverSNOW.C

/*
 * Copyright (c) 2012
 *      Side Effects Software Inc.  All rights reserved.
 *
 * Redistribution and use of Houdini Development Kit samples in source and
 * binary forms, with or without modification, are permitted provided that the
 * following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 * 2. The name of Side Effects Software may not be used to endorse or
 *    promote products derived from this software without specific prior
 *    written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY SIDE EFFECTS SOFTWARE `AS IS' AND ANY EXPRESS
 * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN
 * NO EVENT SHALL SIDE EFFECTS SOFTWARE BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 *----------------------------------------------------------------------------
 */

#include "SIM_SolverSNOW.h"
#include <strstream.h>
#include <UT/UT_DSOVersion.h>
#include <UT/UT_Floor.h>
#include <UT/UT_Vector3.h>
#include <UT/UT_WorkBuffer.h>
#include <GU/GU_Detail.h>
#include <GU/GU_RayIntersect.h>
#include <GU/GU_PrimTriStrip.h>
#include <PRM/PRM_Include.h>
#include <SIM/SIM_Engine.h>
#include <SIM/SIM_Options.h>
#include <SIM/SIM_Object.h>
#include <SIM/SIM_ObjectArray.h>
#include <SIM/SIM_DopDescription.h>
#include <SIM/SIM_Random.h>
#include <SIM/SIM_RandomTwister.h>
#include <SIM/SIM_Position.h>

#define VOXELARRAY_XDIVISIONS   "xdivisions"
#define VOXELARRAY_YDIVISIONS   "ydivisions"
#define VOXELARRAY_ZDIVISIONS   "zdivisions"

static PRM_Name  theBirthRateName(SIM_NAME_BIRTHRATE, "Birth Rate");
static PRM_Name  theDivisionsName(SIM_NAME_DIVISIONS, "Divisions");
static PRM_Name  theOriginalDepthName(SIM_NAME_ORIGINALDEPTH, "Original Depth");

using namespace HDK_Sample;

void
initializeSIM(void *)
{
    IMPLEMENT_DATAFACTORY(SIM_SolverSNOW);
    IMPLEMENT_DATAFACTORY(SNOW_VoxelArray);
}

SNOW_VoxelArray::SNOW_VoxelArray(const SIM_DataFactory *factory)
    : SIM_Geometry(factory),
      myVoxelArray(0)
{
}

SNOW_VoxelArray::~SNOW_VoxelArray()
{
    freeArray();
}

int
SNOW_VoxelArray::getXDivisions() const
{
    return myXDivisions;
}

void
SNOW_VoxelArray::setXDivisions(int divisions)
{
    if( divisions != myXDivisions )
    {
        freeArray();
        myXDivisions = divisions;
    }
}

int
SNOW_VoxelArray::getYDivisions() const
{
    return myYDivisions;
}

void
SNOW_VoxelArray::setYDivisions(int divisions)
{
    if( divisions != myYDivisions )
    {
        freeArray();
        myYDivisions = divisions;
    }
}

int
SNOW_VoxelArray::getZDivisions() const
{
    return myZDivisions;
}

void
SNOW_VoxelArray::setZDivisions(int divisions)
{
    if( divisions != myZDivisions )
    {
        freeArray();
        myZDivisions = divisions;
    }
}

u8
SNOW_VoxelArray::getVoxel(int x, int y, int z) const
{
    if( !myVoxelArray )
        allocateArray();

    return myVoxelArray->getValue(x, y, z);
}

void
SNOW_VoxelArray::setVoxel(u8 voxel, int x, int y, int z)
{
    if( !myVoxelArray )
        allocateArray();

    if (x < 0 || x >= myXDivisions)
        return;
    if (y < 0 || y >= myYDivisions)
        return;
    if (z < 0 || z >= myZDivisions)
        return;

    myVoxelArray->setValue(x, y, z, voxel);
}

GU_ConstDetailHandle
SNOW_VoxelArray::getGeometrySubclass() const
{
    ((SNOW_VoxelArray *)this)->buildGeometryFromArray();
    return myDetailHandle;
}

void
SNOW_VoxelArray::freeArray() const
{
    delete myVoxelArray;
    myVoxelArray = 0;
}

void
SNOW_VoxelArray::allocateArray() const
{
    UT_ASSERT(myVoxelArray == 0);

    myVoxelArray = new UT_VoxelArray<u8>;
    myVoxelArray->size(myXDivisions, myYDivisions, myZDivisions);

    // We want out of bound values to evaluate to wall voxels.
    myVoxelArray->setBorder(UT_VOXELBORDER_CONSTANT, VOXEL_WALL);
}

GEO_Point *
SNOW_VoxelArray::createOrFindPoint(GU_Detail *gdp, int x, int y, int z)
{
    int         idx;
    UT_Thing    thing;
    GEO_Point   *pt;

    idx = (z * (myYDivisions + 1) + y)*(myXDivisions + 1) + x;

    UT_Hash_Int hash(idx);

    if (myPointHash.findSymbol(hash, &thing))
    {
        return (GEO_Point *) thing.value.voidp;
    }

    // Create and add!

    pt = gdp->appendPoint();

    UT_Vector3          v4;
    v4.assign((fpreal) x / (fpreal) (myXDivisions + 1), 
              (fpreal) y / (fpreal) (myYDivisions + 1), 
              (fpreal) z / (fpreal) (myZDivisions + 1));
    pt->getPos() = v4;

    // And add to the hash...
    thing.value.voidp = pt;
    myPointHash.addSymbol(hash, thing);

    return pt;
}

void             
SNOW_VoxelArray::buildFace(GU_Detail *gdp,
                           int x0, int y0, int z0,
                           int x1, int y1, int z1,
                           int x2, int y2, int z2,
                           int x3, int y3, int z3)
{
    // Triangle strip code.
    GU_PrimTriStrip     *tristrip;
    GEO_Point           *pt;

    // Do not append points, as we may scavenge.
    tristrip = GU_PrimTriStrip::build(gdp, 4, 0);

    pt = createOrFindPoint(gdp, x0, y0, z0);
    (*tristrip)(0).setPt(pt);
    pt = createOrFindPoint(gdp, x1, y1, z1);
    (*tristrip)(1).setPt(pt);
    pt = createOrFindPoint(gdp, x3, y3, z3);
    (*tristrip)(2).setPt(pt);
    pt = createOrFindPoint(gdp, x2, y2, z2);
    (*tristrip)(3).setPt(pt);
}

void
SNOW_VoxelArray::buildGeometryFromArray()
{
    if( myDetailHandle.isNull() )
    {
        GU_Detail       *gdp = new GU_Detail();
        int              x, xdiv, y, ydiv, z, zdiv;
        int              xstep, ystep, zstep;

        myDetailHandle.allocateAndSet(gdp);
        xdiv = getXDivisions();
        ydiv = getYDivisions();
        zdiv = getZDivisions();

        // Find the appropriate step value...
        xstep = ystep = zstep = 1;
        if (xdiv > 64)
            xstep = xdiv / 64;
        if (ydiv > 64)
            ystep = ydiv / 64;
        if (zdiv > 64)
            zstep = zdiv / 64;
        
        for( z = 0; z < zdiv; z+=zstep )
        {
            for( y = 0; y < ydiv; y+=ystep )
            {
                for( x = 0; x < xdiv; x+=xstep )
                {
                    if (getVoxel(x, y, z) == VOXEL_SNOW)
                    {
                        // Check each of the cardinal directions
                        // to see if we want to build a face.

                        // We want to render the faces of this cube
                        // that are bordered by an empty unit.
                        // We specify the points as (x,y,z) triplets.
                        // This cube is (x,y,z) to (x+1,y+1,z+1)
                        if (getVoxel(x-xstep, y, z) != VOXEL_SNOW)
                        {
                            buildFace(  gdp, x, y, z,
                                        x, y+ystep, z,
                                        x, y+ystep, z+zstep,
                                        x, y, z+zstep );
                        }
                        if (getVoxel(x+xstep, y, z) != VOXEL_SNOW)
                        {
                            buildFace(  gdp, x+xstep, y, z,
                                        x+xstep, y, z+zstep,
                                        x+xstep, y+ystep, z+zstep,
                                        x+xstep, y+ystep, z );
                        }
                        if (getVoxel(x, y-ystep, z) != VOXEL_SNOW)
                        {
                            buildFace(  gdp, x, y, z,
                                        x, y, z+zstep,
                                        x+xstep, y, z+zstep,
                                        x+xstep, y, z );
                        }
                        if (getVoxel(x, y+ystep, z) != VOXEL_SNOW)
                        {
                            buildFace(  gdp, x, y+ystep, z,
                                        x+xstep, y+ystep, z,
                                        x+xstep, y+ystep, z+zstep,
                                        x, y+ystep, z+zstep );
                        }
                        if (getVoxel(x, y, z-zstep) != VOXEL_SNOW)
                        {
                            buildFace(  gdp, x, y, z,
                                        x+xstep, y, z,
                                        x+xstep, y+ystep, z,
                                        x, y+ystep, z );
                        }
                        if (getVoxel(x, y, z+zstep) != VOXEL_SNOW)
                        {
                            buildFace(  gdp, x, y, z+zstep,
                                        x, y+ystep, z+zstep,
                                        x+xstep, y+ystep, z+zstep,
                                        x+xstep, y, z+zstep );
                        }
                    }
                }
            }
        }

        // Wipe out all the points we allocated.
        myPointHash.clear();
    }
}

void
SNOW_VoxelArray::initializeSubclass()
{
    SIM_Geometry::initializeSubclass();
    freeArray();
    myXDivisions = 0;
    myYDivisions = 0;
    myZDivisions = 0;
    myDetailHandle.clear();
}

void
SNOW_VoxelArray::makeEqualSubclass(const SIM_Data *source)
{
    const SNOW_VoxelArray       *srcvox;

    SIM_Geometry::makeEqualSubclass(source);
    srcvox = SIM_DATA_CASTCONST(source, SNOW_VoxelArray);
    if( srcvox )
    {
        setXDivisions(srcvox->getXDivisions());
        setYDivisions(srcvox->getYDivisions());
        setZDivisions(srcvox->getZDivisions());
        
        if (srcvox->myVoxelArray)
        {
            // Copy over the voxels.
            allocateArray();

            *myVoxelArray  = *srcvox->myVoxelArray;
        }
        else
        {
            // No voxel array, so nothing to copy.
            freeArray();
        }
    }
}

void
SNOW_VoxelArray::saveSubclass(ostream &os) const
{
    SIM_Options          voxeldata;
    int                  x, xdiv, y, ydiv, z, zdiv, value;

    SIM_Geometry::saveSubclass(os);

    xdiv = getXDivisions();
    ydiv = getYDivisions();
    zdiv = getZDivisions();
    voxeldata.setOptionI(VOXELARRAY_XDIVISIONS, xdiv);
    voxeldata.setOptionI(VOXELARRAY_YDIVISIONS, ydiv);
    voxeldata.setOptionI(VOXELARRAY_ZDIVISIONS, zdiv);
    saveOptionPacket(os, classname(), &voxeldata);

    os << "{" << endl;
    for( z = 0; z < zdiv; z++ )
    {
        for( y = 0; y < ydiv; y++ )
        {
            os << "\t";
            for( x = 0; x < xdiv; x++ )
            {
                value = getVoxel(x, y, z);
                os << " " << value;
            }
            os << endl;
        }
    }
    os << "}" << endl;
}

bool
SNOW_VoxelArray::loadSubclass(UT_IStream &is)
{
    SIM_Options          voxeldata;
    UT_WorkBuffer        buf;
    int                  xdiv, ydiv, zdiv, idx, arraysize;
    int                  x, y, z;
    int                  value;
    bool                 result = true;

    if (!SIM_Geometry::loadSubclass(is))
        return false;

    if( loadOptionPacket(is, classname(), &voxeldata) )
    {
        xdiv = voxeldata.getOptionI(VOXELARRAY_XDIVISIONS);
        ydiv = voxeldata.getOptionI(VOXELARRAY_YDIVISIONS);
        zdiv = voxeldata.getOptionI(VOXELARRAY_ZDIVISIONS);
        setXDivisions(xdiv);
        setYDivisions(ydiv);
        setZDivisions(zdiv);
        arraysize = xdiv * ydiv * zdiv;
        idx = 0;
        x = y = z = 0;
        if( is.getLine(buf) && *buf.buffer() == '{' )
        {
            while( is.getLine(buf) && *buf.buffer() != '}' )
            {
                istrstream       bufis(buf.lock(), buf.length());

                while( idx < arraysize && bufis )
                {
                    if( (bufis >> value) )
                    {
                        setVoxel(value, x, y, z);
                        x++;
                        if (x >= xdiv);
                        {
                            x = 0;
                            y++;
                            if (y >= ydiv)
                            {
                                y = 0;
                                z++;
                            }
                        }
                    }
                }
                buf.release();
            }
            UT_ASSERT(idx == arraysize);
        }
    }
    else
        result = false;

    return result;
}

void
SNOW_VoxelArray::handleModificationSubclass(int code)
{
    SIM_Geometry::handleModificationSubclass(code);
    myDetailHandle.clear();
}

void
SNOW_VoxelArray::collapseAllTiles()
{
    if (!myVoxelArray)
        return;

    myVoxelArray->collapseAllTiles();
}

SIM_SolverSNOW::SIM_SolverSNOW(const SIM_DataFactory *factory)
    : BaseClass(factory),
      SIM_OptionsUser(this)
{
}

SIM_SolverSNOW::~SIM_SolverSNOW()
{
}

const SIM_DopDescription *
SIM_SolverSNOW::getSolverSNOWDopDescription()
{
    static PRM_Template  theTemplates[] = {
        PRM_Template(PRM_FLT_J,         1, &theBirthRateName, PRMpointOneDefaults),
        PRM_Template(PRM_INT_J,         1, &theDivisionsName, PRMtenDefaults),
        PRM_Template(PRM_INT_J,         1, &theOriginalDepthName),
        PRM_Template()
    };

    static SIM_DopDescription    theDopDescription(true,
                                                   "hdk_snowsolver",
                                                   "SNOW Solver",
                                                   SIM_SOLVER_DATANAME,
                                                   classname(),
                                                   theTemplates);

    return &theDopDescription;
}

SIM_Random *
SIM_SolverSNOW::createRandomData(SIM_Object *obj) const
{
    SIM_Random  *rand = 0;
    
    // Create the random data as subdata attached to the solver. First
    // we look for any existing SIM_Random. If none is found, we create
    // a SIM_RandomTwister.
    rand = SIM_DATA_GET(*obj, "Random", SIM_Random);
    if( !rand )
        rand = SIM_DATA_CREATE(*obj, "Random", SIM_RandomTwister, 0);

    return rand;
}

bool
SIM_SolverSNOW::brownianize(int &v, int dv, int max, SIM_Random *rand) const
{
    if (dv)
    {
        v += dv;
        if (v < 0 || v >= max)
            return false;
    }
    else
    {
        v += rand_choice(3, rand) - 1;
        if (v < 0)
            v = 0;
        if (v >= max)
            v = max - 1;
    }

    return true;
}

int
SIM_SolverSNOW::clearInDirection(const SNOW_VoxelArray &snow,
                        int sx, int sy, int sz,
                        int dx, int dy, int dz,
                        int &rx, int &ry, int &rz,
                        int maxdist,
                        SIM_Random *rand) const
{
    int                 dist = 0;
    int                 xdiv, ydiv, zdiv;

    xdiv = snow.getXDivisions();
    ydiv = snow.getYDivisions();
    zdiv = snow.getZDivisions();

    while (1)
    {
        // Assume our current location is invalid.  Try one step.
        if (!brownianize(sx, dx, xdiv, rand))
            return maxdist;
        if (!brownianize(sy, dy, ydiv, rand))
            return maxdist;
        if (!brownianize(sz, dz, zdiv, rand))
            return maxdist;
            
        // Now, see if we are suddenly valid.
        if (snow.getVoxel(sx, sy, sz) == VOXEL_EMPTY)
        {
            break;
        }

        dist++;
    }

    rx = sx;
    ry = sy;
    rz = sz;
    return dist;
}

// Assumes that x,y,z had a piece of snow.  Will find a new home
// for this snow.
void
SIM_SolverSNOW::clearSnow(SNOW_VoxelArray &snow,
                          int x, int y, int z, SIM_Random *rand) const
{
    // There are 6 primary axes for snow to be distributed along.
    // The snow tries all 6 directions, and moves in teh one that
    // has the shortest path.

    int                 end_x[6], end_y[6], end_z[6], dist[6];

    int         dxvals[6] = { -1,  0,  0,  1,  0,  0 };
    int         dyvals[6] = {  0, -1,  1,  0,  0,  0 };
    int         dzvals[6] = {  0,  0,  0,  0,  1, -1 };
    int         direction, mindir, mindist = 320000;

    for (direction = 0; direction < 6; direction++)
    {
        dist[direction] =
            clearInDirection(snow,
                        x, y, z, 
                        dxvals[direction], dyvals[direction], dzvals[direction],
                        end_x[direction], end_y[direction], end_z[direction],
                        mindist,
                        rand);

        if (dist[direction] < mindist)
        {
            mindir = direction;
            mindist = dist[direction];
        }
    }

    if (mindist == 320000)
    {
        // Complete failure!
        UT_ASSERT(!"No snow removal possible!");
    }
    else
    {
        // Store in the resulting position...
        snow.setVoxel(VOXEL_SNOW, end_x[mindir], end_y[mindir], end_z[mindir]);
    }
}

int
SIM_SolverSNOW::rand_choice(int numchoice, SIM_Random *rand) const
{
    int         choice;
    
    choice = rand->choice(numchoice);

    return choice;
}

void
SIM_SolverSNOW::fillRow(SNOW_VoxelArray &snow,
                        fpreal startx, fpreal endx,
                        int y, int z,
                        u8 voxeltype,
                        SIM_Random *rand) const
{
    int         xdiv, ydiv, zdiv;
    int         x, sx, ex;

    xdiv = snow.getXDivisions();
    ydiv = snow.getYDivisions();
    zdiv = snow.getZDivisions();

    sx = (int)(startx * xdiv);
    ex = (int)(endx * xdiv);
    if (sx < 0) sx = 0;
    if (sx >= xdiv) return;
    if (ex < 0) return;
    if (ex >= xdiv) ex = xdiv - 1;

    if (voxeltype == VOXEL_OBJECT)
    {
        for (x = sx; x < ex; x++)
        {
            // Set this voxel.  TODO: Move away snow!
            if (snow.getVoxel(x, y, z) == VOXEL_SNOW)
                clearSnow(snow, x, y, z, rand);
            snow.setVoxel(VOXEL_OBJECT, x, y, z);
        }
    }
    else if (voxeltype == VOXEL_SNOW)
    {
        for (x = sx; x < ex; x++)
        {
            snow.setVoxel(VOXEL_SNOW, x, y, z);
        }
    }
}

void
SIM_SolverSNOW::applyGeometry(SNOW_VoxelArray &snow,
                              const GU_ConstDetailHandle &gdh,
                              const UT_DMatrix4 &xform,
                              u8 voxeltype,
                              SIM_Random *rand) const
{
    if( !gdh.isNull() )
    {
        GU_DetailHandleAutoReadLock      gdl(gdh);
        const GU_Detail                 *gdp = gdl.getGdp();
        int                              xdiv, ydiv, zdiv;
        int                              y, z;

        xdiv = snow.getXDivisions();
        ydiv = snow.getYDivisions();
        zdiv = snow.getZDivisions();

        GU_RayIntersect         *isect;
        UT_BoundingBox           bbox;
        GU_RayInfo                       hitinfo;
        int                              numhit, hitnum;
        UT_Matrix4               fxform;

        fxform = xform;

        fxform.invert();
        gdp->getBBox(&bbox);
        bbox.transform(fxform);

        // Find where this is at least valid...
        int                      bminx, bminy, bminz;
        int                      bmaxx, bmaxy, bmaxz;
        UT_Vector3               orig, dir, pos;
        UT_Vector3               xorig, xdir, xpos;

        // Find the range of the bounding box - we only need
        // to search this part of the voxel array.
        bminx = (int)UTfloor(bbox(0, 0) * (xdiv + 1));
        if (bminx < 0) bminx = 0;
        bmaxx = (int)UTceil(bbox(0, 1) * (xdiv + 1));
        if (bmaxx >= xdiv) bmaxx = xdiv-1;
        bminy = (int)UTfloor(bbox(1, 0) * (ydiv + 1));
        if (bminy < 0) bminy = 0;
        bmaxy = (int)UTceil(bbox(1, 1) * (ydiv + 1));
        if (bmaxy >= ydiv) bmaxy = ydiv-1;
        bminz = (int)UTfloor(bbox(2, 0) * (zdiv + 1));
        if (bminz < 0) bminz = 0;
        bmaxz = (int)UTceil(bbox(2, 1) * (zdiv + 1));
        if (bmaxz >= zdiv) bmaxz = zdiv-1;

        // Build the ray intersect cache.
        isect = new GU_RayIntersect(gdp);
        
        // We build downwards so snow tends to
        // compact.
        orig.x() = 0.0;
        dir.assign(1.0, 0.0, 0.0);
        xdir = dir;
        xdir.multiply3(xform);
                
        for (z = bmaxz; z >= bminz; z--)
        {
            orig.z() = (z + 0.5) / (zdiv + 1);
            for (y = bminy; y <= bmaxy; y++)
            {
                orig.y() = (y + 0.5) / (ydiv + 1);
                hitinfo.reset();

                xorig = orig;
                xorig *= xform;

                hitinfo.init(1.0, 0.0, GU_FIND_ALL, 1e-4);

                numhit = isect->sendRay(xorig, xdir, hitinfo);

                // -1 means interrupt from user.
                if (numhit < 0)
                    return;

                // Even if there were no hits, we may still be entirely
                // inside the object.
                numhit = hitinfo.myHitList->entries();

                fpreal          lt, t;
                
                // Now, walk through each hit...
                // First "hit" occurs at position zero.  Last "hit"
                // occurs at position 1.
                lt = 0.0;
                
                for (hitnum = 0; hitnum <= numhit; hitnum++)
                {
                    if (hitnum < numhit)
                        t = (*hitinfo.myHitList)(hitnum).t;
                    else
                        t = 1.0;

                    // Determine if the lt - t segment is inside or not.
                    pos = orig;
                    pos.x() = (t + lt) / 2.0;
                    xpos = pos;
                    xpos *= xform;
                    if (isect->isInsideWinding(xpos, 0))
                    {
                        fillRow(snow, lt, t, y, z, voxeltype, rand);
                    }

                    lt = t;
                }
            }
        }

        delete isect;
    }
}

void
SIM_SolverSNOW::solveForObject(SIM_Object &object,
                              SNOW_VoxelArray &snow,
                              const SIM_Time & /*timestep*/) const
{
    // Birth new snow at top.
    int                  xdiv, ydiv, zdiv;
    int                  x, y, z;
    int                  i, n;
    fpreal               birthrate;
    SIM_Random          *rand = createRandomData(&object);

    xdiv = getDivisions();
    ydiv = getDivisions();
    zdiv = getDivisions();

    birthrate = getBirthRate();

    // Update according to the possibly changed intersection information.
    const SIM_Geometry  *geometry = 0;

    // First, clear out all old intersection information.
    for (z = 0; z < zdiv; z++)
    {
        for (y = 0; y <= ydiv; y++)
        {
            for (x = 0; x <= xdiv; x++)
            {
                if (snow.getVoxel(x, y, z) == VOXEL_OBJECT)
                    snow.setVoxel(VOXEL_EMPTY, x, y, z);
            }
        }
    }

    // Run through each affector looking for source generators...
    SIM_ObjectArray             sourceaffectors;
    SIM_ColliderInfoArray       colliderinfo;
    UT_String                   sourceobjects;

    object.getAffectors(sourceaffectors, "SIM_RelationshipSource");
    n = sourceaffectors.entries();
    for (i = 0; i < n; i++)
    {
        const SIM_Object        *affector = sourceaffectors(i);
        const SIM_Position      *pos = affector->getPosition();
        UT_DMatrix4              xform, worldtogeo;

        geometry = affector->getGeometry();
        // Ignore people that don't have a "geometry" field.
        if (!geometry)
            continue;

        geometry->getTransform(xform);
        xform.invert();
        if (pos)
        {
            pos->getInverseTransform(worldtogeo);
            xform = worldtogeo * xform;
        }

        applyGeometry(snow, geometry->getGeometry(), xform, VOXEL_SNOW, rand);
    }

    // Run through each affector looking for geometry data...
    object.getColliderInfo(colliderinfo);
    n = colliderinfo.entries();
    for (i = 0; i < n; i++)
    {
        const SIM_Object        *affector = colliderinfo(i).getAffector();
        const SIM_Position      *pos = affector->getPosition();
        UT_DMatrix4              xform, worldtogeo;

        geometry = affector->getGeometry();
        // Ignore people that don't have a "geometry" field.
        if (!geometry)
            continue;

        geometry->getTransform(xform);
        xform.invert();
        if (pos)
        {
            pos->getInverseTransform(worldtogeo);
            xform = worldtogeo * xform;
        }

        applyGeometry(snow, geometry->getGeometry(), xform, VOXEL_OBJECT, rand);
    }

    // Birth new snow at the top of the box.
    if (!UTequalZero(birthrate))
        for (y = 0; y < ydiv; y++)
        {
            for (x = 0; x < xdiv; x++)
            {
                if (rand->frandom() < birthrate)
                {
                    snow.setVoxel(VOXEL_SNOW, x, y, zdiv-1);
                }
            }
        }

    int         dxvals[9] = { -1, -1, -1,  0,  0,  0,  1,  1,  1 };
    int         dyvals[9] = { -1,  0,  1, -1,  0,  1, -1,  0,  1 };
    int         validdxidx[9];
    int         numdxidx, dxidx;

    // And move everything down one level...
#if 1
    for (z = 1; z < zdiv; z++)
    {
        // If this snow voxel is set to 1, we want to try and move it down
        // to z-1.
        // We don't want to be too consistent with our direction or we'll
        // induce a strong bias.  Thus we reverse our loops depending
        // on z value.
        int             yend, ystart, yinc;
        int             xend, xstart, xinc;
        
        if (z & 1)
        {
            ystart = 0;
            yend = ydiv;
            yinc = 1;
            xstart = 0;
            xend = xdiv;
            xinc = 1;
        }
        else
        {
            ystart = ydiv-1;
            yend = -1;
            yinc = -1;
            xstart = xdiv-1;
            xend = -1;
            xinc = -1;
        }

        for (y = ystart; y != yend; y += yinc)
        {
            for (x = xstart; x != xend; x += xinc)
            {
                if (snow.getVoxel(x, y, z) == VOXEL_SNOW)
                {
                    // Try all dx combinations.
                    numdxidx = 0;
                    for (dxidx = 0; dxidx < 9; dxidx++)
                    {
                        if (snow.getVoxel(x + dxvals[dxidx],
                                           y + dyvals[dxidx],
                                           z-1) == VOXEL_EMPTY)
                        {
                            validdxidx[numdxidx++] = dxidx;
                        }
                    }

                    if (numdxidx)
                    {
                        dxidx = rand_choice(numdxidx, rand);

                        dxidx = validdxidx[dxidx];
                        
                        // We can successfully move...
                        snow.setVoxel(VOXEL_EMPTY, x, y, z);
                        UT_ASSERT(snow.getVoxel(x + dxvals[dxidx],
                                                y + dyvals[dxidx],
                                                z-1) == VOXEL_EMPTY);
                        snow.setVoxel(VOXEL_SNOW, x + dxvals[dxidx],
                                         y + dyvals[dxidx],
                                         z-1);
                    }
                }
            }
        }
    }
#endif

    // Now we want to auto-collapse anything that is constant.
    snow.collapseAllTiles();
    snow.pubHandleModification();
}

void
SIM_SolverSNOW::setVoxelArrayAttributes(SNOW_VoxelArray *voxelarray) const
{
    if( voxelarray )
    {
        voxelarray->setXDivisions(getDivisions());
        voxelarray->setYDivisions(getDivisions());
        voxelarray->setZDivisions(getDivisions());

        int             x, y, z, div, depth;

        div = getDivisions();

        depth = getOriginalDepth();

        for (z = 0; z < depth; z++)
        {
            for (y = 0; y < div; y++)
            {
                for (x = 0; x < div; x++)
                {
                    voxelarray->setVoxel(VOXEL_SNOW, x, y, z);
                }
            }
        }
        voxelarray->collapseAllTiles();
        voxelarray->pubHandleModification();
    }
}

SIM_Solver::SIM_Result
SIM_SolverSNOW::solveSingleObjectSubclass(SIM_Engine & /*engine*/,
                                          SIM_Object &object,
                                          SIM_ObjectArray &,
                                          const SIM_Time &timestep,
                                          bool)
{
    SNOW_VoxelArray     *snow;
    SIM_Result           result;

    result = SIM_SOLVER_FAIL;
    // First, collect (or create) all the data we need from the object.
    snow = SIM_DATA_GET(object, "SnowValue", SNOW_VoxelArray);
    if (!snow)
    {
        snow = SIM_DATA_CREATE(object, "SnowValue", SNOW_VoxelArray, 0);
        setVoxelArrayAttributes(snow);
    }

    // Once we have our data, do the actual math to update the object.
    if( snow )
    {
        solveForObject(object, *snow, timestep);
        result = SIM_SOLVER_SUCCESS;
    }

    return result;
}