HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_MGOperator.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_MGOperator.h ( UT Library, C++)
7  *
8  * COMMENTS:
9  * This class encapsulates a tiled discrete Laplacian operator that
10  * is algebraically coarsened, allowing for a sparse Poisson solve.
11  */
12 
13 #ifndef __UT_MGOperator__
14 #define __UT_MGOperator__
15 
16 #include "UT_API.h"
17 #include "UT_MatrixSolver.h"
18 #include "UT_VoxelArray.h"
19 
20 /// This function returns true if the given tile has at least one voxel whose
21 /// value exceeds 0.5.
22 template<typename S> UT_API
23 bool ut_tileHasActiveVoxel(const UT_VoxelTile<S>* tile);
24 
25 template<typename T>
27 {
28 public:
29  /// Once initialized, this operator is algebraically coarsened until the number
30  /// of degrees of freedom does not exceed the provided threshold.
31  UT_MGOperatorT(int threshold = 100) : myVoxelCountThreshold(threshold),
32  mySmoothingPasses(2),
33  myDiffusionPasses(1),
34  myAmplificationFactor(1.8),
35  myJacobiOmega(2.0 / 3),
36  myGSOmega(4.0 / 3) {}
38 
39  /// Initializes this operator to the Laplacian matrix with the given spacing and
40  /// boundary conditions. This method also takes care of setting up the coarser
41  /// operators. Tiles in the operator corresponds to tiles of the reference array.
42  /// Active tiles are non-constant or different from inactive.
43  /// All dimensions of the reference array must be divisible by 16.
44  template<typename S>
45  void init(const UT_VoxelArray<S>& ref, S inactive,
46  const UT_Vector3T<T>& spacing,
47  const UT_Vector3I& boundariesNeg,
48  const UT_Vector3I& boundariesPos);
49  /// This version of initialization takes a validity map (whose resolution should
50  /// match that of the finest tile map) and performs the necessary setup. Active
51  /// tiles have a non-zero value in the reference voxel array.
52  /// Resolution of this operator will be 16 times the resolution of the reference
53  /// array.
54  template<typename S>
55  void init(const UT_VoxelArray<S>& ref,
56  const UT_Vector3T<T>& spacing,
57  const UT_Vector3I& boundariesNeg,
58  const UT_Vector3I& boundariesPos);
59  /// This version of initialization accepts a stencil array. Tiles in the operator
60  /// are activated if the corresponding tile in the stencil array has at least one
61  /// voxel whose value exceeds 0.5.
62  /// All dimension of the reference array must be divisible by 16.
63  template<typename S>
64  void initFromStencil(const UT_VoxelArray<S>& stencil,
65  const UT_Vector3T<T>& spacing,
66  const UT_Vector3I& boundariesNeg,
67  const UT_Vector3I& boundariesPos);
68 
69  /// Applies the restriction operator from level l to level l+1.
70  void mgRestrict(const UT_VectorT<T>& src, UT_VectorT<T>& dst, int l) const;
71  /// Applies the prolongation operator from level l+1 to level l.
72  void mgProlong(const UT_VectorT<T>& src, UT_VectorT<T>& dst, int l) const;
73 
74  /// Flattens contents of the given voxel array into the vector.
75  void copyFromVoxels(const UT_VoxelArray<T>& src, UT_VectorT<T>& dst) const;
76  /// Copies the contents of the flattened vector into the voxel array. If topad
77  /// is true, then destination array will have an extra layer of voxels
78  /// compared to the internal tile map, whose entries will be set to 0.
79  void copyToVoxels(const UT_VectorT<T>& src, UT_VoxelArray<T>& dst,
80  bool topad) const;
81 
82  /// Solves the linear system using multigrid.
83  /// The solution is stored in x, which will be initialized to the appropriate
84  /// size. If b has incomplete tiles, x will be padded with an extra layer of
85  /// voxels along every dimension.
86  int solvePoisson(fpreal64 abstol, fpreal64 reltol, int maxiter,
88  UT_ValArray<fpreal64>* resnorminf,
89  UT_ValArray<fpreal64>* resnorm2);
90 
91  /// Calculates b-Lx at the given level and puts the answer in d (note that b
92  /// and d can be the same vector).
93  void subtractApplyLaplacian(const UT_VectorT<T>& x, const UT_VectorT<T>& b,
94  UT_VectorT<T>& d, int level) const;
95  /// Calculates Lx at the given level and puts the answer in b.
96  void applyLaplacian(const UT_VectorT<T>& x, UT_VectorT<T>& b, int level) const;
97 
98  /// Does a full multigrid solve to get an approximation for the solution to
99  /// the equation Lx=b at the given level. Note that x must be pre-initialized,
100  /// but its incoming contents are irrelevant.
101  /// t0 and t1 are used to perform intermediate work and must be pre-allocated.
102  void fullMultigrid(UT_VectorT<T>& x, const UT_VectorT<T>& b, UT_VectorT<T>& t0,
103  UT_VectorT<T>& t1, int level);
104  /// Approximates a solution for the equation Lx=b using a single V-cycle. Note
105  /// that x will be re-initialized, so its incoming contents are irrelevant.
106  /// t is a temporary vector that must be allocated to the same size as x.
107  void vcycle(UT_VectorT<T>& x, UT_VectorT<T>& t, const UT_VectorT<T>& b,
108  int level);
109  /// Does a single iteration of relaxation on the solution vector at the given
110  /// level. For positive levels, uses Jacobi relaxation; for negative levels,
111  /// performs red-black Gauss-Seidel, with redFirst indicating which color to
112  /// smooth first.
113  /// For Jacobi smoothing (non-negative levels) x0 is the vector to relax, x
114  /// will hold the result. These vectors must be distinct! For Gauss-Seidel
115  /// smoothing (negative levels), x holds the input and output solutions, so x0
116  /// is unused.
117  /// If ZERO_GUESS is true, it is assumed that the incoming vector is zero.
118  template<bool ZERO_GUESS = false>
119  void smoothLaplacian(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
120  const UT_VectorT<T>& b, T omega, int level,
121  bool redFirst) const;
122  /// Does a single iteration of homogeneous relaxation (smoothing for the system
123  /// Lx=0). Uses Jacobi relaxation without damping.
124  /// x0 is the vector to relax, x will hold the results; these vectors must be
125  /// distinct!
126  void smoothLaplacianHomogeneous(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
127  int level) const;
128  /// Performs a direct solve at the coarsest level.
129  void directSolve(UT_VectorT<T>& x, const UT_VectorT<T>& b) const;
130 
131  /// Returns the tile map's resolution at the given level.
132  UT_Vector3i getTileMapResolution(int level) const;
133  /// Returns the number of voxels at the given level.
134  exint getNumVoxels(int level) const;
135  /// Returns the number of active tiles at the given level.
136  int getNumTiles(int level) const;
137  /// Returns the number of tile map levels (depth of the operator).
138  int getNumLevels() const;
139 
140  /// Build the tile Laplacian matrix at the specified level. Only works for non-
141  /// negative levels. Returns the number of rows/columns.
142  int buildMatrix(UT_MatrixT<T>& M, int level) const;
143 
144 protected:
145  /// Helper function to create the finest level tile map from the given reference
146  /// voxel array. If tileMap is false, tiles of the array are considered active
147  /// if they are not constant and equal to the inactive argument. If tileMap is
148  /// true, ref is assumed to be the same resolution as the finest tile map, and
149  /// a tile is active if its value in ref is non-zero.
150  /// Spacing and boundary conditions must be set prior to calling this function.
151  template<typename S>
152  void initTileMap0(const UT_VoxelArray<S>& ref, S inactive, bool tileMap);
153  /// Adds a tile map one level coarser than the coarsest in the atlas.
154  void addTileMapLevel();
155 
156  /// Builds the coefficient matrix for the linear system at the last level. Also
157  /// factors it to allow for fast linear solves.
158  void factorFinalLevel();
159  /// Helper function to for initialization. Coarsens the map until the degrees of
160  /// freedom are sufficiently reduced, factors the final level, and allocates the
161  /// scratchpad.
162  void finishInit();
163 
164  /// Returns the off-diagonal value in the matrix between the node at the given
165  /// location and the next one along the supplied axis; works at the provided
166  /// level. Note that axis must be in [0, 5], coresponding to -X, +X, -Y, +Y, -Z,
167  /// +Z.
168  inline T getEdgeWeight(int x, int y, int z, int level, int axis) const;
169  /// Returns the diagonal value in the matrix at the given voxel and level.
170  inline T getDiagonalWeight(int x, int y, int z, int level) const;
171 
172  /// Helper function for mgRestrict() and mgProlong(), for negative levels.
173  /// If RESTRICT is true, src is restricted to a coarser level and stored in dst.
174  /// If RESTRICT is false, src is prolongated to a finer level and stored in dst.
175  /// LEVEL must be pre-negated (so it is positive).
176  template<bool RESTRICT, int LEVEL>
177  inline void mgRestrictProlongN(const UT_VectorT<T>& src, UT_VectorT<T>& dst)
178  const;
179  /// Helper function for mgRestrict(), for non-negative levels.
180  /// If RESTRICT is true, src is restricted to a coarser level and stored in dst.
181  /// If RESTRICT is false, src is prolongated to a finer level and stored in dst.
182  template<bool RESTRICT>
183  inline void mgRestrictProlongP(const UT_VectorT<T>& src, UT_VectorT<T>& dst,
184  int level) const;
185 
186  /// Helper function for applying the Laplacian operator. LEVEL must be pre-negated
187  /// (so it is positive). If SUBTRACT is true computes b-Lx, otherwise calculates
188  /// Lx.
189  template<int LEVEL, bool SUBTRACT>
190  inline void applyLaplacianInternalN(const UT_VectorT<T>& x, const UT_VectorT<T>& b,
191  UT_VectorT<T>& d) const;
192  /// Helper function for subtractApplyLaplacian(), for non-negative levels. If
193  /// SUBTRACT is true, computes b-Lx, otherwise calculates Lx.
194  template<bool SUBTRACT>
195  inline void applyLaplacianInternalP(const UT_VectorT<T>& x, const UT_VectorT<T>& b,
196  UT_VectorT<T>& d, int level) const;
197 
198  /// Performs a single step of red-black Gauss-Seidel relaxation for voxels of the
199  /// given parity. Note that this helper function is intended for negative levels
200  /// (multiple voxels per tile). The template parameter must be pre-negated (so it
201  /// is positive).
202  /// If ZERO_GUESS is true, FIRST_PASS signals whether or not one of the colors has
203  /// already been processed. The difference is that on the second pass, values of
204  /// the previous color must be used, but colors of the current color are still
205  /// assumed to be zero.
206  template<int LEVEL, bool ZERO_GUESS, bool FIRST_PASS>
207  inline void smoothLaplacianInternalGS(UT_VectorT<T>& x, const UT_VectorT<T>& b,
208  T omega, int parity) const;
209  /// Performs a single step of Jacobi relaxation. This helper function is intended
210  /// for negative levels (the template parameter must be pre-negated so it is
211  /// positive.
212  /// CAUTION: x and x0 must be distinct!
213  /// b vector is unused if relaxation is homogeneous. Likewise, omega is unused (no
214  /// damping).
215  template<int LEVEL, bool HOMOGENEOUS>
216  inline void smoothLaplacianInternalJN(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
217  const UT_VectorT<T>& b, T omega) const;
218  /// Performs a single step of Jacobi relaxation. This helper function is intended
219  /// for non-negative levels.
220  /// CAUTION: x and x0 must be distinct!
221  /// b vector is unused if relaxation is homogeneous. Likewise, omega is unused (no
222  /// damping).
223  template<bool HOMOGENEOUS, bool ZERO_GUESS>
224  inline void smoothLaplacianInternalJP(UT_VectorT<T>& x, const UT_VectorT<T>& x0,
225  const UT_VectorT<T>& b, T omega,
226  int level) const;
227 
228  /// Puts pointers to the blocks from the given vector into the data array.
229  /// Border compensations for adjusting the diagonal value are also placed in the
230  /// border array.
231  template<int LEVEL>
232  void getTilePointers(int tile, const UT_VectorT<T>& b, T* data[7], T border[6])
233  const;
234 
235 public:
236  /// This struct contains the coordinates of the current tile, as well as its
237  /// Laplacian coefficients.
238  /// Coefficient order convention: -X, +X, -Y, +Y, -Z, +Z, diagonal.
240  {
243  };
244 
245  /// This struct holds the tile map (which gives the linear tile number for each
246  /// tile location), its inverse map (which produces location of a tile from its
247  /// linear index), and the coefficients for the row corresponding to the tile.
248  struct utMGTileMap
249  {
252  };
253 
254 protected:
255  /// Internal coefficients for the Laplacian operator; first three hold the
256  /// weights along each axis, followed by the diagonal weight.
258  /// Conditions at the lower/left/back boundaries of the map.
260  /// Conditions at the upper/right/front boundaries of the map.
262 
263  /// Array of tile map levels.
265 
266  /// The factored coefficient matrix for the linear system at the coarsest tile
267  /// level.
270 
271  /// The tile map is coarsened until its number of voxels is smaller than this
272  /// value.
274 
275  /// Array of pre-allocated temporary vectors for use in the Poisson solve (2 per
276  /// level).
278 
279 public:
280  /// Number of pre- and post-smoothing passes to perform in a V-cycle.
282  /// Number of homogeneous smoothing passes to perform after prolongation in a
283  /// V- and full multigrid cycles.
285  /// Amplification factor after prolongation (for V- and full multigrid cycles).
287  /// Damping parameter for Jacobi smoothing passes.
289  /// Overrelaxation parameter for Gauss-Seidel smoothing passes.
291 };
292 
296 
299 
300 #endif
301 
UT_MGOperatorT(int threshold=100)
Definition: UT_MGOperator.h:31
GLenum src
Definition: glew.h:2410
GLenum GLint ref
Definition: glew.h:1845
UT_Array< utMGTileMap > myAtlas
Array of tile map levels.
#define UT_API
Definition: UT_API.h:13
UT_Vector3I myBoundariesPos
Conditions at the upper/right/front boundaries of the map.
GLdouble l
Definition: glew.h:9122
GLint GLfloat GLint stencil
Definition: glew.h:2167
3D Vector class.
4D Vector class.
Definition: UT_Vector4.h:163
UT_API bool ut_tileHasActiveVoxel(const UT_VoxelTile< S > *tile)
GLdouble GLdouble z
Definition: glew.h:1559
int mySmoothingPasses
Number of pre- and post-smoothing passes to perform in a V-cycle.
UT_MGOperatorT< fpreal32 > UT_MGOperatorF
GLint GLint GLint GLint GLint x
Definition: glew.h:1252
GLint GLint GLint GLint GLint GLint y
Definition: glew.h:1252
int64 exint
Definition: SYS_Types.h:120
GLint GLenum GLsizei GLint GLsizei const void * data
Definition: glew.h:1379
double fpreal64
Definition: SYS_Types.h:196
GLint GLint GLsizei GLsizei GLsizei GLint border
Definition: glew.h:1254
UT_VectorT< T > myCoarseCholeskyD
GLenum GLenum dst
Definition: glew.h:2410
UT_VoxelArray< int > tileMap
const int myVoxelCountThreshold
T myAmplificationFactor
Amplification factor after prolongation (for V- and full multigrid cycles).
UT_Vector4T< T > myInternalCoeff
UT_MatrixT< T > myCoarseCholeskyL
UT_Vector3I myBoundariesNeg
Conditions at the lower/left/back boundaries of the map.
GLuint GLfloat x0
Definition: glew.h:12681
GLdouble GLdouble GLdouble b
Definition: glew.h:9122
UT_Array< utMGCoordAndCoeff > invTileMapMat
UT_EXTERN_TEMPLATE(UT_MGOperatorT< fpreal32 >)
UT_Array< UT_VectorT< T > > myScratchpad
UT_MGOperatorT< fpreal64 > UT_MGOperatorD
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
Definition: glew.h:12681
T myGSOmega
Overrelaxation parameter for Gauss-Seidel smoothing passes.
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
Definition: glew.h:12681
UT_FixedVector< T, 7, true > coeff
GLint level
Definition: glew.h:1252
GLdouble GLdouble t
Definition: glew.h:1398
UT_MGOperatorT< fpreal > UT_MGOperatorR
T myJacobiOmega
Damping parameter for Jacobi smoothing passes.