1 /*
2 Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
3 All rights reserved.
4
5 Redistribution and use in source and binary forms, with or without modification,
6 are permitted provided that the following conditions are met:
7
8 Redistributions of source code must retain the above copyright notice, this list of
9 conditions and the following disclaimer. Redistributions in binary form must reproduce
10 the above copyright notice, this list of conditions and the following disclaimer
11 in the documentation and/or other materials provided with the distribution.
12
13 Neither the name of the Johns Hopkins University nor the names of its contributors
14 may be used to endorse or promote products derived from this software without specific
15 prior written permission.
16
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
18 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
19 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
20 SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
21 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
22 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
23 BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
25 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
26 DAMAGE.
27 */
28 // [COMMENTS]
29 // -- Throughout the code, should make a distinction between indices and offsets
30 // -- Make an instance of _evaluate that samples the finite-elements correctly (specifically, to handle the boundaries)
31 // -- Make functions like depthAndOffset parity dependent (ideally all "depth"s should be relative to the B-Slpline resolution
32 // -- Make all points relative to the unit-cube, regardless of degree parity
33 // -- It's possible that for odd degrees, the iso-surfacing will fail because the leaves in the SortedTreeNodes do not form a partition of space
34 // -- [MAYBE] Treat normal field as a sum of delta functions, rather than a smoothed signal (again, so that high degrees aren't forced to generate smooth reconstructions)
35 // -- [MAYBE] Make the degree of the B-Spline with which the normals are splatted independent of the degree of the FEM system. (This way, higher degree systems aren't forced to generate smoother normal fields.)
36 // -- [MAYBE] Remove the isValidFEM/isValidSpace functions since the octree supports all degrees/boundary types (up to the max degree for which finalizedBrooded... was called)
37
38 // [TODO]
39 // -- Currently, the implementation assumes that the boundary constraints are the same for vector fields and scalar fields
40 // -- Modify the setting of the flags so that only the subset of the broods that are needed
41
42 #ifndef MULTI_GRID_OCTREE_DATA_INCLUDED
43 #define MULTI_GRID_OCTREE_DATA_INCLUDED
44
45 #define NEW_CODE
46 #define FAST_SET_UP // If enabled, kernel density estimation is done aglomeratively
47
48 #define POINT_DATA_RES 0 // Specifies the resolution of the subgrid storing points with each voxel (0==1 but is faster)
49
50 #define DATA_DEGREE 1 // The order of the B-Spline used to splat in data for color interpolation
51 #define WEIGHT_DEGREE 2 // The order of the B-Spline used to splat in the weights for density estimation
52 #define NORMAL_DEGREE 2 // The order of the B-Spline used to splat int the normals for constructing the Laplacian constraints
53 //#define MAX_MEMORY_GB 15 // The maximum memory the application is allowed to use
54 #define MAX_MEMORY_GB 0
55
56 #include <unordered_map>
57 #ifdef _OPENMP
58 #include <omp.h>
59 #endif // _OPENMP
60 #include "BSplineData.h"
61 #include "PointStream.h"
62 #include "Geometry.h"
63 #include "Octree.h"
64 #include "SparseMatrix.h"
65
66 #ifndef _OPENMP
omp_get_num_procs(void)67 int omp_get_num_procs( void ){ return 1; }
omp_get_thread_num(void)68 int omp_get_thread_num( void ){ return 0; }
69 #endif // _OPENMP
70
71 #define DERIVATIVES( Degree ) ( ( Degree>1 ) ? 2 : ( Degree==1 ? 1 : 0 ) )
72
73 class TreeNodeData
74 {
75 public:
76 enum
77 {
78 SPACE_FLAG = 1 ,
79 FEM_FLAG = 2 ,
80 GHOST_FLAG = 1<<7
81 };
82 int nodeIndex;
83 char flags;
84
setGhostFlag(bool f)85 void setGhostFlag( bool f ){ if( f ) flags |= GHOST_FLAG ; else flags &= ~GHOST_FLAG; }
getGhostFlag(void)86 bool getGhostFlag( void ) const { return ( flags & GHOST_FLAG )!=0; }
87 TreeNodeData( void );
88 ~TreeNodeData( void );
89 };
90
91 class VertexData
92 {
93 typedef OctNode< TreeNodeData > TreeOctNode;
94 public:
95 static const int VERTEX_COORDINATE_SHIFT = ( sizeof( long long ) * 8 ) / 3;
96 static long long EdgeIndex( const TreeOctNode* node , int eIndex , int maxDepth , int index[DIMENSION] );
97 static long long EdgeIndex( const TreeOctNode* node , int eIndex , int maxDepth );
98 static long long FaceIndex( const TreeOctNode* node , int fIndex , int maxDepth,int index[DIMENSION] );
99 static long long FaceIndex( const TreeOctNode* node , int fIndex , int maxDepth );
100 static long long CornerIndex( const TreeOctNode* node , int cIndex , int maxDepth , int index[DIMENSION] );
101 static long long CornerIndex( const TreeOctNode* node , int cIndex , int maxDepth );
102 static long long CenterIndex( const TreeOctNode* node , int maxDepth , int index[DIMENSION] );
103 static long long CenterIndex( const TreeOctNode* node , int maxDepth );
104 static long long CornerIndex( int depth , const int offSet[DIMENSION] , int cIndex , int maxDepth , int index[DIMENSION] );
105 static long long CenterIndex( int depth , const int offSet[DIMENSION] , int maxDepth , int index[DIMENSION] );
106 static long long CornerIndexKey( const int index[DIMENSION] );
107 };
108
109 // This class stores the octree nodes, sorted by depth and then by z-slice.
110 // To support primal representations, the initializer takes a function that
111 // determines if a node should be included/indexed in the sorted list.
112 // [NOTE] Indexing of nodes is _GLOBAL_
113 class SortedTreeNodes
114 {
115 typedef OctNode< TreeNodeData > TreeOctNode;
116 protected:
117 Pointer( Pointer( int ) ) _sliceStart;
118 int _levels;
119 public:
120 Pointer( TreeOctNode* ) treeNodes;
begin(int depth)121 int begin( int depth ) const{ return _sliceStart[depth][0]; }
end(int depth)122 int end( int depth ) const{ return _sliceStart[depth][(size_t)1<<depth]; }
begin(int depth,int slice)123 int begin( int depth , int slice ) const{ return _sliceStart[depth][slice ] ; }
end(int depth,int slice)124 int end( int depth , int slice ) const{ if(depth<0||depth>=_levels||slice<0||slice>=(1<<depth)) printf( "uh oh\n" ) ; return _sliceStart[depth][slice+1]; }
size(void)125 int size( void ) const { return _sliceStart[_levels-1][(size_t)1<<(_levels-1)]; }
size(int depth)126 int size( int depth ) const { if(depth<0||depth>=_levels) printf( "uhoh\n" ); return _sliceStart[depth][(size_t)1<<depth] - _sliceStart[depth][0]; }
size(int depth,int slice)127 int size( int depth , int slice ) const { return _sliceStart[depth][slice+1] - _sliceStart[depth][slice]; }
levels(void)128 int levels( void ) const { return _levels; }
129
130 SortedTreeNodes( void );
131 ~SortedTreeNodes( void );
132 void set( TreeOctNode& root , std::vector< int >* map );
133 void set( TreeOctNode& root );
134
135 template< int Indices >
136 struct _Indices
137 {
138 int idx[Indices];
_Indices_Indices139 _Indices( void ){ memset( idx , -1 , sizeof( int ) * Indices ); }
140 int& operator[] ( int i ) { return idx[i]; }
141 const int& operator[] ( int i ) const { return idx[i]; }
142 };
143 typedef _Indices< Square::CORNERS > SquareCornerIndices;
144 typedef _Indices< Square::EDGES > SquareEdgeIndices;
145 typedef _Indices< Square::FACES > SquareFaceIndices;
146
147 struct SliceTableData
148 {
149 Pointer( SquareCornerIndices ) cTable;
150 Pointer( SquareEdgeIndices ) eTable;
151 Pointer( SquareFaceIndices ) fTable;
152 int cCount , eCount , fCount , nodeOffset , nodeCount;
SliceTableDataSliceTableData153 SliceTableData( void ){ fCount = eCount = cCount = 0 , cTable = NullPointer( SquareCornerIndices ) , eTable = NullPointer( SquareEdgeIndices ) , fTable = NullPointer( SquareFaceIndices ) , _cMap = _eMap = _fMap = NullPointer( int ); }
~SliceTableDataSliceTableData154 ~SliceTableData( void ){ clear(); }
155 #ifdef BRUNO_LEVY_FIX
clearSliceTableData156 void clear( void ){ DeletePointer( cTable ) ; DeletePointer( eTable ) ; DeletePointer( fTable ) ; DeletePointer( _cMap ) ; DeletePointer( _eMap ) ; DeletePointer( _fMap ) ; fCount = eCount = cCount = 0; }
157 #else // !BRUNO_LEVY_FIX
clearSliceTableData158 void clear( void ){ DeletePointer( cTable ) ; DeletePointer( eTable ) ; DeletePointer( fTable ) ; fCount = eCount = cCount = 0; }
159 #endif // BRUNO_LEVY_FIX
160 SquareCornerIndices& cornerIndices( const TreeOctNode* node );
161 SquareCornerIndices& cornerIndices( int idx );
162 const SquareCornerIndices& cornerIndices( const TreeOctNode* node ) const;
163 const SquareCornerIndices& cornerIndices( int idx ) const;
164 SquareEdgeIndices& edgeIndices( const TreeOctNode* node );
165 SquareEdgeIndices& edgeIndices( int idx );
166 const SquareEdgeIndices& edgeIndices( const TreeOctNode* node ) const;
167 const SquareEdgeIndices& edgeIndices( int idx ) const;
168 SquareFaceIndices& faceIndices( const TreeOctNode* node );
169 SquareFaceIndices& faceIndices( int idx );
170 const SquareFaceIndices& faceIndices( const TreeOctNode* node ) const;
171 const SquareFaceIndices& faceIndices( int idx ) const;
172 protected:
173 Pointer( int ) _cMap;
174 Pointer( int ) _eMap;
175 Pointer( int ) _fMap;
176 friend class SortedTreeNodes;
177 };
178 struct XSliceTableData
179 {
180 Pointer( SquareCornerIndices ) eTable;
181 Pointer( SquareEdgeIndices ) fTable;
182 int fCount , eCount , nodeOffset , nodeCount;
XSliceTableDataXSliceTableData183 XSliceTableData( void ){ fCount = eCount = 0 , eTable = NullPointer( SquareCornerIndices ) , fTable = NullPointer( SquareEdgeIndices ) , _eMap = _fMap = NullPointer( int ); }
~XSliceTableDataXSliceTableData184 ~XSliceTableData( void ){ clear(); }
185 #ifdef BRUNO_LEVY_FIX
clearXSliceTableData186 void clear( void ) { DeletePointer( fTable ) ; DeletePointer( eTable ) ; DeletePointer( _eMap ) ; DeletePointer( _fMap ) ; fCount = eCount = 0; }
187 #else // !BRUNO_LEVY_FIX
clearXSliceTableData188 void clear( void ) { DeletePointer( fTable ) ; DeletePointer( eTable ) ; fCount = eCount = 0; }
189 #endif // BRUNO_LEVY_FIX
190 SquareCornerIndices& edgeIndices( const TreeOctNode* node );
191 SquareCornerIndices& edgeIndices( int idx );
192 const SquareCornerIndices& edgeIndices( const TreeOctNode* node ) const;
193 const SquareCornerIndices& edgeIndices( int idx ) const;
194 SquareEdgeIndices& faceIndices( const TreeOctNode* node );
195 SquareEdgeIndices& faceIndices( int idx );
196 const SquareEdgeIndices& faceIndices( const TreeOctNode* node ) const;
197 const SquareEdgeIndices& faceIndices( int idx ) const;
198 protected:
199 Pointer( int ) _eMap;
200 Pointer( int ) _fMap;
201 friend class SortedTreeNodes;
202 };
203 void setSliceTableData ( SliceTableData& sData , int depth , int offset , int threads ) const;
204 void setXSliceTableData( XSliceTableData& sData , int depth , int offset , int threads ) const;
205 };
206
207 template< int Degree >
208 struct PointSupportKey : public OctNode< TreeNodeData >::NeighborKey< BSplineSupportSizes< Degree >::SupportEnd , -BSplineSupportSizes< Degree >::SupportStart >
209 {
210 static const int LeftRadius = BSplineSupportSizes< Degree >::SupportEnd;
211 static const int RightRadius = -BSplineSupportSizes< Degree >::SupportStart;
212 static const int Size = LeftRadius + RightRadius + 1;
213 };
214 template< int Degree >
215 struct ConstPointSupportKey : public OctNode< TreeNodeData >::ConstNeighborKey< BSplineSupportSizes< Degree >::SupportEnd , -BSplineSupportSizes< Degree >::SupportStart >
216 {
217 static const int LeftRadius = BSplineSupportSizes< Degree >::SupportEnd;
218 static const int RightRadius = -BSplineSupportSizes< Degree >::SupportStart;
219 static const int Size = LeftRadius + RightRadius + 1;
220 };
221
222 template< class Real , bool HasGradients >
223 struct SinglePointData
224 {
225 Point3D< Real > position;
226 Real weight;
227 Real value , _value;
228 SinglePointData operator + ( const SinglePointData& p ) const { return SinglePointData( position + p.position , value + p.value , weight + p.weight ); }
229 SinglePointData& operator += ( const SinglePointData& p ){ position += p.position ; weight += p.weight , value += p.value ; return *this; }
230 SinglePointData operator * ( Real s ) const { return SinglePointData( position*s , weight*s , value*s ); }
231 SinglePointData& operator *= ( Real s ){ position *= s , weight *= s , value *= s ; return *this; }
232 SinglePointData operator / ( Real s ) const { return SinglePointData( position/s , weight/s , value/s ); }
233 SinglePointData& operator /= ( Real s ){ position /= s , weight /= s , value /= s ; return *this; }
SinglePointDataSinglePointData234 SinglePointData( void ) : position( Point3D< Real >() ) , weight(0) , value(0) , _value(0) { ; }
SinglePointDataSinglePointData235 SinglePointData( Point3D< Real > p , Real v , Real w ) { position = p , value = v , weight = w , _value = (Real)0; }
236 };
237 template< class Real >
238 struct SinglePointData< Real , true > : public SinglePointData< Real , false >
239 {
240 using SinglePointData< Real , false >::position;
241 using SinglePointData< Real , false >::weight;
242 using SinglePointData< Real , false >::value;
243 using SinglePointData< Real , false >::_value;
244 Point3D< Real > gradient , _gradient;
245 SinglePointData operator + ( const SinglePointData& p ) const { return SinglePointData( position + p.position , weight + p.weight , value + p.value , gradient + p.gradient ); }
246 SinglePointData& operator += ( const SinglePointData& p ){ position += p.position , weight += p.weight , value += p.value , gradient += p.gradient ; return *this; }
247 SinglePointData operator * ( Real s ) const { return SinglePointData( position*s , weight*s , value*s , gradient*s ); }
248 SinglePointData& operator *= ( Real s ){ position *= s , weight *= s , value *= s , gradient *= s ; return *this; }
249 SinglePointData operator / ( Real s ) const { return SinglePointData( position/s , weight/s , value/s , gradient/s ); }
250 SinglePointData& operator /= ( Real s ){ position /= s , weight /= s , value /= s , gradient /= s ; return *this; }
251 SinglePointData( void ) : SinglePointData< Real , false >() , gradient( Point3D< Real >() ) , _gradient( Point3D< Real >() ) { ; }
252 SinglePointData( Point3D< Real > p , Real v , Point3D< Real > g , Real w ) : SinglePointData< Real , false >( p , v , w ) { gradient = g , _gradient = Point3D< Real >(); }
253 };
254
255 #if POINT_DATA_RES
256 template< class Real , bool HasGradients >
257 struct PointData
258 {
259 static const int RES = POINT_DATA_RES;
260 static const int SAMPLES = RES * RES * RES;
261
262 SinglePointData< Real , HasGradients > points[SAMPLES];
263 SinglePointData< Real , HasGradients >& operator[] ( int idx ) { return points[idx]; }
264 const SinglePointData< Real , HasGradients >& operator[] ( int idx ) const { return points[idx]; }
265
266 static void SetIndices( Point3D< Real > p , Point3D< Real > c , Real w , int x[3] )
267 {
268 for( int d=0 ; d<3 ; d++ ) x[d] = std::max< int >( 0 , std::min< int >( RES-1 , int( floor( ( p[d]-( c[d]-w/2 ) ) / w * RES ) ) ) );
269 }
270
271 void addPoint( SinglePointData< Real , HasGradients > p , Point3D< Real > center , Real width )
272 {
273 int x[3];
274 SetIndices( p.position , center , width , x );
275 points[ x[0]+x[1]*RES+x[2]*RES*RES ] += p;
276 }
277
278 PointData operator + ( const PointData& p ) const { PointData _p ; for( int c=0 ; c<SAMPLES ; c++ ) _p.points[c] = points[c] + _p.points[c] ; return _p; }
279 PointData& operator += ( const PointData& p ){ for( int c=0 ; c<SAMPLES ; c++ ) points[c] += p.points[c] ; return *this; }
280 PointData operator * ( Real s ) const { PointData _p ; for( int c=0 ; c<SAMPLES ; c++ ) _p.points[c] = points[c] * s ; return _p; }
281 PointData& operator *= ( Real s ){ for( int c=0 ; c<SAMPLES ; c++ ) points[c] *= s ; return *this; }
282 PointData operator / ( Real s ) const { PointData _p ; for( int c=0 ; c<SAMPLES ; c++ ) _p.points[c] = points[c] / s ; return _p; }
283 PointData& operator /= ( Real s ){ for( int c=0 ; c<SAMPLES ; c++ ) points[c] /= s ; return *this; }
284 };
285 #else // !POINT_DATA_RES
286 template< class Real , bool HasGradients > using PointData = SinglePointData< Real , HasGradients >;
287 #endif // POINT_DATA_RES
288
289 template< class Data , int Degree >
290 struct SparseNodeData
291 {
292 size_t size( void ) const { return _data.size(); }
293 const Data& operator[] ( int idx ) const { return _data[idx]; }
294 Data& operator[] ( int idx ) { return _data[idx]; }
295 void reserve( size_t sz ){ if( sz>_indices.size() ) _indices.resize( sz , -1 ); }
296 Data* operator()( const OctNode< TreeNodeData >* node ){ return ( node->nodeData.nodeIndex<0 || node->nodeData.nodeIndex>=(int)_indices.size() || _indices[ node->nodeData.nodeIndex ]<0 ) ? NULL : &_data[ _indices[ node->nodeData.nodeIndex ] ]; }
297 const Data* operator()( const OctNode< TreeNodeData >* node ) const { return ( node->nodeData.nodeIndex<0 || node->nodeData.nodeIndex>=(int)_indices.size() || _indices[ node->nodeData.nodeIndex ]<0 ) ? NULL : &_data[ _indices[ node->nodeData.nodeIndex ] ]; }
298 Data& operator[]( const OctNode< TreeNodeData >* node )
299 {
300 if( node->nodeData.nodeIndex>=(int)_indices.size() ) _indices.resize( node->nodeData.nodeIndex+1 , -1 );
301 if( _indices[ node->nodeData.nodeIndex ]==-1 )
302 {
303 _indices[ node->nodeData.nodeIndex ] = (int)_data.size();
304 _data.push_back( Data() );
305 }
306 return _data[ _indices[ node->nodeData.nodeIndex ] ];
307 }
308 void remapIndices( const std::vector< int >& map )
309 {
310 std::vector< int > temp = _indices;
311 _indices.resize( map.size() );
312 for( size_t i=0 ; i<map.size() ; i++ )
313 if( map[i]<(int)temp.size() ) _indices[i] = temp[ map[i] ];
314 else _indices[i] = -1;
315 }
316 template< class _Data , int _Degree > friend struct SparseNodeData;
317 template< class _Data , int _Degree >
318 void init( const SparseNodeData< _Data , _Degree >& snd ){ _indices = snd._indices , _data.resize( snd._data.size() ); }
319 void remove( const OctNode< TreeNodeData >* node ){ if( node->nodeData.nodeIndex<(int)_indices.size() && node->nodeData.nodeIndex>=0 ) _indices[ node->nodeData.nodeIndex ] = -1; }
320 protected:
321 std::vector< int > _indices;
322 std::vector< Data > _data;
323 };
324 template< class Data , int Degree >
325 struct DenseNodeData
326 {
327 DenseNodeData( void ){ _data = NullPointer( Data ) ; _sz = 0; }
328 DenseNodeData( size_t sz ){ _sz = sz ; if( sz ) _data = NewPointer< Data >( sz ) ; else _data = NullPointer( Data ); }
329 DenseNodeData( const DenseNodeData& d ) : DenseNodeData() { _resize( d._sz ) ; if( _sz ) memcpy( _data , d._data , sizeof(Data) * _sz ); }
330 DenseNodeData( DenseNodeData&& d ){ _data = d._data , _sz = d._sz ; d._data = NullPointer( Data ) , d._sz = 0; }
331 DenseNodeData& operator = ( const DenseNodeData& d ){ _resize( d._sz ) ; if( _sz ) memcpy( _data , d._data , sizeof(Data) * _sz ) ; return *this; }
332 DenseNodeData& operator = ( DenseNodeData&& d ){ size_t __sz = _sz ; Pointer( Data ) __data = _data ; _data = d._data , _sz = d._sz ; d._data = __data , d._sz = __sz ; return *this; }
333 ~DenseNodeData( void ){ DeletePointer( _data ) ; _sz = 0; }
334
335 Data& operator[] ( int idx ) { return _data[idx]; }
336 const Data& operator[] ( int idx ) const { return _data[idx]; }
337 size_t size( void ) const { return _sz; }
338 Data& operator[]( const OctNode< TreeNodeData >* node ) { return _data[ node->nodeData.nodeIndex ]; }
339 Data* operator()( const OctNode< TreeNodeData >* node ) { return ( node==NULL || node->nodeData.nodeIndex>=(int)_sz ) ? NULL : &_data[ node->nodeData.nodeIndex ]; }
340 const Data* operator()( const OctNode< TreeNodeData >* node ) const { return ( node==NULL || node->nodeData.nodeIndex>=(int)_sz ) ? NULL : &_data[ node->nodeData.nodeIndex ]; }
341 int index( const OctNode< TreeNodeData >* node ) const { return ( !node || node->nodeData.nodeIndex<0 || node->nodeData.nodeIndex>=(int)this->_data.size() ) ? -1 : node->nodeData.nodeIndex; }
342 protected:
343 size_t _sz;
344 void _resize( size_t sz ){ DeletePointer( _data ) ; if( sz ) _data = NewPointer< Data >( sz ) ; else _data = NullPointer( Data ) ; _sz = sz; }
345 Pointer( Data ) _data;
346 };
347
348 // This is may be necessary in case the memory usage is larger than what fits on the stack
349 template< class C , int N > struct Stencil
350 {
351 Stencil( void ){ _values = NewPointer< C >( N * N * N ); }
352 ~Stencil( void ){ DeletePointer( _values ); }
353 C& operator()( int i , int j , int k ){ return _values[ i*N*N + j*N + k ]; }
354 const C& operator()( int i , int j , int k ) const { return _values[ i*N*N + j*N + k ]; }
355 protected:
356 Pointer( C ) _values;
357 };
358
359 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 >
360 class SystemCoefficients
361 {
362 typedef typename BSplineIntegrationData< Degree1 , BType1 , Degree2 , BType2 >::FunctionIntegrator FunctionIntegrator;
363 static const int OverlapSize = BSplineOverlapSizes< Degree1 , Degree2 >::OverlapSize;
364 static const int OverlapStart = BSplineOverlapSizes< Degree1 , Degree2 >::OverlapStart;
365 static const int OverlapEnd = BSplineOverlapSizes< Degree1 , Degree2 >::OverlapEnd;
366 public:
367 typedef typename BSplineIntegrationData< Degree1 , BType1 , Degree2 , BType2 >::FunctionIntegrator::template Integrator< DERIVATIVES( Degree1 ) , DERIVATIVES( Degree2 ) > Integrator;
368 typedef typename BSplineIntegrationData< Degree1 , BType1 , Degree2 , BType2 >::FunctionIntegrator::template ChildIntegrator< DERIVATIVES( Degree1 ) , DERIVATIVES( Degree2 ) > ChildIntegrator;
369
370 // The FEMSystemFunctor is a class that takes an object of type Integrator/ChildIntegrator, as well as a pair of indices of octree nodes
371 // and returns the corresponding system coefficient.
372 template< class _FEMSystemFunctor > static void SetCentralSystemStencil ( const _FEMSystemFunctor& F , const Integrator& integrator , Stencil< double , OverlapSize >& stencil );
373 template< class _FEMSystemFunctor > static void SetCentralSystemStencils( const _FEMSystemFunctor& F , const ChildIntegrator& integrator , Stencil< double , OverlapSize > stencils[2][2][2] );
374 template< bool Reverse , class _FEMSystemFunctor > static void SetCentralConstraintStencil ( const _FEMSystemFunctor& F , const Integrator& integrator , Stencil< double , OverlapSize >& stencil );
375 template< bool Reverse , class _FEMSystemFunctor > static void SetCentralConstraintStencils( const _FEMSystemFunctor& F , const ChildIntegrator& integrator , Stencil< double , OverlapSize > stencils[2][2][2] );
376 template< bool Reverse , class _FEMSystemFunctor > static void SetCentralConstraintStencil ( const _FEMSystemFunctor& F , const Integrator& integrator , Stencil< Point3D< double > , OverlapSize >& stencil );
377 template< bool Reverse , class _FEMSystemFunctor > static void SetCentralConstraintStencils( const _FEMSystemFunctor& F , const ChildIntegrator& integrator , Stencil< Point3D< double > , OverlapSize > stencils[2][2][2] );
378 };
379
380 template< int FEMDegree , BoundaryType BType >
381 struct FEMSystemFunctor
382 {
383 double massWeight , lapWeight , biLapWeight;
384 FEMSystemFunctor( double mWeight=0 , double lWeight=0 , double bWeight=0 ) : massWeight( mWeight ) , lapWeight( lWeight ) , biLapWeight( bWeight ) { ; }
385 double integrate( const typename SystemCoefficients< FEMDegree , BType , FEMDegree , BType >:: Integrator& integrator , const int off1[] , const int off2[] ) const { return _integrate( integrator , off1 , off2 ); }
386 double integrate( const typename SystemCoefficients< FEMDegree , BType , FEMDegree , BType >::ChildIntegrator& integrator , const int off1[] , const int off2[] ) const { return _integrate( integrator , off1 , off2 ); }
387 bool vanishesOnConstants( void ) const { return massWeight==0; }
388 protected:
389 template< class I > double _integrate( const I& integrator , const int off1[] , const int off2[] ) const;
390 };
391 template< int SFDegree , BoundaryType SFBType , int FEMDegree , BoundaryType FEMBType >
392 struct FEMSFConstraintFunctor
393 {
394 double massWeight , lapWeight , biLapWeight;
395 FEMSFConstraintFunctor( double mWeight=0 , double lWeight=0 , double bWeight=0 ) : massWeight( mWeight ) , lapWeight( lWeight ) , biLapWeight( bWeight ) { ; }
396 template< bool Reverse >
397 double integrate( const typename SystemCoefficients< Reverse ? FEMDegree : SFDegree , Reverse ? FEMBType : SFBType , Reverse ? SFDegree : FEMDegree , Reverse ? SFBType : FEMBType >:: Integrator& integrator , const int off1[] , const int off2[] ) const { return _integrate< Reverse >( integrator , off1 , off2 ); }
398 template< bool Reverse >
399 double integrate( const typename SystemCoefficients< Reverse ? FEMDegree : SFDegree , Reverse ? FEMBType : SFBType , Reverse ? SFDegree : FEMDegree , Reverse ? SFBType : FEMBType >::ChildIntegrator& integrator , const int off1[] , const int off2[] ) const { return _integrate< Reverse >( integrator , off1 , off2 ); }
400 protected:
401 template< bool Reverse , class I > double _integrate( const I& integrator , const int off1[] , const int off[2] ) const;
402 };
403 template< int VFDegree , BoundaryType VFBType , int FEMDegree , BoundaryType FEMBType >
404 struct FEMVFConstraintFunctor
405 {
406 double lapWeight , biLapWeight;
407 FEMVFConstraintFunctor( double lWeight=0 , double bWeight=0 ) : lapWeight( lWeight ) , biLapWeight( bWeight ) { ; }
408 template< bool Reverse >
409 Point3D< double > integrate( const typename SystemCoefficients< Reverse ? FEMDegree : VFDegree , Reverse ? FEMBType : VFBType , Reverse ? VFDegree : FEMDegree , Reverse ? VFBType : FEMBType >:: Integrator& integrator , const int off1[] , const int off2[] ) const { return _integrate< Reverse >( integrator , off1 , off2 ); }
410 template< bool Reverse >
411 Point3D< double > integrate( const typename SystemCoefficients< Reverse ? FEMDegree : VFDegree , Reverse ? FEMBType : VFBType , Reverse ? VFDegree : FEMDegree , Reverse ? VFBType : FEMBType >::ChildIntegrator& integrator , const int off1[] , const int off2[] ) const { return _integrate< Reverse >( integrator , off1 , off2 ); }
412 protected:
413 template< bool Reverse , class I > Point3D< double > _integrate( const I& integrator , const int off1[] , const int off[2] ) const;
414 };
415
416 inline void SetGhostFlag( OctNode< TreeNodeData >* node , bool flag ){ if( node && node->parent ) node->parent->nodeData.setGhostFlag( flag ); }
417 inline bool GetGhostFlag( const OctNode< TreeNodeData >* node ){ return node==NULL || node->parent==NULL || node->parent->nodeData.getGhostFlag( ); }
418 inline bool IsActiveNode( const OctNode< TreeNodeData >* node ){ return !GetGhostFlag( node ); }
419
420 template< class Real >
421 class Octree
422 {
423 typedef OctNode< TreeNodeData > TreeOctNode;
424 static int _NodeCount;
425 static void _NodeInitializer( TreeOctNode& node ){ node.nodeData.nodeIndex = _NodeCount++; }
426 public:
427 #if 0
428 struct LocalDepth
429 {
430 LocalDepth( int d=0 ) : _d(d) { ; }
431 operator int&() { return _d; }
432 operator int () const { return _d; }
433 protected:
434 int _d;
435 };
436 struct LocalOffset
437 {
438 LocalOffset( const int* off=NULL ){ if( off ) memcpy( _off , off , sizeof(_off) ) ; else memset( _off , 0 , sizeof( _off ) ); }
439 operator int*() { return _off; }
440 operator const int*() const { return _off; }
441 protected:
442 int _off[3];
443 };
444 #else
445 typedef int LocalDepth;
446 typedef int LocalOffset[3];
447 #endif
448
449 static void ResetNodeCount( void ){ _NodeCount = 0 ; }
450 static int NodeCount( void ){ return _NodeCount; }
451 template< int FEMDegree , BoundaryType BType > void functionIndex( const TreeOctNode* node , int idx[3] ) const;
452
453 struct PointSample{ const TreeOctNode* node ; ProjectiveData< OrientedPoint3D< Real > , Real > sample; };
454
455 typedef typename TreeOctNode:: NeighborKey< 1 , 1 > AdjacenctNodeKey;
456 typedef typename TreeOctNode::ConstNeighborKey< 1 , 1 > ConstAdjacenctNodeKey;
457
458 template< int FEMDegree , BoundaryType BType > bool isValidFEMNode( const TreeOctNode* node ) const;
459 bool isValidSpaceNode( const TreeOctNode* node ) const;
460 TreeOctNode* leaf( Point3D< Real > p );
461 const TreeOctNode* leaf( Point3D< Real > p ) const;
462
463 template< bool HasGradients >
464 struct InterpolationInfo
465 {
466 SparseNodeData< PointData< Real , HasGradients > , 0 > iData;
467 Real valueWeight , gradientWeight;
468 InterpolationInfo( const class Octree< Real >& tree , const std::vector< PointSample >& samples , Real pointValue , int adaptiveExponent , Real v , Real g ) : valueWeight(v) , gradientWeight(g)
469 { iData = tree._densifyInterpolationInfo< HasGradients >( samples , pointValue , adaptiveExponent ); }
470 PointData< Real , HasGradients >* operator()( const OctNode< TreeNodeData >* node ){ return iData(node); }
471 const PointData< Real , HasGradients >* operator()( const OctNode< TreeNodeData >* node ) const { return iData(node); }
472 };
473
474 template< int DensityDegree > struct DensityEstimator : public SparseNodeData< Real , DensityDegree >
475 {
476 DensityEstimator( int kernelDepth ) : _kernelDepth( kernelDepth ){ ; }
477 int kernelDepth( void ) const { return _kernelDepth; }
478 protected:
479 int _kernelDepth;
480 };
481 protected:
482 bool _isValidSpaceNode( const TreeOctNode* node ) const { return !GetGhostFlag( node ) && ( node->nodeData.flags & TreeNodeData::SPACE_FLAG ); }
483 bool _isValidFEMNode( const TreeOctNode* node ) const { return !GetGhostFlag( node ) && ( node->nodeData.flags & TreeNodeData::FEM_FLAG ); }
484
485 TreeOctNode* _tree;
486 TreeOctNode* _spaceRoot;
487 SortedTreeNodes _sNodes;
488 LocalDepth _fullDepth , _maxDepth;
489
490 static bool _InBounds( Point3D< Real > p );
491
492 int _depthOffset;
493 int _localToGlobal( LocalDepth d ) const { return d + _depthOffset; }
494 LocalDepth _localDepth( const TreeOctNode* node ) const { return node->depth() - _depthOffset; }
495 LocalDepth _localMaxDepth( const TreeOctNode* tree ) const { return tree->maxDepth() - _depthOffset; }
496 int _localInset( LocalDepth d ) const { return _depthOffset<=1 ? 0 : 1<<( d + _depthOffset - 1 ); }
497 void _localDepthAndOffset( const TreeOctNode* node , LocalDepth& d , LocalOffset& off ) const
498 {
499 node->depthAndOffset( d , off ) ; d -= _depthOffset;
500 int inset = _localInset( d );
501 off[0] -= inset , off[1] -= inset , off[2] -= inset;
502 }
503 template< int FEMDegree , BoundaryType BType > static int _BSplineBegin( LocalDepth depth ){ return BSplineEvaluationData< FEMDegree , BType >::Begin( depth ); }
504 template< int FEMDegree , BoundaryType BType > static int _BSplineEnd ( LocalDepth depth ){ return BSplineEvaluationData< FEMDegree , BType >::End ( depth ); }
505 template< int FEMDegree , BoundaryType BType >
506 bool _outOfBounds( const TreeOctNode* node ) const
507 {
508 if( !node ) return true;
509 LocalDepth d ; LocalOffset off;
510 _localDepthAndOffset( node , d , off );
511 return d<0 || BSplineEvaluationData< FEMDegree , BType >::OutOfBounds( d , off[0] ) || BSplineEvaluationData< FEMDegree , BType >::OutOfBounds( d , off[1] ) || BSplineEvaluationData< FEMDegree , BType >::OutOfBounds( d , off[2] );
512 }
513 int _sNodesBegin( LocalDepth d ) const { return _sNodes.begin( _localToGlobal( d ) ); }
514 int _sNodesEnd ( LocalDepth d ) const { return _sNodes.end ( _localToGlobal( d ) ); }
515 int _sNodesSize ( LocalDepth d ) const { return _sNodes.size ( _localToGlobal( d ) ); }
516 int _sNodesBegin( LocalDepth d , int slice ) const { return _sNodes.begin( _localToGlobal( d ) , slice + _localInset( d ) ); }
517 int _sNodesEnd ( LocalDepth d , int slice ) const { return _sNodes.end ( _localToGlobal( d ) , slice + _localInset( d ) ); }
518 int _sNodesSize ( LocalDepth d , int slice ) const { return _sNodes.size ( _localToGlobal( d ) , slice + _localInset( d ) ); }
519
520 template< int FEMDegree > static bool _IsInteriorlySupported( LocalDepth depth , const LocalOffset off )
521 {
522 if( depth>=0 )
523 {
524 int begin , end;
525 BSplineSupportSizes< FEMDegree >::InteriorSupportedSpan( depth , begin , end );
526 return ( off[0]>=begin && off[0]<end && off[1]>=begin && off[1]<end && off[2]>=begin && off[2]<end );
527 }
528 else return false;
529 }
530 template< int FEMDegree > bool _isInteriorlySupported( const TreeOctNode* node ) const
531 {
532 if( !node ) return false;
533 LocalDepth d ; LocalOffset off;
534 _localDepthAndOffset( node , d , off );
535 return _IsInteriorlySupported< FEMDegree >( d , off );
536 }
537 template< int FEMDegree1 , int FEMDegree2 > static bool _IsInteriorlyOverlapped( LocalDepth depth , const LocalOffset off )
538 {
539 if( depth>=0 )
540 {
541 int begin , end;
542 BSplineIntegrationData< FEMDegree1 , BOUNDARY_NEUMANN , FEMDegree2 , BOUNDARY_NEUMANN >::InteriorOverlappedSpan( depth , begin , end );
543 return ( off[0]>=begin && off[0]<end && off[1]>=begin && off[1]<end && off[2]>=begin && off[2]<end );
544 }
545 else return false;
546 }
547 template< int FEMDegree1 , int FEMDegree2 > bool _isInteriorlyOverlapped( const TreeOctNode* node ) const
548 {
549 if( !node ) return false;
550 LocalDepth d ; LocalOffset off;
551 _localDepthAndOffset( node , d , off );
552 return _IsInteriorlyOverlapped< FEMDegree1 , FEMDegree2 >( d , off );
553 }
554 void _startAndWidth( const TreeOctNode* node , Point3D< Real >& start , Real& width ) const
555 {
556 LocalDepth d ; LocalOffset off;
557 _localDepthAndOffset( node , d , off );
558 if( d>=0 ) width = Real( 1.0 / (1<< d ) );
559 else width = Real( 1.0 * (1<<(-d)) );
560 for( int dd=0 ; dd<DIMENSION ; dd++ ) start[dd] = Real( off[dd] ) * width;
561 }
562 void _centerAndWidth( const TreeOctNode* node , Point3D< Real >& center , Real& width ) const
563 {
564 int d , off[3];
565 _localDepthAndOffset( node , d , off );
566 width = Real( 1.0 / (1<<d) );
567 for( int dd=0 ; dd<DIMENSION ; dd++ ) center[dd] = Real( off[dd] + 0.5 ) * width;
568 }
569 int _childIndex( const TreeOctNode* node , Point3D< Real > p ) const
570 {
571 Point3D< Real > c ; Real w;
572 _centerAndWidth( node , c , w );
573 return ( p[0]<c[0] ? 0 : 1 ) | ( p[1]<c[1] ? 0 : 2 ) | ( p[2]<c[2] ? 0 : 4 );
574 }
575
576 template< int Degree , BoundaryType BType > void _setFullDepth( TreeOctNode* node , LocalDepth depth ) const;
577 template< int Degree , BoundaryType BType > void _setFullDepth( LocalDepth depth );
578
579 template< int LeftRadius , int RightRadius >
580 static typename TreeOctNode::ConstNeighbors< LeftRadius + RightRadius + 1 >& _neighbors( TreeOctNode::ConstNeighborKey< LeftRadius , RightRadius >& key , const TreeOctNode* node ){ return key.neighbors[ node->depth() ]; }
581 template< int LeftRadius , int RightRadius >
582 static typename TreeOctNode::Neighbors< LeftRadius + RightRadius + 1 >& _neighbors( TreeOctNode::NeighborKey< LeftRadius , RightRadius >& key , const TreeOctNode* node ){ return key.neighbors[ node->depth() ]; }
583 template< int LeftRadius , int RightRadius >
584 static const typename TreeOctNode::template Neighbors< LeftRadius + RightRadius + 1 >& _neighbors( const typename TreeOctNode::template NeighborKey< LeftRadius , RightRadius >& key , const TreeOctNode* node ){ return key.neighbors[ node->depth() ]; }
585 template< int LeftRadius , int RightRadius >
586 static const typename TreeOctNode::template ConstNeighbors< LeftRadius + RightRadius + 1 >& _neighbors( const typename TreeOctNode::template ConstNeighborKey< LeftRadius , RightRadius >& key , const TreeOctNode* node ){ return key.neighbors[ node->depth() ]; }
587
588 public:
589 LocalDepth depth( const TreeOctNode* node ) const { return _localDepth( node ); }
590 void depthAndOffset( const TreeOctNode* node , LocalDepth& depth , LocalOffset& offset ) const { _localDepthAndOffset( node , depth , offset ); }
591
592 int nodesBegin( LocalDepth d ) const { return _sNodes.begin( _localToGlobal( d ) ); }
593 int nodesEnd ( LocalDepth d ) const { return _sNodes.end ( _localToGlobal( d ) ); }
594 int nodesSize ( LocalDepth d ) const { return _sNodes.size ( _localToGlobal( d ) ); }
595 int nodesBegin( LocalDepth d , int slice ) const { return _sNodes.begin( _localToGlobal( d ) , slice + _localInset( d ) ); }
596 int nodesEnd ( LocalDepth d , int slice ) const { return _sNodes.end ( _localToGlobal( d ) , slice + _localInset( d ) ); }
597 int nodesSize ( LocalDepth d , int slice ) const { return _sNodes.size ( _localToGlobal( d ) , slice + _localInset( d ) ); }
598 const TreeOctNode* node( int idx ) const { return _sNodes.treeNodes[idx]; }
599 protected:
600
601 ////////////////////////////////////
602 // System construction code //
603 // MultiGridOctreeData.System.inl //
604 ////////////////////////////////////
605 template< int FEMDegree >
606 void _setMultiColorIndices( int start , int end , std::vector< std::vector< int > >& indices ) const;
607 struct _SolverStats
608 {
609 double evaluateTime , systemTime , solveTime;
610 double bNorm2 , inRNorm2 , outRNorm2;
611 };
612 template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients >
613 int _solveSystemGS( const FEMSystemFunctor& F , const BSplineData< FEMDegree , BType >& bsData , InterpolationInfo< HasGradients >* interpolationInfo , LocalDepth depth , DenseNodeData< Real , FEMDegree >& solution , DenseNodeData< Real , FEMDegree >& constraints , DenseNodeData< Real , FEMDegree >& metSolutionConstraints , int iters , bool coarseToFine , _SolverStats& stats , bool computeNorms );
614 template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients >
615 int _solveSystemCG( const FEMSystemFunctor& F , const BSplineData< FEMDegree , BType >& bsData , InterpolationInfo< HasGradients >* interpolationInfo , LocalDepth depth , DenseNodeData< Real , FEMDegree >& solution , DenseNodeData< Real , FEMDegree >& constraints , DenseNodeData< Real , FEMDegree >& metSolutionConstraints , int iters , bool coarseToFine , _SolverStats& stats , bool computeNorms , double accuracy );
616 template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients >
617 int _setMatrixRow( const FEMSystemFunctor& F , const InterpolationInfo< HasGradients >* interpolationInfo , const typename TreeOctNode::Neighbors< BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize >& neighbors , Pointer( MatrixEntry< Real > ) row , int offset , const typename BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::FunctionIntegrator::template Integrator< DERIVATIVES( FEMDegree ) , DERIVATIVES( FEMDegree ) >& integrator , const Stencil< double , BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize >& stencil , const BSplineData< FEMDegree , BType >& bsData ) const;
618 template< int FEMDegree , BoundaryType BType >
619 int _getMatrixRowSize( const typename TreeOctNode::Neighbors< BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize >& neighbors ) const;
620
621 template< int FEMDegree1 , int FEMDegree2 > static void _SetParentOverlapBounds( const TreeOctNode* node , int& startX , int& endX , int& startY , int& endY , int& startZ , int& endZ );
622 // Updates the constraints @(depth) based on the solution coefficients @(depth-1)
623
624 template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients >
625 void _updateConstraintsFromCoarser( const FEMSystemFunctor& F , const InterpolationInfo< HasGradients >* interpolationInfo , const typename TreeOctNode::Neighbors< BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize >& neighbors , const typename TreeOctNode::Neighbors< BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize >& pNeighbors , TreeOctNode* node , DenseNodeData< Real , FEMDegree >& constraints , const DenseNodeData< Real , FEMDegree >& metSolution , const typename BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::FunctionIntegrator::template ChildIntegrator< DERIVATIVES( FEMDegree ) , DERIVATIVES( FEMDegree ) >& childIntegrator , const Stencil< double , BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize >& stencil , const BSplineData< FEMDegree , BType >& bsData ) const;
626
627 // evaluate the points @(depth) using coefficients @(depth-1)
628 template< int FEMDegree , BoundaryType BType , bool HasGradients >
629 void _setPointValuesFromCoarser( InterpolationInfo< HasGradients >& interpolationInfo , LocalDepth highDepth , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& upSampledCoefficients );
630
631 // Updates the cumulative integral constraints @(depth-1) based on the change in solution coefficients @(depth)
632 template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor >
633 void _updateCumulativeIntegralConstraintsFromFiner( const FEMSystemFunctor& F ,
634 const BSplineData< FEMDegree , BType >& bsData , LocalDepth highDepth , const DenseNodeData< Real , FEMDegree >& fineSolution , DenseNodeData< Real , FEMDegree >& cumulativeConstraints ) const;
635 // Updates the cumulative interpolation constraints @(depth-1) based on the change in solution coefficient @(depth)
636 template< int FEMDegree , BoundaryType BType , bool HasGradients >
637 void _updateCumulativeInterpolationConstraintsFromFiner( const InterpolationInfo< HasGradients >& interpolationInfo ,
638 const BSplineData< FEMDegree , BType >& bsData , LocalDepth highDepth , const DenseNodeData< Real , FEMDegree >& fineSolution , DenseNodeData< Real , FEMDegree >& cumulativeConstraints ) const;
639
640 template< int FEMDegree , BoundaryType BType >
641 Real _coarserFunctionValue( Point3D< Real > p , const PointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& upSampledCoefficients ) const;
642 template< int FEMDegree , BoundaryType BType >
643 Point3D< Real > _coarserFunctionGradient( Point3D< Real > p , const PointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& upSampledCoefficients ) const;
644 template< int FEMDegree , BoundaryType BType >
645 Real _finerFunctionValue( Point3D< Real > p , const PointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& coefficients ) const;
646 template< int FEMDegree , BoundaryType BType >
647 Point3D< Real > _finerFunctionGradient( Point3D< Real > p , const PointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& coefficients ) const;
648 template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients >
649 int _getSliceMatrixAndUpdateConstraints( const FEMSystemFunctor& F , const InterpolationInfo< HasGradients >* interpolationInfo , SparseMatrix< Real >& matrix , DenseNodeData< Real , FEMDegree >& constraints , typename BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::FunctionIntegrator::template Integrator< DERIVATIVES( FEMDegree ) , DERIVATIVES( FEMDegree ) >& integrator , typename BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::FunctionIntegrator::template ChildIntegrator< DERIVATIVES( FEMDegree ) , DERIVATIVES( FEMDegree ) >& childIntegrator , const BSplineData< FEMDegree , BType >& bsData , LocalDepth depth , int slice , const DenseNodeData< Real , FEMDegree >& metSolution , bool coarseToFine );
650 template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients >
651 int _getMatrixAndUpdateConstraints( const FEMSystemFunctor& F , const InterpolationInfo< HasGradients >* interpolationInfo , SparseMatrix< Real >& matrix , DenseNodeData< Real , FEMDegree >& constraints , typename BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::FunctionIntegrator::template Integrator< DERIVATIVES( FEMDegree ) , DERIVATIVES( FEMDegree ) >& integrator , typename BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::FunctionIntegrator::template ChildIntegrator< DERIVATIVES( FEMDegree ) , DERIVATIVES( FEMDegree ) >& childIntegrator , const BSplineData< FEMDegree , BType >& bsData , LocalDepth depth , const DenseNodeData< Real , FEMDegree >& metSolution , bool coarseToFine );
652
653 // Down samples constraints @(depth) to constraints @(depth-1)
654 template< class C , int FEMDegree , BoundaryType BType > void _downSample( LocalDepth highDepth , DenseNodeData< C , FEMDegree >& constraints ) const;
655 // Up samples coefficients @(depth-1) to coefficients @(depth)
656 template< class C , int FEMDegree , BoundaryType BType > void _upSample( LocalDepth highDepth , DenseNodeData< C , FEMDegree >& coefficients ) const;
657 template< class C , int FEMDegree , BoundaryType BType > static void _UpSample( LocalDepth highDepth , ConstPointer( C ) lowCoefficients , Pointer( C ) highCoefficients , int threads );
658 public:
659 template< class C , int FEMDegree , BoundaryType BType > DenseNodeData< C , FEMDegree > coarseCoefficients( const DenseNodeData< C , FEMDegree >& coefficients ) const;
660 template< class C , int FEMDegree , BoundaryType BType > DenseNodeData< C , FEMDegree > coarseCoefficients( const SparseNodeData< C , FEMDegree >& coefficients ) const;
661 protected:
662
663 /////////////////////////////////////////////
664 // Code for splatting point-sample data //
665 // MultiGridOctreeData.WeightedSamples.inl //
666 /////////////////////////////////////////////
667 template< int WeightDegree >
668 void _addWeightContribution( DensityEstimator< WeightDegree >& densityWeights , TreeOctNode* node , Point3D< Real > position , PointSupportKey< WeightDegree >& weightKey , Real weight=Real(1.0) );
669 template< int WeightDegree , class PointSupportKey >
670 Real _getSamplesPerNode( const DensityEstimator< WeightDegree >& densityWeights , const TreeOctNode* node , Point3D< Real > position , PointSupportKey& weightKey ) const;
671 template< int WeightDegree , class PointSupportKey >
672 void _getSampleDepthAndWeight( const DensityEstimator< WeightDegree >& densityWeights , const TreeOctNode* node , Point3D< Real > position , PointSupportKey& weightKey , Real& depth , Real& weight ) const;
673 template< int WeightDegree , class PointSupportKey >
674 void _getSampleDepthAndWeight( const DensityEstimator< WeightDegree >& densityWeights , Point3D< Real > position , PointSupportKey& weightKey , Real& depth , Real& weight ) const;
675 template< bool CreateNodes , int DataDegree , class V > void _splatPointData( TreeOctNode* node , Point3D< Real > point , V v , SparseNodeData< V , DataDegree >& data , PointSupportKey< DataDegree >& dataKey );
676 template< bool CreateNodes , int WeightDegree , int DataDegree , class V > Real _splatPointData( const DensityEstimator< WeightDegree >& densityWeights , Point3D< Real > point , V v , SparseNodeData< V , DataDegree >& data , PointSupportKey< WeightDegree >& weightKey , PointSupportKey< DataDegree >& dataKey , LocalDepth minDepth , LocalDepth maxDepth , int dim=DIMENSION );
677 template< bool CreateNodes , int WeightDegree , int DataDegree , class V > Real _multiSplatPointData( const DensityEstimator< WeightDegree >* densityWeights , TreeOctNode* node , Point3D< Real > point , V v , SparseNodeData< V , DataDegree >& data , PointSupportKey< WeightDegree >& weightKey , PointSupportKey< DataDegree >& dataKey , int dim=DIMENSION );
678 template< class V , int DataDegree , BoundaryType BType , class Coefficients > V _evaluate( const Coefficients& coefficients , Point3D< Real > p , const BSplineData< DataDegree , BType >& bsData , const ConstPointSupportKey< DataDegree >& dataKey ) const;
679 public:
680 template< class V , int DataDegree , BoundaryType BType > Pointer( V ) voxelEvaluate( const DenseNodeData< V , DataDegree >& coefficients , int& res , Real isoValue=0.f , LocalDepth depth=-1 , bool primal=false );
681
682 template< int NormalDegree >
683 struct HasNormalDataFunctor
684 {
685 const SparseNodeData< Point3D< Real > , NormalDegree >& normalInfo;
686 HasNormalDataFunctor( const SparseNodeData< Point3D< Real > , NormalDegree >& ni ) : normalInfo( ni ){ ; }
687 bool operator() ( const TreeOctNode* node ) const
688 {
689 const Point3D< Real >* n = normalInfo( node );
690 if( n )
691 {
692 const Point3D< Real >& normal = *n;
693 if( normal[0]!=0 || normal[1]!=0 || normal[2]!=0 ) return true;
694 }
695 if( node->children ) for( int c=0 ; c<Cube::CORNERS ; c++ ) if( (*this)( node->children + c ) ) return true;
696 return false;
697 }
698 };
699 struct TrivialHasDataFunctor{ bool operator() ( const TreeOctNode* node ) const{ return true; } };
700
701 // [NOTE] The input/output for this method is pre-scaled by weight
702 template< bool HasGradients > bool _setInterpolationInfoFromChildren( TreeOctNode* node , SparseNodeData< PointData< Real , HasGradients > , 0 >& iInfo ) const;
703 template< bool HasGradients > SparseNodeData< PointData< Real , HasGradients > , 0 > _densifyInterpolationInfo( const std::vector< PointSample >& samples , Real pointValue , int adaptiveExponent ) const;
704
705 template< int FEMDegree , BoundaryType BType > void _setValidityFlags( void );
706 template< class HasDataFunctor > void _clipTree( const HasDataFunctor& f );
707
708 template< int FEMDegree , BoundaryType BType > SparseNodeData< Real , 0 > leafValues ( const DenseNodeData< Real , FEMDegree >& coefficients ) const;
709 template< int FEMDegree , BoundaryType BType > SparseNodeData< Point3D< Real > , 0 > leafGradients( const DenseNodeData< Real , FEMDegree >& coefficients ) const;
710
711 ////////////////////////////////////
712 // Evaluation Methods //
713 // MultiGridOctreeData.Evaluation //
714 ////////////////////////////////////
715 static const int CHILDREN = Cube::CORNERS;
716 template< int FEMDegree , BoundaryType BType >
717 struct _Evaluator
718 {
719 typename BSplineEvaluationData< FEMDegree , BType >::Evaluator evaluator;
720 typename BSplineEvaluationData< FEMDegree , BType >::ChildEvaluator childEvaluator;
721 Stencil< double , BSplineSupportSizes< FEMDegree >::SupportSize > cellStencil;
722 Stencil< double , BSplineSupportSizes< FEMDegree >::SupportSize > cellStencils [CHILDREN];
723 Stencil< double , BSplineSupportSizes< FEMDegree >::SupportSize > edgeStencil [Cube::EDGES ];
724 Stencil< double , BSplineSupportSizes< FEMDegree >::SupportSize > edgeStencils [CHILDREN][Cube::EDGES ];
725 Stencil< double , BSplineSupportSizes< FEMDegree >::SupportSize > faceStencil [Cube::FACES ];
726 Stencil< double , BSplineSupportSizes< FEMDegree >::SupportSize > faceStencils [CHILDREN][Cube::FACES ];
727 Stencil< double , BSplineSupportSizes< FEMDegree >::SupportSize > cornerStencil [Cube::CORNERS];
728 Stencil< double , BSplineSupportSizes< FEMDegree >::SupportSize > cornerStencils[CHILDREN][Cube::CORNERS];
729
730 Stencil< Point3D< double > , BSplineSupportSizes< FEMDegree >::SupportSize > dCellStencil;
731 Stencil< Point3D< double > , BSplineSupportSizes< FEMDegree >::SupportSize > dCellStencils [CHILDREN];
732 Stencil< Point3D< double > , BSplineSupportSizes< FEMDegree >::SupportSize > dEdgeStencil [Cube::EDGES ];
733 Stencil< Point3D< double > , BSplineSupportSizes< FEMDegree >::SupportSize > dEdgeStencils [CHILDREN][Cube::EDGES ];
734 Stencil< Point3D< double > , BSplineSupportSizes< FEMDegree >::SupportSize > dFaceStencil [Cube::FACES ];
735 Stencil< Point3D< double > , BSplineSupportSizes< FEMDegree >::SupportSize > dFaceStencils [CHILDREN][Cube::FACES ];
736 Stencil< Point3D< double > , BSplineSupportSizes< FEMDegree >::SupportSize > dCornerStencil [Cube::CORNERS];
737 Stencil< Point3D< double > , BSplineSupportSizes< FEMDegree >::SupportSize > dCornerStencils[CHILDREN][Cube::CORNERS];
738
739 void set( LocalDepth depth );
740 _Evaluator( void ){ _bsData = NULL; }
741 ~_Evaluator( void ){ if( _bsData ) delete _bsData , _bsData = NULL; }
742 protected:
743 BSplineData< FEMDegree , BType >* _bsData;
744 friend Octree;
745 };
746 template< class V , int FEMDegree , BoundaryType BType >
747 V _getCenterValue( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , const DenseNodeData< V , FEMDegree >& solution , const DenseNodeData< V , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const;
748 template< class V , int FEMDegree , BoundaryType BType >
749 V _getCornerValue( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , int corner , const DenseNodeData< V , FEMDegree >& solution , const DenseNodeData< V , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const;
750 template< class V , int FEMDegree , BoundaryType BType >
751 V _getEdgeValue ( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , int edge , const DenseNodeData< V , FEMDegree >& solution , const DenseNodeData< V , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const;
752 template< class V , int FEMDegree , BoundaryType BType >
753 V _getValue ( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , Point3D< Real > p , const DenseNodeData< V , FEMDegree >& solution , const DenseNodeData< V , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator ) const;
754
755 template< int FEMDegree , BoundaryType BType >
756 std::pair< Real , Point3D< Real > > _getCenterValueAndGradient( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const;
757 template< int FEMDegree , BoundaryType BType >
758 std::pair< Real , Point3D< Real > > _getCornerValueAndGradient( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , int corner , const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const;
759 template< int FEMDegree , BoundaryType BType >
760 std::pair< Real , Point3D< Real > > _getEdgeValueAndGradient ( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , int edge , const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const;
761 template< int FEMDegree , BoundaryType BType >
762 std::pair< Real , Point3D< Real > > _getValueAndGradient ( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , Point3D< Real > p , const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator ) const;
763
764 public:
765 template< int Degree , BoundaryType BType >
766 class MultiThreadedEvaluator
767 {
768 const Octree* _tree;
769 int _threads;
770 std::vector< ConstPointSupportKey< Degree > > _neighborKeys;
771 _Evaluator< Degree , BType > _evaluator;
772 const DenseNodeData< Real , Degree >& _coefficients;
773 DenseNodeData< Real , Degree > _coarseCoefficients;
774 public:
775 MultiThreadedEvaluator( const Octree* tree , const DenseNodeData< Real , Degree >& coefficients , int threads=1 );
776 Real value( Point3D< Real > p , int thread=0 , const TreeOctNode* node=NULL );
777 std::pair< Real , Point3D< Real > > valueAndGradient( Point3D< Real > , int thread=0 , const TreeOctNode* node=NULL );
778 };
779
780 ////////////////////////////////////////
781 // Iso-Surfacing Methods //
782 // MultiGridOctreeData.IsoSurface.inl //
783 ////////////////////////////////////////
784 protected:
785 struct _IsoEdge
786 {
787 long long edges[2];
788 _IsoEdge( void ){ edges[0] = edges[1] = 0; }
789 _IsoEdge( long long v1 , long long v2 ){ edges[0] = v1 , edges[1] = v2; }
790 long long& operator[]( int idx ){ return edges[idx]; }
791 const long long& operator[]( int idx ) const { return edges[idx]; }
792 };
793 struct _FaceEdges
794 {
795 _IsoEdge edges[2];
796 int count;
797 };
798 template< class Vertex >
799 struct _SliceValues
800 {
801 typename SortedTreeNodes::SliceTableData sliceData;
802 Pointer( Real ) cornerValues ; Pointer( Point3D< Real > ) cornerGradients ; Pointer( char ) cornerSet;
803 Pointer( long long ) edgeKeys ; Pointer( char ) edgeSet;
804 Pointer( _FaceEdges ) faceEdges ; Pointer( char ) faceSet;
805 Pointer( char ) mcIndices;
806 std::unordered_map< long long, std::vector< _IsoEdge > > faceEdgeMap;
807 std::unordered_map< long long, std::pair< int, Vertex > > edgeVertexMap;
808 std::unordered_map< long long, long long > vertexPairMap;
809
810 _SliceValues( void );
811 ~_SliceValues( void );
812 void reset( bool nonLinearFit );
813 protected:
814 int _oldCCount , _oldECount , _oldFCount , _oldNCount;
815 };
816 template< class Vertex >
817 struct _XSliceValues
818 {
819 typename SortedTreeNodes::XSliceTableData xSliceData;
820 Pointer( long long ) edgeKeys ; Pointer( char ) edgeSet;
821 Pointer( _FaceEdges ) faceEdges ; Pointer( char ) faceSet;
822 std::unordered_map< long long, std::vector< _IsoEdge > > faceEdgeMap;
823 std::unordered_map< long long, std::pair< int, Vertex > > edgeVertexMap;
824 std::unordered_map< long long, long long > vertexPairMap;
825
826 _XSliceValues( void );
827 ~_XSliceValues( void );
828 void reset( void );
829 protected:
830 int _oldECount , _oldFCount;
831 };
832 template< class Vertex >
833 struct _SlabValues
834 {
835 protected:
836 _XSliceValues< Vertex > _xSliceValues[2];
837 _SliceValues< Vertex > _sliceValues[2];
838 public:
839 _SliceValues< Vertex >& sliceValues( int idx ){ return _sliceValues[idx&1]; }
840 const _SliceValues< Vertex >& sliceValues( int idx ) const { return _sliceValues[idx&1]; }
841 _XSliceValues< Vertex >& xSliceValues( int idx ){ return _xSliceValues[idx&1]; }
842 const _XSliceValues< Vertex >& xSliceValues( int idx ) const { return _xSliceValues[idx&1]; }
843 };
844 template< class Vertex , int FEMDegree , BoundaryType BType >
845 void _setSliceIsoCorners( const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , Real isoValue , LocalDepth depth , int slice , std::vector< _SlabValues< Vertex > >& sValues , const _Evaluator< FEMDegree , BType >& evaluator , int threads );
846 template< class Vertex , int FEMDegree , BoundaryType BType >
847 void _setSliceIsoCorners( const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , Real isoValue , LocalDepth depth , int slice , int z , std::vector< _SlabValues< Vertex > >& sValues , const _Evaluator< FEMDegree , BType >& evaluator , int threads );
848 template< int WeightDegree , int ColorDegree , BoundaryType BType , class Vertex >
849 void _setSliceIsoVertices( const BSplineData< ColorDegree , BType >* colorBSData , const DensityEstimator< WeightDegree >* densityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > , Real > , ColorDegree >* colorData , Real isoValue , LocalDepth depth , int slice , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< _SlabValues< Vertex > >& sValues , int threads );
850 template< int WeightDegree , int ColorDegree , BoundaryType BType , class Vertex >
851 void _setSliceIsoVertices( const BSplineData< ColorDegree , BType >* colorBSData , const DensityEstimator< WeightDegree >* densityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > , Real > , ColorDegree >* colorData , Real isoValue , LocalDepth depth , int slice , int z , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< _SlabValues< Vertex > >& sValues , int threads );
852 template< int WeightDegree , int ColorDegree , BoundaryType BType , class Vertex >
853 void _setXSliceIsoVertices( const BSplineData< ColorDegree , BType >* colorBSData , const DensityEstimator< WeightDegree >* densityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > , Real > , ColorDegree >* colorData , Real isoValue , LocalDepth depth , int slab , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< _SlabValues< Vertex > >& sValues , int threads );
854 template< class Vertex >
855 void _setSliceIsoEdges( LocalDepth depth , int slice , std::vector< _SlabValues< Vertex > >& slabValues , int threads );
856 template< class Vertex >
857 void _setSliceIsoEdges( LocalDepth depth , int slice , int z , std::vector< _SlabValues< Vertex > >& slabValues , int threads );
858 template< class Vertex >
859 void _setXSliceIsoEdges( LocalDepth depth , int slice , std::vector< _SlabValues< Vertex > >& slabValues , int threads );
860 template< class Vertex >
861 void _copyFinerSliceIsoEdgeKeys( LocalDepth depth , int slice , std::vector< _SlabValues< Vertex > >& sValues , int threads );
862 template< class Vertex >
863 void _copyFinerSliceIsoEdgeKeys( LocalDepth depth , int slice , int z , std::vector< _SlabValues< Vertex > >& sValues , int threads );
864 template< class Vertex >
865 void _copyFinerXSliceIsoEdgeKeys( LocalDepth depth , int slab , std::vector< _SlabValues< Vertex > >& sValues , int threads );
866
867 template< class Vertex >
868 void _setIsoSurface( LocalDepth depth , int offset , const _SliceValues< Vertex >& bValues , const _SliceValues< Vertex >& fValues , const _XSliceValues< Vertex >& xValues , CoredMeshData< Vertex >& mesh , bool polygonMesh , bool addBarycenter , int& vOffset , int threads );
869
870 template< class Vertex >
871 static int _addIsoPolygons( CoredMeshData< Vertex >& mesh , std::vector< std::pair< int , Vertex > >& polygon , bool polygonMesh , bool addBarycenter , int& vOffset );
872
873 template< int WeightDegree , int ColorDegree , BoundaryType BType , class Vertex >
874 bool _getIsoVertex( const BSplineData< ColorDegree , BType >* colorBSData , const DensityEstimator< WeightDegree >* densityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > , Real > , ColorDegree >* colorData , Real isoValue , ConstPointSupportKey< WeightDegree >& weightKey , ConstPointSupportKey< ColorDegree >& colorKey , const TreeOctNode* node , int edgeIndex , int z , const _SliceValues< Vertex >& sValues , Vertex& vertex );
875 template< int WeightDegree , int ColorDegree , BoundaryType BType , class Vertex >
876 bool _getIsoVertex( const BSplineData< ColorDegree , BType >* colorBSData , const DensityEstimator< WeightDegree >* densityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > , Real > , ColorDegree >* colorData , Real isoValue , ConstPointSupportKey< WeightDegree >& weightKey , ConstPointSupportKey< ColorDegree >& colorKey , const TreeOctNode* node , int cornerIndex , const _SliceValues< Vertex >& bValues , const _SliceValues< Vertex >& fValues , Vertex& vertex );
877
878 void _init( TreeOctNode* node , LocalDepth maxDepth , bool (*Refine)( LocalDepth d , LocalOffset off ) );
879
880 double _maxMemoryUsage , _localMemoryUsage;
881 public:
882 int threads;
883 double maxMemoryUsage( void ) const { return _maxMemoryUsage; }
884 double localMemoryUsage( void ) const { return _localMemoryUsage; }
885 void resetLocalMemoryUsage( void ){ _localMemoryUsage = 0; }
886 double memoryUsage( void );
887
888 Octree( void );
889
890 void init( LocalDepth maxDepth , bool (*Refine)( LocalDepth d , LocalOffset off ) );
891 template< class Data >
892 int init( OrientedPointStream< Real >& pointStream , LocalDepth maxDepth , bool useConfidence , std::vector< PointSample >& samples , std::vector< ProjectiveData< Data , Real > >* sampleData );
893 template< int DensityDegree >
894 typename Octree::template DensityEstimator< DensityDegree >* setDensityEstimator( const std::vector< PointSample >& samples , LocalDepth splatDepth , Real samplesPerNode );
895 template< int NormalDegree , int DensityDegree >
896 SparseNodeData< Point3D< Real > , NormalDegree > setNormalField( const std::vector< PointSample >& samples , const DensityEstimator< DensityDegree >& density , Real& pointWeightSum , bool forceNeumann );
897 template< int DataDegree , bool CreateNodes , int DensityDegree , class Data >
898 SparseNodeData< ProjectiveData< Data , Real > , DataDegree > setDataField( const std::vector< PointSample >& samples , std::vector< ProjectiveData< Data , Real > >& sampleData , const DensityEstimator< DensityDegree >* density );
899 template< int MaxDegree , int FEMDegree , BoundaryType FEMBType , class HasDataFunctor > void inalizeForBroodedMultigrid( LocalDepth fullDepth , const HasDataFunctor& F , std::vector< int >* map=NULL );
900
901 // Generate an empty set of constraints
902 template< int FEMDegree > DenseNodeData< Real , FEMDegree > initDenseNodeData( void );
903
904 // Add finite-elements constraints (derived from a sparse scalar field)
905 template< int FEMDegree , BoundaryType FEMBType , int SFDegree , BoundaryType SFBType , class FEMSFConstraintFunctor > void addFEMConstraints( const FEMSFConstraintFunctor& F , const SparseNodeData< Real , SFDegree >& sfCoefficients , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxDepth )
906 { return _addFEMConstraints< FEMDegree , FEMBType , SFDegree , SFBType , FEMSFConstraintFunctor , const SparseNodeData< Real , SFDegree > , Real , double >( F , sfCoefficients , constraints , maxDepth ); }
907 // Add finite-elements constraints (derived from a dense scalar field)
908 template< int FEMDegree , BoundaryType FEMBType , int SFDegree , BoundaryType SFBType , class FEMSFConstraintFunctor > void addFEMConstraints( const FEMSFConstraintFunctor& F , const DenseNodeData< Real , SFDegree >& sfCoefficients , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxDepth )
909 { return _addFEMConstraints< FEMDegree , FEMBType , SFDegree , SFBType , FEMSFConstraintFunctor , const DenseNodeData< Real , SFDegree > , Real , double >( F , sfCoefficients , constraints , maxDepth ); }
910 // Add finite-elements constraints (derived from a sparse vector field)
911 template< int FEMDegree , BoundaryType FEMBType , int VFDegree , BoundaryType VFBType , class FEMVFConstraintFunctor > void addFEMConstraints( const FEMVFConstraintFunctor& F , const SparseNodeData< Point3D< Real > , VFDegree >& vfCoefficients , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxDepth )
912 { return _addFEMConstraints< FEMDegree , FEMBType , VFDegree , VFBType , FEMVFConstraintFunctor , const SparseNodeData< Point3D< Real > , VFDegree > , Point3D< Real > , Point3D< double > >( F , vfCoefficients , constraints , maxDepth ); }
913 // Add finite-elements constraints (derived from a dense vector field)
914 template< int FEMDegree , BoundaryType FEMBType , int VFDegree , BoundaryType VFBType , class FEMVFConstraintFunctor > void addFEMConstraints( const FEMVFConstraintFunctor& F , const DenseNodeData< Point3D< Real > , VFDegree >& vfCoefficients , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxDepth )
915 { return _addFEMConstraints< FEMDegree , FEMBType , VFDegree , VFBType , FEMVFConstraintFunctor , const DenseNodeData< Point3D< Real > , VFDegree > , Point3D< Real > , Point3D< double > >( F , vfCoefficients , constraints , maxDepth ); }
916 // Add interpolation constraints
917 template< int FEMDegree , BoundaryType FEMBType , bool HasGradients > void addInterpolationConstraints( const InterpolationInfo< HasGradients >& interpolationInfo , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxDepth );
918
919 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 , class DotFunctor > double dot( const DotFunctor& F , const SparseNodeData< Real , Degree1 >& coefficients1 , const SparseNodeData< Real , Degree2 >& coefficients2 ) const
920 { return _dot< Degree1 , BType1 , Degree2 , BType2 , DotFunctor , false >( F , (const InterpolationInfo< false >*)NULL , coefficients1 , coefficients2 ); }
921 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 , class DotFunctor > double dot( const DotFunctor& F , const SparseNodeData< Real , Degree1 >& coefficients1 , const DenseNodeData< Real , Degree2 >& coefficients2 ) const
922 { return _dot< Degree1 , BType1 , Degree2 , BType2 , DotFunctor , false >( F , (const InterpolationInfo< false >*)NULL , coefficients1 , coefficients2 ); }
923 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 , class DotFunctor > double dot( const DotFunctor& F , const DenseNodeData< Real , Degree1 >& coefficients1 , const SparseNodeData< Real , Degree2 >& coefficients2 ) const
924 { return _dot< Degree1 , BType1 , Degree2 , BType2 , DotFunctor , false >( F , (const InterpolationInfo< false >*)NULL , coefficients1 , coefficients2 ); }
925 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 , class DotFunctor > double dot( const DotFunctor& F , const DenseNodeData< Real , Degree1 >& coefficients1 , const DenseNodeData< Real , Degree2 >& coefficients2 ) const
926 { return _dot< Degree1 , BType1 , Degree2 , BType2 , DotFunctor , false >( F , (const InterpolationInfo< false >*)NULL , coefficients1 , coefficients2 ); }
927
928 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 , class DotFunctor , bool HasGradients > double dot( const DotFunctor& F , const InterpolationInfo< HasGradients >* iInfo , const SparseNodeData< Real , Degree1 >& coefficients1 , const SparseNodeData< Real , Degree2 >& coefficients2 ) const
929 { return _dot< Degree1 , BType1 , Degree2 , BType2 , DotFunctor , HasGradients >( F , iInfo , coefficients1 , coefficients2 ); }
930 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 , class DotFunctor , bool HasGradients > double dot( const DotFunctor& F , const InterpolationInfo< HasGradients >* iInfo , const SparseNodeData< Real , Degree1 >& coefficients1 , const DenseNodeData< Real , Degree2 >& coefficients2 ) const
931 { return _dot< Degree1 , BType1 , Degree2 , BType2 , DotFunctor , HasGradients >( F , iInfo , coefficients1 , coefficients2 ); }
932 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 , class DotFunctor , bool HasGradients > double dot( const DotFunctor& F , const InterpolationInfo< HasGradients >* iInfo , const DenseNodeData< Real , Degree1 >& coefficients1 , const SparseNodeData< Real , Degree2 >& coefficients2 ) const
933 { return _dot< Degree1 , BType1 , Degree2 , BType2 , DotFunctor , HasGradients >( F , iInfo , coefficients1 , coefficients2 ); }
934 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 , class DotFunctor , bool HasGradients > double dot( const DotFunctor& F , const InterpolationInfo< HasGradients >* iInfo , const DenseNodeData< Real , Degree1 >& coefficients1 , const DenseNodeData< Real , Degree2 >& coefficients2 ) const
935 { return _dot< Degree1 , BType1 , Degree2 , BType2 , DotFunctor , HasGradients >( F , iInfo , coefficients1 , coefficients2 ); }
936
937 template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients >
938 void setSystemMatrix( const FEMSystemFunctor& F , const InterpolationInfo< HasGradients >* interpolationInfo , LocalDepth depth , SparseMatrix< Real >& matrix ) const;
939
940 // Solve the linear system
941 struct SolverInfo
942 {
943 // How to solve
944 LocalDepth cgDepth;
945 int iters;
946 double cgAccuracy , lowResIterMultiplier;
947 // What to output
948 bool verbose , showResidual;
949
950 SolverInfo( void ) : cgDepth(0) , iters(1), cgAccuracy(0) , lowResIterMultiplier(0) , verbose(false) , showResidual(false) { ; }
951 };
952 template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients >
953 DenseNodeData< Real , FEMDegree > solveSystem( const FEMSystemFunctor& F , InterpolationInfo< HasGradients >* iData , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxSolveDepth , const SolverInfo& solverInfo );
954
955 template< int FEMDegree , BoundaryType BType , int WeightDegree , int ColorDegree , class Vertex >
956 void getMCIsoSurface( const DensityEstimator< WeightDegree >* densityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > , Real > , ColorDegree >* colorData , const DenseNodeData< Real , FEMDegree >& solution , Real isoValue , CoredMeshData< Vertex >& mesh , bool nonLinearFit=true , bool addBarycenter=false , bool polygonMesh=false );
957
958
959 const TreeOctNode& tree( void ) const{ return *_tree; }
960 size_t leaves( void ) const { return _tree->leaves(); }
961 size_t nodes( void ) const { int count = 0 ; for( const TreeOctNode* n=_tree->nextNode() ; n ; n=_tree->nextNode( n ) ) if( IsActiveNode( n ) ) count++ ; return count; }
962 size_t ghostNodes( void ) const { int count = 0 ; for( const TreeOctNode* n=_tree->nextNode() ; n ; n=_tree->nextNode( n ) ) if( !IsActiveNode( n ) ) count++ ; return count; }
963 inline size_t validSpaceNodes( void ) const { int count = 0 ; for( const TreeOctNode* n=_tree->nextNode() ; n ; n=_tree->nextNode( n ) ) if( isValidSpaceNode( n ) ) count++ ; return count; }
964 inline size_t validSpaceNodes( LocalDepth d ) const { int count = 0 ; for( const TreeOctNode* n=_tree->nextNode() ; n ; n=_tree->nextNode( n ) ) if( _localDepth(n)==d && isValidSpaceNode( n ) ) count++ ; return count; }
965 template< int Degree , BoundaryType BType > size_t validFEMNodes( void ) const { int count = 0 ; for( const TreeOctNode* n=_tree->nextNode() ; n ; n=_tree->nextNode( n ) ) if( isValidFEMNode< Degree , BType >( n ) ) count++ ; return count; }
966 template< int Degree , BoundaryType BType > size_t validFEMNodes( LocalDepth d ) const { int count = 0 ; for( const TreeOctNode* n=_tree->nextNode() ; n ; n=_tree->nextNode( n ) ) if( _localDepth(n)==d && isValidFEMNode< Degree , BType >( n ) ) count++ ; return count; }
967 LocalDepth depth( void ) const { return _localMaxDepth( _tree ); }
968 void resetNodeIndices( void ){ _NodeCount = 0 ; for( TreeOctNode* node=_tree->nextNode() ; node ; node=_tree->nextNode( node ) ) _NodeInitializer( *node ) , node->nodeData.flags=0; }
969
970 protected:
971 template< class D > static bool _IsZero( const D& d );
972 template< class D > static Real _Dot( const D& d1 , const D& d2 );
973 template< int FEMDegree , BoundaryType FEMBType , int CDegree , BoundaryType CBType , class FEMConstraintFunctor , class Coefficients , class D , class _D >
974 void _addFEMConstraints( const FEMConstraintFunctor& F , const Coefficients& coefficients , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxDepth );
975 template< int FEMDegree1 , BoundaryType FEMBType1 , int FEMDegree2 , BoundaryType FEMBType2 , class DotFunctor , bool HasGradients , class Coefficients1 , class Coefficients2 >
976 double _dot( const DotFunctor& F , const InterpolationInfo< HasGradients >* iInfo , const Coefficients1& coefficients1 , const Coefficients2& coefficients2 ) const;
977 };
978 template< class Real > int Octree< Real >::_NodeCount = 0;
979
980
981 template< class Real > void Reset( void ){ Octree< Real >::ResetNodeCount(); }
982
983
984 #include "MultiGridOctreeData.inl"
985 #include "MultiGridOctreeData.SortedTreeNodes.inl"
986 #include "MultiGridOctreeData.WeightedSamples.inl"
987 #include "MultiGridOctreeData.System.inl"
988 #include "MultiGridOctreeData.IsoSurface.inl"
989 #include "MultiGridOctreeData.Evaluation.inl"
990 #endif // MULTI_GRID_OCTREE_DATA_INCLUDED
991