GAS/GAS_SPH.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:        GAS_SPH.h ( GSA Library, C++)
 *
 * COMMENTS:
 *      Defines a spherical particle hydronymics search structure.
 *      This provides a toolbox of things to do with SPH weighting functions.
 */

#ifndef __GAS_SPH__
#define __GAS_SPH__

#include "GAS_API.h"

#include <GEO/GEO_AttributeHandle.h>

#include <UT/UT_FloatArray.h>
#include <UT/UT_IntArray.h>
#include <UT/UT_PtrArray.h>
#include <UT/UT_RefArray.h>
#include <UT/UT_Vector3Array.h>

class GAS_SubSolver;
class GU_Detail;
class GB_PointGroup;
class GEO_Point;
class GEO_PointTree;
class SIM_Engine;
class UT_HashTable;
class UT_FloatArray;
class UT_String;

class GAS_API GAS_SPH
{
public:
         GAS_SPH();
        ~GAS_SPH();

    enum sphWeight
    {
        POLY6,
        SPIKY,
        VISCOSITY,
        SPIKY2,
        SPIKY5,
        POLY6_UNSCALED,
        SPIKY_UNSCALED
    };

    // A static function which builds an appropriate warning message
    // for missing particle fluid attributes
    static void missingAttribWarning(SIM_Engine &engine, GAS_SubSolver *solver,
                                const UT_String &objName,
                                const UT_String &geometryName,
                                const UT_String &attribName);

    // Returns false if initalization failed, due to missing attributes,
    // etc.  If fullinit == true, then we initialize cells by adding the
    // indices of all points in the cell and all points in neighbour cells.
    // If not, then on each query we have to examine multiple neighbour cells.
    // radiusScale is used to force an increase of the radius size used for
    // density queries and such.  needVolume indicates whether or not mass
    // and density attributes are required.  They are not need, for intance,
    // when surfacing a particle fluid.
    bool        initialize(const GU_Detail *gdp, bool fullinit = true,
                        const GB_PointGroup *excludegroup = NULL,
                        fpreal radiusScale = 1.0, bool needVolume = true);


    // Perform a simple sampling operation.
    // The color variable, if present, will be set to the total weight
    // of the kernels without taking the attribute into account.
    // This is useful for efficiently calculating (pi+pj)/2, for example,
    // as the constant term can be pulled out and multiplied by the color.
    fpreal      sampleF(const UT_Vector3 &p, 
                        GEO_AttributeHandle gah,
                        sphWeight wfunc,
                        fpreal *color = 0,
                        bool normalize = false,
                        bool volumeScale = true);
    fpreal      sampleF(const UT_Vector3 &p,
                        UT_FloatArray &floats,
                        sphWeight wfunc,
                        bool volumescale = true,
                        bool normalize = false);
    UT_Vector3  sampleV3(const UT_Vector3 &p, 
                        GEO_AttributeHandle gah,
                        sphWeight wfunc,
                        fpreal *color = 0,
                        bool normalize = false);

    // A threadsafe velocity sampling function.  This assumes that
    // initVelocityData has been invoked on this sph object.
    // excludePoint is a point to exclude from the velocity
    // calculation (in case we want to sample the velocity around
    // a point, but ignore the contribution of that point).
    UT_Vector3  sampleVelocity(const UT_Vector3 &p,
                        sphWeight wfunc,
                        int excludePoint = -1);

    // Find the XSPH velocity for the given point.  We take "velocities"
    // in an input array in order to support a more general notion
    // of velocity.
    UT_Vector3 xsphVelocity(int ptnum, sphWeight wfunc,
                        const UT_Vector3Array &velocities,
                        fpreal xsph_constant);

    // Calculates the density field at a point.  This is special
    // because it ignores particle density and just sums the mass
    // weighting...
    fpreal      sampleDensity(const UT_Vector3 &p,
                        sphWeight wfunc,
                        const bool massScale = true);

    // A convenience function which samples two densities (using the
    // two given kernels) and also populates the given array with
    // a list of points close to the given point.
    void         sampleDoubleDensity(const UT_Vector3 &p,
                        sphWeight wfunc1, sphWeight wfunc2,
                        fpreal &density1, fpreal &density2,
                        UT_PtrArray<GEO_Point *> &ptlist,
                        const bool massScale = true);

    // Calculate the surface density at a point using the given
    // surface distance attribute.  This uses the surface definition
    // from the "Adaptively Sampled Particle Fluids" paper.
    fpreal      sampleSurfaceDensity(const UT_Vector3 &p,
                        GEO_AttributeHandle &distance_gah);

    // This function is similar to the one above, but uses a fixed
    // distance from each particle.  If distanceModifier is true then
    // the given surface distance is just treated as a multiplier of
    // each particles interaction radius.
    fpreal      sampleSurfaceDensity(const UT_Vector3 &p,
                        fpreal surfaceDistance,
                        const bool distanceModifier = false);

    // Calculate the gradient of the surface density field at
    // the given point.
    UT_Vector3  sampleSurfaceDensityGrad(const UT_Vector3 &p,
                        GEO_AttributeHandle &distance_gah);

    // This function is similar to the one above, but since a fixed
    // distance from each particle is used, the vector ends up being
    // somewhat different.
    UT_Vector3  sampleSurfaceDensityGrad(const UT_Vector3 &p);

    // Computes the color field - 1 inside the area handled by
    // particles, 0 outside.  This uses an implied attribute that
    // is 1 wherever particles are and zero everywhere else.
    fpreal      sampleColor(const UT_Vector3 &p,
                        sphWeight wfunc);

    // Functions to compute the gradient and laplacian of the
    // color field.  Used in the calculation of surface tension
    // forces.
    UT_Vector3  colorGradient(const UT_Vector3 &p,
                        sphWeight wfunc);

