UT_MGOperator.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.
 * 
 * NAME:        UT_MGOperator.h ( UT Library, C++)
 *
 * COMMENTS:
 *      This class encapsulates a tiled discrete Laplacian operator that
 *      is algebraically coarsened, allowing for a sparse Poisson solve.
 */

#ifndef __UT_MGOperator__
#define __UT_MGOperator__

#include "UT_API.h"
#include "UT_MatrixSolver.h"
#include "UT_VoxelArray.h"
#include "UT_FixedVector.h"

/// This function returns true if the given tile has at least one voxel whose
/// value exceeds 0.5.
template<typename S> UT_API
bool ut_tileHasActiveVoxel(const UT_VoxelTile<S>* tile);

template<typename T>
class UT_API UT_MGOperatorT
{
public:
    /// Once initialized, this operator is algebraically coarsened until the number
    /// of degrees of freedom does not exceed the provided threshold.
    UT_MGOperatorT(int threshold = 100) : myVoxelCountThreshold(threshold),
                                          mySmoothingPasses(2),
                                          myDiffusionPasses(1),
                                          myAmplificationFactor(1.8),
                                          myJacobiOmega(2.0 / 3),
                                          myGSOmega(4.0 / 3) {}

    /// Initializes this operator to the Laplacian matrix with the given spacing and
    /// boundary conditions. This method also takes care of setting up the coarser
    /// operators. Tiles in the operator corresponds to tiles of the reference array.
    /// Active tiles are non-constant or different from inactive.
    /// All dimensions of the reference array must be divisible by 16.
    template<typename S>
    void init(const UT_VoxelArray<S>& ref, S inactive,
              const UT_Vector3T<T>& spacing,
              const UT_Vector3I& boundariesNeg,
              const UT_Vector3I& boundariesPos);
    /// This version of initialization takes a validity map (whose resolution should
    /// match that of the finest tile map) and performs the necessary setup. Active
    /// tiles have a non-zero value in the reference voxel array.
    /// Resolution of this operator will be 16 times the resolution of the reference
    /// array.
    template<typename S>
    void init(const UT_VoxelArray<S>& ref,
              const UT_Vector3T<T>& spacing,
              const UT_Vector3I& boundariesNeg,
              const UT_Vector3I& boundariesPos);
    /// This version of initialization accepts a stencil array. Tiles in the operator
    /// are activated if the corresponding tile in the stencil array has at least one
    /// voxel whose value exceeds 0.5.
    /// All dimension of the reference array must be divisible by 16.
    template<typename S>
    void initFromStencil(const UT_VoxelArray<S>& stencil,
                         const UT_Vector3T<T>& spacing,
                         const UT_Vector3I& boundariesNeg,
                         const UT_Vector3I& boundariesPos);

    /// Applies the restriction operator from level l to level l+1.
    void mgRestrict(const UT_VectorT<T>& src, UT_VectorT<T>& dst, int l) const;
    /// Applies the prolongation operator from level l+1 to level l.
    void mgProlong(const UT_VectorT<T>& src, UT_VectorT<T>& dst, int l) const;

    /// Flattens contents of the given voxel array into the vector.
    void copyFromVoxels(const UT_VoxelArray<T>& src, UT_VectorT<T>& dst) const;
    /// Copies the contents of the flattened vector into the voxel array. If topad
    /// is true, then destination array will have an extra layer of voxels
    /// compared to the internal tile map, whose entries will be set to 0.
    void copyToVoxels(const UT_VectorT<T>& src, UT_VoxelArray<T>& dst,
                      bool topad) const;

    /// Solves the linear system using multigrid.
    /// The solution is stored in x, which will be initialized to the appropriate
    /// size. If b has incomplete tiles, x will be padded with an extra layer of
    /// voxels along every dimension.
    /// If startFMG is true, the first iteration will be full-multigrid;
    /// otherwise, only v-cycles are used.    
    int solvePoisson(fpreal64 abstol, fpreal64 reltol, int maxiter,
                     UT_VoxelArray<T>& x, const UT_VoxelArray<T>& b,
                     UT_ValArray<fpreal64>* resnorminf,
                     UT_ValArray<fpreal64>* resnorm2,
                     bool startFMG = true);

    /// Calculates b-Lx at the given level and puts the answer in d (note that b
    /// and d can be the same vector).
    void subtractApplyLaplacian(const UT_VectorT<T>& x, const UT_VectorT<T>& b,
                                UT_VectorT<T>& d, int level) const;
    /// Calculates Lx at the given level and puts the answer in b.
    void applyLaplacian(const UT_VectorT<T>& x, UT_VectorT<T>& b, int level) const;

    /// Does a full multigrid solve to get an approximation for the solution to
    /// the equation Lx=b at the given level. Note that x must be pre-initialized,
    /// but its incoming contents are irrelevant.
    /// t0 and t1 are used to perform intermediate work and must be pre-allocated.
    void fullMultigrid(UT_VectorT<T>& x, const UT_VectorT<T>& b, UT_VectorT<T>& t0,
                       UT_VectorT<T>& t1, int level);
    /// Approximates a solution for the equation Lx=b using a single V-cycle. Note
    /// that x will be re-initialized, so its incoming contents are irrelevant.
    /// t is a temporary vector that must be allocated to the same size as x.
    void vcycle(UT_VectorT<T>& x, UT_VectorT<T>& t, const UT_VectorT<T>& b,
                int level);
    /// Does a single iteration of relaxation on the solution vector at the given
    /// level. For positive levels, uses Jacobi relaxation; for negative levels,
    /// performs red-black Gauss-Seidel, with redFirst indicating which color to
    /// smooth first.
    /// For Jacobi smoothing (non-negative levels) x0 is the vector to relax, x
    /// will hold the result. These vectors must be distinct! For Gauss-Seidel
    /// smoothing (negative levels), x holds the input and output solutions, so x0
    /// is unused.
    /// If ZERO_GUESS is true, it is assumed that the incoming vector is zero.
    template<bool ZERO_GUESS = false>
    void smoothLaplacian(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
                         const UT_VectorT<T>& b, T omega, int level,
                         bool redFirst) const;
    /// Does a single iteration of homogeneous relaxation (smoothing for the system
    /// Lx=0). Uses Jacobi relaxation without damping.
    /// x0 is the vector to relax, x will hold the results; these vectors must be
    /// distinct!
    void smoothLaplacianHomogeneous(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
                                    int level) const;
    /// Performs a direct solve at the coarsest level.
    void directSolve(UT_VectorT<T>& x, const UT_VectorT<T>& b) const;

    /// Returns the tile map's resolution at the given level.
    UT_Vector3i getTileMapResolution(int level) const;
    /// Returns the number of voxels at the given level.
    exint getNumVoxels(int level) const;
    /// Returns the number of active tiles at the given level.
    int getNumTiles(int level) const;
    /// Returns the number of tile map levels (depth of the operator).
    int getNumLevels() const;

    /// Build the tile Laplacian matrix at the specified level. Only works for non-
    /// negative levels. Returns the number of rows/columns.
    int buildMatrix(UT_MatrixT<T>& M, int level) const;

protected:
    /// Helper function to create the finest level tile map from the given reference
    /// voxel array. If tileMap is false, tiles of the array are considered active
    /// if they are not constant and equal to the inactive argument. If tileMap is
    /// true, ref is assumed to be the same resolution as the finest tile map, and
    /// a tile is active if its value in ref is non-zero.
    /// Spacing and boundary conditions must be set prior to calling this function.
    template<typename S>
    void initTileMap0(const UT_VoxelArray<S>& ref, S inactive, bool tileMap);
    /// Adds a tile map one level coarser than the coarsest in the atlas.
    void addTileMapLevel();

    /// Builds the coefficient matrix for the linear system at the last level. Also
    /// factors it to allow for fast linear solves.
    void factorFinalLevel();
    /// Helper function to for initialization. Coarsens the map until the degrees of
    /// freedom are sufficiently reduced, factors the final level, and allocates the
    /// scratchpad.
    void finishInit();

    /// Returns the off-diagonal value in the matrix between the node at the given
    /// location and the next one along the supplied axis; works at the provided
    /// level. Note that axis must be in [0, 5], coresponding to -X, +X, -Y, +Y, -Z,
    /// +Z.
    inline T getEdgeWeight(int x, int y, int z, int level, int axis) const;
    /// Returns the diagonal value in the matrix at the given voxel and level.
    inline T getDiagonalWeight(int x, int y, int z, int level) const;