    // Computes an unscaled color gradient (eg. no scaling by
    // particle size, density, mass, etc.)
    UT_Vector3  colorGradientUnscaled(const UT_Vector3 &p,
                        bool scale = false);

    fpreal      colorLaplacian(const UT_Vector3 &p,
                        sphWeight wfunc,
                        GEO_AttributeHandle *scale_gah = NULL,
                        fpreal point_scale = 1.0);

    UT_Vector3  gradient(const UT_Vector3 &p,
                        GEO_AttributeHandle gah,
                        sphWeight wfunc,
                        UT_Vector3 *color = 0,
                        bool debugprint = false,
                        GEO_AttributeHandle *scale_gah = NULL,
                        fpreal point_scale = 1.0);

    // This function takes the gradient of a Vector3 attribute
    // along the given axis.
    UT_Vector3  gradient(const UT_Vector3 &p,
                        GEO_AttributeHandle gah,
                        sphWeight wfunc,
                        int axis, UT_Vector3 *color = 0,
                        bool normalize = false);

    // Samples the laplacian of a float[1] attribute.
    fpreal      laplacian(const UT_Vector3 &p,
                        GEO_AttributeHandle gah,
                        sphWeight wfunc,
                        fpreal *color = 0);

    // Samples the laplacian of a float[n] attribute
    // where 1 <= n <= 4.
    UT_FloatArray laplacianF(const UT_Vector3 &p,
                        GEO_AttributeHandle gah,
                        sphWeight wfunc, int size,
                        fpreal *color = 0);

    UT_Vector3  laplacianV3(const UT_Vector3 &p,
                        GEO_AttributeHandle gah,
                        sphWeight wfunc,
                        fpreal *color = 0,
                        GEO_AttributeHandle *scale_gah = NULL,
                        fpreal point_scale = 1.0);

    // Returns all points close to the given point.
    void        findClosePoints(const UT_Vector3 &p,
                        UT_PtrArray<GEO_Point *> &ptlist);

    // Returns the maximum number of elements in a cell.
    int         getMaxElements() const { return myMaxElements; }

    // Initializes information to allow fast velocity queries to
    // be made to this structure (without having to pass an
    // attribute handle in).  Also has the option of initializing
    // a KD tree of the points in the structure's gdp.  The purpose
    // of this is so that if no close points can be found during
    // a velocity query, the velocity of the nearest neighbour can
    // be taken instead.
    bool        initVelocityData(bool initKD = false);

    // Get the width of a particle cell.
    fpreal      getCellWidth() { return myCellWidth; }

    // This function instructs the gdp to exclude all points
    // with fewer than numNeighbours neighbours from all further
    // calculations.  It also returns an array of boolean values,
    // indicating what points were excluded.
    const UT_RefArray<bool> &
                excludeLonePoints(bool exclude, int numNeighbours = 2);

protected:
    // A class that defines a cubic grid cell.
    class gas_SPHGridCell {
        public:
            gas_SPHGridCell() {}
            gas_SPHGridCell(int i, int j, int k, int index)
                    : myIndices(15), myXindex(i), myYindex(j),
                      myZindex(k), myIndex(index), myNumEntries(0)
            {
            }
            virtual ~gas_SPHGridCell()
            {
            }

            UT_IntArray  myIndices;

            int          myXindex;
            int          myYindex;
            int          myZindex;

            int          myIndex;

            int          myNumEntries;
    };

    typedef UT_RefArray<gas_SPHGridCell>     gas_CellArray;

    // These take the radius squared and determine the weight 
    // or gradient of the weight for the given function.
    // h is the radius of the kernel to sample.
    fpreal64    weight(fpreal64 r2, fpreal64 h, sphWeight wfunc);
    fpreal64    dweight(fpreal64 r, fpreal64 h, sphWeight wfunc);
    fpreal64    d2weight(fpreal64 r, fpreal64 h, sphWeight wfunc);
    fpreal64    scaleddweight(fpreal64 r, fpreal64 h, sphWeight wfunc);

    // Places a point in the appropriate grid cell.
    void                 assignPoint(const GEO_Point *ppt);

    // Helper function used to get a grid cell from the hash
    // table, or generate a new one if necessary.
    gas_SPHGridCell *    getCell(const UT_Vector3 &p, bool create = false);
    gas_SPHGridCell *    getCell(int i, int j, int k, bool create = false);

    // Store the maximum number of cells in any dimension,
    // along with the minimum x, y and z coordinates.
    int                  myMaxCells;
    int                  myMaxCells2;
    fpreal               myXmin, myYmin, myZmin;
    fpreal               myCellWidth;
    fpreal               myCellWidth2;

    // Stores the maximum number of elements in a cell.
    int                  myMaxElements;

    const GU_Detail     *myGdp;

    UT_HashTable        *myGridSet;
    gas_CellArray        myCellArray;

    bool                 myFullInit;

    GEO_AttributeHandle  myRadGah, myMassGah, myDensityGah;

    // We can optionally store a list of the mass, density and
    // velocity at each point in our geometry.  We do this to
    // allow multithreaded sampling of the same sph structure.
    UT_FloatArray        myMasses;
    UT_FloatArray        myDensities;
    UT_FloatArray        myRadii;
    UT_Vector3Array      myVelocities;

    // Also store a point tree to be used for velocity evaluations.
    // We need this so that a velocity can be computed regardless of
    // the given point in space.
    GEO_PointTree       *myPointTree;

    // Radius scale used for weight queries
    fpreal               myRadiusScale;

    // A set of flags indicating whether or not points should
    // be excluded from all calculations.
    UT_RefArray<bool>    myExcludePoint;
};

#endif

---