    /// Helper function for mgRestrict() and mgProlong(), for negative levels.
    /// If RESTRICT is true, src is restricted to a coarser level and stored in dst.
    /// If RESTRICT is false, src is prolongated to a finer level and stored in dst.
    /// LEVEL must be pre-negated (so it is positive).
    template<bool RESTRICT, int LEVEL>
    inline void mgRestrictProlongN(const UT_VectorT<T>& src, UT_VectorT<T>& dst)
        const;
    /// Helper function for mgRestrict(), for non-negative levels.
    /// If RESTRICT is true, src is restricted to a coarser level and stored in dst.
    /// If RESTRICT is false, src is prolongated to a finer level and stored in dst.
    template<bool RESTRICT>
    inline void mgRestrictProlongP(const UT_VectorT<T>& src, UT_VectorT<T>& dst,
                                   int level) const;

    /// Helper function for applying the Laplacian operator. LEVEL must be pre-negated
    /// (so it is positive). If SUBTRACT is true computes b-Lx, otherwise calculates
    /// Lx.
    template<int LEVEL, bool SUBTRACT>
    inline void applyLaplacianInternalN(const UT_VectorT<T>& x, const UT_VectorT<T>& b,
                                        UT_VectorT<T>& d) const;
    /// Helper function for subtractApplyLaplacian(), for non-negative levels. If
    /// SUBTRACT is true, computes b-Lx, otherwise calculates Lx.
    template<bool SUBTRACT>
    inline void applyLaplacianInternalP(const UT_VectorT<T>& x, const UT_VectorT<T>& b,
                                        UT_VectorT<T>& d, int level) const;

    /// Performs a single step of red-black Gauss-Seidel relaxation for voxels of the
    /// given parity. Note that this helper function is intended for negative levels
    /// (multiple voxels per tile). The template parameter must be pre-negated (so it
    /// is positive).
    /// If ZERO_GUESS is true, FIRST_PASS signals whether or not one of the colors has
    /// already been processed. The difference is that on the second pass, values of
    /// the previous color must be used, but colors of the current color are still
    /// assumed to be zero.
    /// This generic function implements the relaxation for all types and negative
    /// levels.
    template<int LEVEL, bool ZERO_GUESS, bool FIRST_PASS>
    inline void smoothLaplacianInternalGS(UT_VectorT<T>& x, const UT_VectorT<T>& b,
                                          T omega, int parity) const;
    /// Performs a single step of red-black Gauss-Seidel relaxation for voxels of the
    /// given parity at level -1. This is a specialized method that either calls the
    /// generic function or does the work itself using SIMD.
    template<bool ZERO_GUESS, bool FIRST_PASS>
    inline void smoothLaplacianInternalGS_V1(UT_VectorT<T>& x, const UT_VectorT<T>& b,
                                             T omega, int parity) const;
    /// Performs a single step of red-black Gauss-Seidel relaxation for voxels of the
    /// given parity at level -2. This is a specialized method that either calls the
    /// generic function or does the work itself using SIMD.
    template<bool ZERO_GUESS, bool FIRST_PASS>
    inline void smoothLaplacianInternalGS_V2(UT_VectorT<T>& x, const UT_VectorT<T>& b,
                                             T omega, int parity) const;
    /// Performs a single step of red-black Gauss-Seidel relaxation for voxels of the
    /// given parity at level -3 or -4. This is a specialized method that either calls
    /// the generic function or does the work itself using SIMD.
    template<int LEVEL, bool ZERO_GUESS, bool FIRST_PASS>
    inline void smoothLaplacianInternalGS_V34(UT_VectorT<T>& x, const UT_VectorT<T>& b,
                                              T omega, int parity) const;
    /// Performs a single step of Jacobi relaxation. This helper function is intended
    /// for negative levels (the template parameter must be pre-negated so it is
    /// positive.
    /// CAUTION: x and x0 must be distinct!
    /// b vector is unused if relaxation is homogeneous. Likewise, omega is unused (no
    /// damping).
    /// This generic function implements the relaxation for all types and negative
    /// levels.
    template<int LEVEL, bool HOMOGENEOUS>
    inline void smoothLaplacianInternalJN(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
                                          const UT_VectorT<T>& b, T omega) const;
    /// Performs a single step of Jacobi relaxation at level -1. This is a specialized
    /// method that either calls the generic function or does the work itself using SIMD.
    template<bool HOMOGENEOUS>
    inline void smoothLaplacianInternalJN_V1(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
                                             const UT_VectorT<T>& b, T omega) const;
    /// Performs a single step of Jacobi relaxation at levels -2, -3, or -4. This is a
    /// specialized method that either calls the generic function or does the work itself
    /// using SIMD.
    template<int LEVEL, bool HOMOGENEOUS>
    inline void smoothLaplacianInternalJN_V234(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
                                               const UT_VectorT<T>& b, T omega) const;
    /// Performs a single step of Jacobi relaxation. This helper function is intended
    /// for non-negative levels.
    /// CAUTION: x and x0 must be distinct!
    /// b vector is unused if relaxation is homogeneous. Likewise, omega is unused (no
    /// damping).
    template<bool HOMOGENEOUS, bool ZERO_GUESS>
    inline void smoothLaplacianInternalJP(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
                                          const UT_VectorT<T>& b, T omega,
                                          int level) const;

    /// Puts pointers to the blocks from the given vector into the data array.
    /// Border compensations for adjusting the diagonal value are also placed in the
    /// border array.
    template<int LEVEL>
    void getTilePointers(int tile, const UT_VectorT<T>& b, T* data[7], T border[6])
        const;

    /// Multithreaded helper for copyToVoxels. dst should be set to the correct size
    /// before calling this.
    THREADED_METHOD3_CONST(UT_MGOperatorT, dst.numTiles() > 500,
                           copyToVoxelsHelper,
                           const UT_VectorT<T>&, src,
                           UT_VoxelArray<T>&, dst,
                           bool, topad);
    void copyToVoxelsHelperPartial(const UT_VectorT<T>& src, UT_VoxelArray<T>& dst,
                                   bool toPad, const UT_JobInfo& info) const;
                                   

public:
    /// This struct contains the coordinates of the current tile, as well as its
    /// Laplacian coefficients.
    /// Coefficient order convention: -X, +X, -Y, +Y, -Z, +Z, diagonal.
    struct utMGCoordAndCoeff
    {
        UT_Vector3i                     coord;
        UT_FixedVector<T, 7>            coeff;
    };

    /// This struct holds the tile map (which gives the linear tile number for each
    /// tile location), its inverse map (which produces location of a tile from its
    /// linear index), and the coefficients for the row corresponding to the tile.
    struct utMGTileMap
    {
        UT_VoxelArray<int>              tileMap;
        UT_Array<utMGCoordAndCoeff>     invTileMapMat;
    };

protected:
    /// Internal coefficients for the Laplacian operator; first three hold the
    /// weights along each axis, followed by the diagonal weight.
    UT_Vector4T<T>                      myInternalCoeff;
    /// Conditions at the lower/left/back boundaries of the map.
    UT_Vector3I                         myBoundariesNeg;
    /// Conditions at the upper/right/front boundaries of the map.
    UT_Vector3I                         myBoundariesPos;

    /// Array of tile map levels.
    UT_Array<utMGTileMap>               myAtlas;

    /// The factored coefficient matrix for the linear system at the coarsest tile
    /// level.
    UT_MatrixT<T>                       myCoarseCholeskyL;
    UT_VectorT<T>                       myCoarseCholeskyD;

    /// The tile map is coarsened until its number of voxels is smaller than this
    /// value.
    const int                           myVoxelCountThreshold;

    /// Array of pre-allocated temporary vectors for use in the Poisson solve (2 per
    /// level).
    UT_Array<UT_VectorT<T>>             myScratchpad;

public:
    /// Number of pre- and post-smoothing passes to perform in a V-cycle.
    int                                 mySmoothingPasses;
    /// Number of homogeneous smoothing passes to perform after prolongation in a
    /// V- and full multigrid cycles.
    int                                 myDiffusionPasses;
    /// Amplification factor after prolongation (for V- and full multigrid cycles).
    T                                   myAmplificationFactor;
    /// Damping parameter for Jacobi smoothing passes.
    T                                   myJacobiOmega;
    /// Overrelaxation parameter for Gauss-Seidel smoothing passes.
    T                                   myGSOmega;
};

typedef UT_MGOperatorT<fpreal>          UT_MGOperatorR;
typedef UT_MGOperatorT<fpreal32>        UT_MGOperatorF;
typedef UT_MGOperatorT<fpreal64>        UT_MGOperatorD;

UT_EXTERN_TEMPLATE(UT_MGOperatorT<fpreal32>);
UT_EXTERN_TEMPLATE(UT_MGOperatorT<fpreal64>);

#endif