diff --git a/src/meshlabplugins/filter_screened_poisson/Src/Array.inl b/src/meshlabplugins/filter_screened_poisson/Src/Array.inl index 511b9a082..6247f5282 100755 --- a/src/meshlabplugins/filter_screened_poisson/Src/Array.inl +++ b/src/meshlabplugins/filter_screened_poisson/Src/Array.inl @@ -188,7 +188,7 @@ public: // [WARNING] Chaning szC and szD to size_t causes some really strange behavior. long long szC = sizeof( C ); long long szD = sizeof( D ); - data = (C*)&a[0]; + data = (C*)a.data; min = ( a.minimum() * szD ) / szC; max = ( a.maximum() * szD ) / szC; if( min*szC!=a.minimum()*szD || max*szC!=a.maximum()*szD ) @@ -306,6 +306,7 @@ public: return (*this); } inline Array& operator ++ ( void ) { return (*this) += 1; } + inline Array operator++( int ){ Array< C > temp = (*this) ; (*this) +=1 ; return temp; } Array operator - ( int idx ) const { return (*this) + (-idx); } Array operator - ( long long idx ) const { return (*this) + (-idx); } Array operator - ( unsigned int idx ) const { return (*this) + (-idx); } @@ -315,6 +316,7 @@ public: Array& operator -= ( unsigned int idx ) { return (*this) += (-idx); } Array& operator -= ( unsigned long long idx ) { return (*this) += (-idx); } Array& operator -- ( void ) { return (*this) -= 1; } + inline Array operator--( int ){ Array< C > temp = (*this) ; (*this) -=1 ; return temp; } long long operator - ( const Array& a ) const { return ( long long )( data - a.data ); } void Free( void ) @@ -501,6 +503,7 @@ public: return (*this); } inline ConstArray& operator ++ ( void ) { return (*this) += 1; } + inline ConstArray operator++( int ){ ConstArray< C > temp = (*this) ; (*this) +=1 ; return temp; } ConstArray operator - ( int idx ) const { return (*this) + (-idx); } ConstArray operator - ( long long idx ) const { return (*this) + (-idx); } ConstArray operator - ( unsigned int idx ) const { return (*this) + (-idx); } @@ -510,6 +513,7 @@ public: ConstArray& operator -= ( unsigned int idx ) { return (*this) += (-idx); } ConstArray& operator -= ( unsigned long long idx ) { return (*this) += (-idx); } ConstArray& operator -- ( void ) { return (*this) -= 1; } + inline ConstArray operator--( int ){ ConstArray< C > temp = (*this) ; (*this) -=1 ; return temp; } long long operator - ( const ConstArray& a ) const { return ( long long )( data - a.data ); } long long operator - ( const Array< C >& a ) const { return ( long long )( data - a.pointer() ); } diff --git a/src/meshlabplugins/filter_screened_poisson/Src/BinaryNode.h b/src/meshlabplugins/filter_screened_poisson/Src/BinaryNode.h index 96f306172..7c5a17153 100755 --- a/src/meshlabplugins/filter_screened_poisson/Src/BinaryNode.h +++ b/src/meshlabplugins/filter_screened_poisson/Src/BinaryNode.h @@ -29,9 +29,6 @@ DAMAGE. #ifndef BINARY_NODE_INCLUDED #define BINARY_NODE_INCLUDED -#define MSVC_2010_FIX 1 - - class BinaryNode { public: @@ -40,39 +37,34 @@ public: static inline int CumulativeCenterCount( int maxDepth ) { return (1<<(maxDepth+1))-1; } static inline int CumulativeCornerCount( int maxDepth ) { return (1<<(maxDepth+1))+maxDepth; } static inline int CenterIndex( int depth , int offSet ) { return (1< static inline Real CornerIndexPosition(int index,int maxDepth){ return Real(index)/(1< static inline Real Width(int depth){ return Real(1.0/(1< static inline void CenterAndWidth( int depth , int offset , Real& center , Real& width ) - { - width=Real (1.0/(1< static inline Real Width( int depth ){ return Real(1.0/(1< static inline void CenterAndWidth( int depth , int offset , Real& center , Real& width ){ width = Real (1.0/(1< static inline void CornerAndWidth( int depth , int offset , Real& corner , Real& width ){ width = Real(1.0/(1< static inline void CenterAndWidth( int idx , Real& center , Real& width ) - { - int depth , offset; - DepthAndOffset( idx , depth , offset ); - CenterAndWidth( depth , offset , center , width ); - } - static inline void DepthAndOffset( int idx , int& depth , int& offset ) - { - int i=idx+1; -#if MSVC_2010_FIX - depth = 0; -#else // !MSVC_2010_FIX - depth = -1; -#endif // MSVC_2010_FIX - while( i ) - { - i >>= 1; - depth++; - } -#if MSVC_2010_FIX - depth--; -#endif // MSVC_2010_FIX - offset = ( idx+1 ) - (1< static inline void CornerAndWidth( int idx , Real& corner , Real& width ) + { + int depth , offset; + CornerDepthAndOffset( idx , depth , offset ); + CornerAndWidth( depth , offset , corner , width ); + } + static inline void CenterDepthAndOffset( int idx , int& depth , int& offset ) + { + offset = idx , depth = 0; + while( offset>=(1<=( (1< FunctionData::FunctionData(void) diff --git a/src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.Evaluation.inl b/src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.Evaluation.inl new file mode 100644 index 000000000..c965effbe --- /dev/null +++ b/src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.Evaluation.inl @@ -0,0 +1,1151 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. Redistributions in binary form must reproduce +the above copyright notice, this list of conditions and the following disclaimer +in the documentation and/or other materials provided with the distribution. + +Neither the name of the Johns Hopkins University nor the names of its contributors +may be used to endorse or promote products derived from this software without specific +prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +template< class Real > +template< int FEMDegree , BoundaryType BType> +void Octree< Real >::_Evaluator< FEMDegree , BType >::set( LocalDepth depth ) +{ + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + + BSplineEvaluationData< FEMDegree , BType >::SetEvaluator( evaluator , depth ); + if( depth>0 ) BSplineEvaluationData< FEMDegree , BType >::SetChildEvaluator( childEvaluator , depth-1 ); + int center = ( 1<>1; + + // First set the stencils for the current depth + for( int x=-LeftPointSupportRadius ; x<=RightPointSupportRadius ; x++ ) for( int y=-LeftPointSupportRadius ; y<=RightPointSupportRadius ; y++ ) for( int z=-LeftPointSupportRadius ; z<=RightPointSupportRadius ; z++ ) + { + int fIdx[] = { center+x , center+y , center+z }; + + // The cell stencil + { + double vv[3] , dv[3]; + for( int dd=0 ; dd( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] ); + } + + //// The face stencil + for( int f=0 ; f( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] ); + } + + //// The edge stencil + for( int e=0 ; e( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] ); + } + + //// The corner stencil + for( int c=0 ; c( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] ); + } + } + + // Now set the stencils for the parents + for( int child=0 ; child( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] ); + } + + //// The face stencil + for( int f=0 ; f( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] ); + } + + //// The edge stencil + for( int e=0 ; e( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] ); + } + + //// The corner stencil + for( int c=0 ; c( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] ); + } + } + } + if( _bsData ) delete _bsData; + _bsData = new BSplineData< FEMDegree , BType >( depth ); +} +template< class Real > +template< class V , int FEMDegree , BoundaryType BType > +V Octree< Real >::_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 +{ + static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize; + static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart; + + if( IsActiveNode( node->children ) ) fprintf( stderr , "[WARNING] getValue assumes leaf node\n" ); + V value(0); + + while( GetGhostFlag( node ) ) + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + + for( int i=0 ; i _s ; Real _w; + _startAndWidth( _n , _s , _w ); + int _fIdx[3]; + functionIndex< FEMDegree , BType >( _n , _fIdx ); + for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) ); + value += + solution[ _n->nodeData.nodeIndex ] * + (Real) + ( + evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * + evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * + evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) + ); + } + } + node = node->parent; + } + + LocalDepth d = _localDepth( node ); + + for( int dd=0 ; dd<3 ; dd++ ) + if ( p[dd]==0 ) p[dd] = (Real)(0.+1e-6); + else if( p[dd]==1 ) p[dd] = (Real)(1.-1e-6); + + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + + for( int i=0 ; i _s ; Real _w; + _startAndWidth( _n , _s , _w ); + int _fIdx[3]; + functionIndex< FEMDegree , BType >( _n , _fIdx ); + for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) ); + value += + solution[ _n->nodeData.nodeIndex ] * + (Real) + ( + evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * + evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * + evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) + ); + } + } + if( d>0 ) + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + for( int i=0 ; i _s ; Real _w; + _startAndWidth( _n , _s , _w ); + int _fIdx[3]; + functionIndex< FEMDegree , BType >( _n , _fIdx ); + for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) ); + value += + coarseSolution[ _n->nodeData.nodeIndex ] * + (Real) + ( + evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * + evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * + evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) + ); + } + } + } + } + return value; +} +template< class Real > +template< int FEMDegree , BoundaryType BType > +std::pair< Real , Point3D< Real > > Octree< 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 +{ + static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize; + static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart; + + if( IsActiveNode( node->children ) ) fprintf( stderr , "[WARNING] _getValueAndGradient assumes leaf node\n" ); + Real value(0); + Point3D< Real > gradient; + + while( GetGhostFlag( node ) ) + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + + for( int i=0 ; i _s; Real _w; + _startAndWidth( _n , _s , _w ); + int _fIdx[3]; + functionIndex< FEMDegree , BType >( _n , _fIdx ); + for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) ); + value += + solution[ _n->nodeData.nodeIndex ] * + (Real) + ( + evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) + ); + gradient += + Point3D< Real > + ( + evaluator._bsData->dBaseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) , + evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->dBaseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) , + evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->dBaseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) + ) * solution[ _n->nodeData.nodeIndex ]; + } + } + node = node->parent; + } + + + LocalDepth d = _localDepth( node ); + + for( int dd=0 ; dd<3 ; dd++ ) + if ( p[dd]==0 ) p[dd] = (Real)(0.+1e-6); + else if( p[dd]==1 ) p[dd] = (Real)(1.-1e-6); + + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + + for( int i=0 ; i _s ; Real _w; + _startAndWidth( _n , _s , _w ); + int _fIdx[3]; + functionIndex< FEMDegree , BType >( _n , _fIdx ); + for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) ); + value += + solution[ _n->nodeData.nodeIndex ] * + (Real) + ( + evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) + ); + gradient += + Point3D< Real > + ( + evaluator._bsData->dBaseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) , + evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->dBaseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) , + evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->dBaseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) + ) * solution[ _n->nodeData.nodeIndex ]; + } + } + if( d>0 ) + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + for( int i=0 ; i _s ; Real _w; + _startAndWidth( _n , _s , _w ); + int _fIdx[3]; + functionIndex< FEMDegree , BType >( _n , _fIdx ); + for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) ); + value += + coarseSolution[ _n->nodeData.nodeIndex ] * + (Real) + ( + evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) + ); + gradient += + Point3D< Real > + ( + evaluator._bsData->dBaseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) , + evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->dBaseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) , + evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->dBaseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) + ) * coarseSolution[ _n->nodeData.nodeIndex ]; + } + } + } + } + return std::pair< Real , Point3D< Real > >( value , gradient ); +} +template< class Real > +template< class V , int FEMDegree , BoundaryType BType > +V Octree< Real >::_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 +{ + static const int SupportSize = BSplineEvaluationData< FEMDegree , BType >::SupportSize; + static const int LeftPointSupportRadius = BSplineEvaluationData< FEMDegree , BType >::SupportEnd; + static const int RightPointSupportRadius = - BSplineEvaluationData< FEMDegree , BType >::SupportStart; + + if( IsActiveNode( node->children ) ) fprintf( stderr , "[WARNING] getCenterValue assumes leaf node\n" ); + V value(0); + LocalDepth d = _localDepth( node ); + + if( isInterior ) + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + for( int i=0 ; inodeData.nodeIndex ] * Real( evaluator.cellStencil( i , j , k ) ); + } + if( d>0 ) + { + int _corner = int( node - node->parent->children ); + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + for( int i=0 ; inodeData.nodeIndex] * Real( evaluator.cellStencils[_corner]( i , j , k ) ); + } + } + } + else + { + LocalOffset cIdx; + _localDepthAndOffset( node , d , cIdx ); + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + + for( int i=0 ; inodeData.nodeIndex ] * + Real( + evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false ) + ); + } + } + if( d>0 ) + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + for( int i=0 ; inodeData.nodeIndex ] * + Real( + evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false ) + ); + } + } + } + } + return value; +} +template< class Real > +template< int FEMDegree , BoundaryType BType > +std::pair< Real , Point3D< Real > > Octree< 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 +{ + static const int SupportSize = BSplineEvaluationData< FEMDegree , BType >::SupportSize; + static const int LeftPointSupportRadius = BSplineEvaluationData< FEMDegree , BType >::SupportEnd; + static const int RightPointSupportRadius = - BSplineEvaluationData< FEMDegree , BType >::SupportStart; + + if( IsActiveNode( node->children ) ) fprintf( stderr , "[WARNING] getCenterValueAndGradient assumes leaf node\n" ); + Real value(0); + Point3D< Real > gradient; + LocalDepth d = _localDepth( node ); + + if( isInterior ) + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + for( int i=0 ; inodeData.nodeIndex ]; + gradient += Point3D< Real >( evaluator.dCellStencil( i , j , k ) ) * solution[ n->nodeData.nodeIndex ]; + } + } + if( d>0 ) + { + int _corner = int( node - node->parent->children ); + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + for( int i=0 ; inodeData.nodeIndex]; + gradient += Point3D< Real >( evaluator.dCellStencils[_corner]( i , j , k ) ) * coarseSolution[n->nodeData.nodeIndex]; + } + } + } + } + else + { + LocalOffset cIdx; + _localDepthAndOffset( node , d , cIdx ); + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + + for( int i=0 ; inodeData.nodeIndex ]; + gradient += + Point3D< Real > + ( + evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , true ) * evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false ) * evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false ) , + evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false ) * evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , true ) * evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false ) , + evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false ) * evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false ) * evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , true ) + ) * solution[ n->nodeData.nodeIndex ]; + } + } + if( d>0 ) + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + for( int i=0 ; inodeData.nodeIndex ]; + gradient += + Point3D< Real > + ( + evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , true ) * evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false ) * evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false ) , + evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false ) * evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , true ) * evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false ) , + evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false ) * evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false ) * evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , true ) + ) * coarseSolution[ n->nodeData.nodeIndex ]; + } + } + } + } + return std::pair< Real , Point3D< Real > >( value , gradient ); +} +template< class Real > +template< class V , int FEMDegree , BoundaryType BType > +V Octree< Real >::_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 +{ + static const int SupportSize = BSplineEvaluationData< FEMDegree , BType >::SupportSize; + static const int LeftPointSupportRadius = BSplineEvaluationData< FEMDegree , BType >::SupportEnd; + static const int RightPointSupportRadius = -BSplineEvaluationData< FEMDegree , BType >::SupportStart; + V value(0); + LocalDepth d ; LocalOffset cIdx; + _localDepthAndOffset( node , d , cIdx ); + int startX = 0 , endX = SupportSize , startY = 0 , endY = SupportSize , startZ = 0 , endZ = SupportSize; + int orientation , i1 , i2; + Cube::FactorEdgeIndex( edge , orientation , i1 , i2 ); + switch( orientation ) + { + case 0: + cIdx[1] += i1 , cIdx[2] += i2; + if( i1 ) startY++ ; else endY--; + if( i2 ) startZ++ ; else endZ--; + break; + case 1: + cIdx[0] += i1 , cIdx[2] += i2; + if( i1 ) startX++ ; else endX--; + if( i2 ) startZ++ ; else endZ--; + break; + case 2: + cIdx[0] += i1 , cIdx[1] += i2; + if( i1 ) startX++ ; else endX--; + if( i2 ) startY++ ; else endY--; + break; + } + + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , d ); + for( int x=startX ; xnodeData.nodeIndex ] * evaluator.edgeStencil[edge]( x , y , z ); + else + { + LocalDepth _d ; LocalOffset fIdx; + _localDepthAndOffset( _node , _d , fIdx ); + switch( orientation ) + { + case 0: + value += + solution[ _node->nodeData.nodeIndex ] * + Real( + evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false ) + ); + break; + case 1: + value += + solution[ _node->nodeData.nodeIndex ] * + Real( + evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false ) + ); + break; + case 2: + value += + solution[ _node->nodeData.nodeIndex ] * + Real( + evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false ) + ); + break; + } + } + } + } + } + if( d>0 ) + { + int _corner = int( node - node->parent->children ); + int _cx , _cy , _cz; + Cube::FactorCornerIndex( _corner , _cx , _cy , _cz ); + // If the corner/child indices don't match, then the sample position is in the interior of the + // coarser cell and so the full support resolution should be used. + switch( orientation ) + { + case 0: + if( _cy!=i1 ) startY = 0 , endY = SupportSize; + if( _cz!=i2 ) startZ = 0 , endZ = SupportSize; + break; + case 1: + if( _cx!=i1 ) startX = 0 , endX = SupportSize; + if( _cz!=i2 ) startZ = 0 , endZ = SupportSize; + break; + case 2: + if( _cx!=i1 ) startX = 0 , endX = SupportSize; + if( _cy!=i2 ) startY = 0 , endY = SupportSize; + break; + } + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + for( int x=startX ; xnodeData.nodeIndex ] * evaluator.edgeStencils[_corner][edge]( x , y , z ); + else + { + LocalDepth _d ; LocalOffset fIdx; + _localDepthAndOffset( _node , _d , fIdx ); + switch( orientation ) + { + case 0: + value += + coarseSolution[ _node->nodeData.nodeIndex ] * + Real( + evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false ) + ); + break; + case 1: + value += + coarseSolution[ _node->nodeData.nodeIndex ] * + Real( + evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false ) + ); + break; + case 2: + value += + coarseSolution[ _node->nodeData.nodeIndex ] * + Real( + evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false ) + ); + break; + } + } + } + } + } + return Real( value ); +} +template< class Real > +template< int FEMDegree , BoundaryType BType > +std::pair< Real , Point3D< Real > > Octree< 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 +{ + static const int SupportSize = BSplineEvaluationData< FEMDegree , BType >::SupportSize; + static const int LeftPointSupportRadius = BSplineEvaluationData< FEMDegree , BType >::SupportEnd; + static const int RightPointSupportRadius = -BSplineEvaluationData< FEMDegree , BType >::SupportStart; + double value = 0; + Point3D< double > gradient; + LocalDepth d ; LocalOffset cIdx; + _localDepthAndOffset( node , d , cIdx ); + + int startX = 0 , endX = SupportSize , startY = 0 , endY = SupportSize , startZ = 0 , endZ = SupportSize; + int orientation , i1 , i2; + Cube::FactorEdgeIndex( edge , orientation , i1 , i2 ); + switch( orientation ) + { + case 0: + cIdx[1] += i1 , cIdx[2] += i2; + if( i1 ) startY++ ; else endY--; + if( i2 ) startZ++ ; else endZ--; + break; + case 1: + cIdx[0] += i1 , cIdx[2] += i2; + if( i1 ) startX++ ; else endX--; + if( i2 ) startZ++ ; else endZ--; + break; + case 2: + cIdx[0] += i1 , cIdx[1] += i2; + if( i1 ) startX++ ; else endX--; + if( i2 ) startY++ ; else endY--; + break; + } + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + for( int x=startX ; xnodeData.nodeIndex ]; + gradient += evaluator.dEdgeStencil[edge]( x , y , z ) * solution[ _node->nodeData.nodeIndex ]; + } + else + { + LocalDepth _d ; LocalOffset fIdx; + _localDepthAndOffset( _node , _d , fIdx ); + + double vv[3] , dv[3]; + switch( orientation ) + { + case 0: + vv[0] = evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false ); + vv[1] = evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false ); + vv[2] = evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false ); + dv[0] = evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , true ); + dv[1] = evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , true ); + dv[2] = evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , true ); + break; + case 1: + vv[0] = evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false ); + vv[1] = evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false ); + vv[2] = evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false ); + dv[0] = evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , true ); + dv[1] = evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , true ); + dv[2] = evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , true ); + break; + case 2: + vv[0] = evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false ); + vv[1] = evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false ); + vv[2] = evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false ); + dv[0] = evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , true ); + dv[1] = evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , true ); + dv[2] = evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , true ); + break; + } + value += solution[ _node->nodeData.nodeIndex ] * vv[0] * vv[1] * vv[2]; + gradient += Point3D< double >( dv[0]*vv[1]*vv[2] , vv[0]*dv[1]*vv[2] , vv[0]*vv[1]*dv[2] ) * solution[ _node->nodeData.nodeIndex ]; + } + } + } + } + if( d>0 ) + { + int _corner = int( node - node->parent->children ); + int _cx , _cy , _cz; + Cube::FactorCornerIndex( _corner , _cx , _cy , _cz ); + // If the corner/child indices don't match, then the sample position is in the interior of the + // coarser cell and so the full support resolution should be used. + switch( orientation ) + { + case 0: + if( _cy!=i1 ) startY = 0 , endY = SupportSize; + if( _cz!=i2 ) startZ = 0 , endZ = SupportSize; + break; + case 1: + if( _cx!=i1 ) startX = 0 , endX = SupportSize; + if( _cz!=i2 ) startZ = 0 , endZ = SupportSize; + break; + case 2: + if( _cx!=i1 ) startX = 0 , endX = SupportSize; + if( _cy!=i2 ) startY = 0 , endY = SupportSize; + break; + } + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + for( int x=startX ; xnodeData.nodeIndex ]; + gradient += evaluator.dEdgeStencils[_corner][edge]( x , y , z ) * coarseSolution[ _node->nodeData.nodeIndex ]; + } + else + { + LocalDepth _d ; LocalOffset fIdx; + _localDepthAndOffset( _node , _d , fIdx ); + double vv[3] , dv[3]; + switch( orientation ) + { + case 0: + vv[0] = evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false ); + vv[1] = evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false ); + vv[2] = evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false ); + dv[0] = evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , true ); + dv[1] = evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , true ); + dv[2] = evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , true ); + break; + case 1: + vv[0] = evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false ); + vv[1] = evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false ); + vv[2] = evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false ); + dv[0] = evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , true ); + dv[1] = evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , true ); + dv[2] = evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , true ); + break; + case 2: + vv[0] = evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false ); + vv[1] = evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false ); + vv[2] = evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false ); + dv[0] = evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , true ); + dv[1] = evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , true ); + dv[2] = evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , true ); + break; + } + value += coarseSolution[ _node->nodeData.nodeIndex ] * vv[0] * vv[1] * vv[2]; + gradient += Point3D< double >( dv[0]*vv[1]*vv[2] , vv[0]*dv[1]*vv[2] , vv[0]*vv[1]*dv[2] ) * coarseSolution[ _node->nodeData.nodeIndex ]; + } + } + } + } + return std::pair< Real , Point3D< Real > >( Real( value ) , Point3D< Real >( gradient ) ); +} + +template< class Real > +template< class V , int FEMDegree , BoundaryType BType > +V Octree< Real >::_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 +{ + static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart; + + V value(0); + LocalDepth d ; LocalOffset cIdx; + _localDepthAndOffset( node , d , cIdx ); + + int cx , cy , cz; + int startX = 0 , endX = SupportSize , startY = 0 , endY = SupportSize , startZ = 0 , endZ = SupportSize; + Cube::FactorCornerIndex( corner , cx , cy , cz ); + cIdx[0] += cx , cIdx[1] += cy , cIdx[2] += cz; + { + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + if( cx==0 ) endX--; + else startX++; + if( cy==0 ) endY--; + else startY++; + if( cz==0 ) endZ--; + else startZ++; + if( isInterior ) + for( int x=startX ; xnodeData.nodeIndex ] * Real( evaluator.cornerStencil[corner]( x , y , z ) ); + } + else + for( int x=startX ; xnodeData.nodeIndex ] * + Real( + evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false ) + ); + } + } + } + if( d>0 ) + { + int _corner = int( node - node->parent->children ); + int _cx , _cy , _cz; + Cube::FactorCornerIndex( _corner , _cx , _cy , _cz ); + // If the corner/child indices don't match, then the sample position is in the interior of the + // coarser cell and so the full support resolution should be used. + if( cx!=_cx ) startX = 0 , endX = SupportSize; + if( cy!=_cy ) startY = 0 , endY = SupportSize; + if( cz!=_cz ) startZ = 0 , endZ = SupportSize; + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + if( isInterior ) + for( int x=startX ; xnodeData.nodeIndex ] * Real( evaluator.cornerStencils[_corner][corner]( x , y , z ) ); + } + else + for( int x=startX ; xnodeData.nodeIndex ] * + Real( + evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false ) * + evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false ) * + evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false ) + ); + } + } + } + return Real( value ); +} +template< class Real > +template< int FEMDegree , BoundaryType BType > +std::pair< Real , Point3D< Real > > Octree< 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 +{ + static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart; + + double value = 0; + Point3D< double > gradient; + LocalDepth d ; LocalOffset cIdx; + _localDepthAndOffset( node , d , cIdx ); + + int cx , cy , cz; + int startX = 0 , endX = SupportSize , startY = 0 , endY = SupportSize , startZ = 0 , endZ = SupportSize; + Cube::FactorCornerIndex( corner , cx , cy , cz ); + cIdx[0] += cx , cIdx[1] += cy , cIdx[2] += cz; + { + if( cx==0 ) endX--; + else startX++; + if( cy==0 ) endY--; + else startY++; + if( cz==0 ) endZ--; + else startZ++; + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node ); + if( isInterior ) + for( int x=startX ; xnodeData.nodeIndex ] * evaluator.cornerStencil[corner]( x , y , z ) , gradient += evaluator.dCornerStencil[corner]( x , y , z ) * solution[ _node->nodeData.nodeIndex ]; + } + else + for( int x=startX ; xnodeData.nodeIndex ] * v[0] * v[1] * v[2]; + gradient += Point3D< double >( dv[0]*v[1]*v[2] , v[0]*dv[1]*v[2] , v[0]*v[1]*dv[2] ) * solution[ _node->nodeData.nodeIndex ]; + } + } + } + if( d>0 ) + { + int _corner = int( node - node->parent->children ); + int _cx , _cy , _cz; + Cube::FactorCornerIndex( _corner , _cx , _cy , _cz ); + if( cx!=_cx ) startX = 0 , endX = SupportSize; + if( cy!=_cy ) startY = 0 , endY = SupportSize; + if( cz!=_cz ) startZ = 0 , endZ = SupportSize; + const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent ); + if( isInterior ) + for( int x=startX ; xnodeData.nodeIndex ] * evaluator.cornerStencils[_corner][corner]( x , y , z ) , gradient += evaluator.dCornerStencils[_corner][corner]( x , y , z ) * coarseSolution[ _node->nodeData.nodeIndex ]; + } + else + for( int x=startX ; xnodeData.nodeIndex ] * v[0] * v[1] * v[2]; + gradient += Point3D< double >( dv[0]*v[1]*v[2] , v[0]*dv[1]*v[2] , v[0]*v[1]*dv[2] ) * coarseSolution[ _node->nodeData.nodeIndex ]; + } + } + } + return std::pair< Real , Point3D< Real > >( Real( value ) , Point3D< Real >( gradient ) ); +} +template< class Real > +template< int Degree , BoundaryType BType > +Octree< Real >::MultiThreadedEvaluator< Degree , BType >::MultiThreadedEvaluator( const Octree< Real >* tree , const DenseNodeData< Real , Degree >& coefficients , int threads ) : _coefficients( coefficients ) , _tree( tree ) +{ + _threads = std::max< int >( 1 , threads ); + _neighborKeys.resize( _threads ); + _coarseCoefficients = _tree->template coarseCoefficients< Real , Degree , BType >( _coefficients ); + _evaluator.set( _tree->_maxDepth ); + for( int t=0 ; t<_threads ; t++ ) _neighborKeys[t].set( tree->_localToGlobal( _tree->_maxDepth ) ); +} +template< class Real > +template< int Degree , BoundaryType BType > +Real Octree< Real >::MultiThreadedEvaluator< Degree , BType >::value( Point3D< Real > p , int thread , const TreeOctNode* node ) +{ + if( !node ) node = _tree->leaf( p ); + ConstPointSupportKey< Degree >& nKey = _neighborKeys[thread]; + nKey.getNeighbors( node ); + return _tree->template _getValue< Real , Degree >( nKey , node , p , _coefficients , _coarseCoefficients , _evaluator ); +} +template< class Real > +template< int Degree , BoundaryType BType > +std::pair< Real , Point3D< Real > > Octree< Real >::MultiThreadedEvaluator< Degree , BType >::valueAndGradient( Point3D< Real > p , int thread , const TreeOctNode* node ) +{ + if( !node ) node = _tree->leaf( p ); + ConstPointSupportKey< Degree >& nKey = _neighborKeys[thread]; + nKey.getNeighbors( node ); + return _tree->template _getValueAndGradient< Degree >( nKey , node , p , _coefficients , _coarseCoefficients , _evaluator ); +} diff --git a/src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.System.inl b/src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.System.inl new file mode 100644 index 000000000..f1cea1956 --- /dev/null +++ b/src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.System.inl @@ -0,0 +1,2274 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. Redistributions in binary form must reproduce +the above copyright notice, this list of conditions and the following disclaimer +in the documentation and/or other materials provided with the distribution. + +Neither the name of the Johns Hopkins University nor the names of its contributors +may be used to endorse or promote products derived from this software without specific +prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +template< class Real , int Degree , bool HasGradients > +struct _ConstraintCalculator_ +{ + static inline Real _CalculateConstraint_( const PointData< Real , HasGradients >& p , const Polynomial< Degree >& px , const Polynomial< Degree >& py , const Polynomial< Degree >& pz , const Polynomial< Degree >& dpx , const Polynomial< Degree >& dpy , const Polynomial< Degree >& dpz , Real valueWeight , Real gradientWeight ); + static inline Real _CalculateConstraint_( const PointData< Real , HasGradients >& p , const Polynomial< Degree >& px , const Polynomial< Degree >& py , const Polynomial< Degree >& pz , const Polynomial< Degree >& dpx , const Polynomial< Degree >& dpy , const Polynomial< Degree >& dpz ); +#if POINT_DATA_RES + static inline void _CalculateCoarser_( int c , PointData< Real , HasGradients >& p , Real value , Point3D< Real > gradient , Real valueWeight , Real gradientWeight ); +#else // !POINT_DATA_RES + static inline void _CalculateCoarser_( PointData< Real , HasGradients >& p , Real value , Point3D< Real > gradient , Real valueWeight , Real gradientWeight ); +#endif // POINT_DATA_RES + +}; +template< class Real , int Degree > +struct _ConstraintCalculator_< Real , Degree , false > +{ + static inline Real _CalculateConstraint_( const PointData< Real , false >& p , const Polynomial< Degree >& px , const Polynomial< Degree >& py , const Polynomial< Degree >& pz , const Polynomial< Degree >& dpx , const Polynomial< Degree >& dpy , const Polynomial< Degree >& dpz , Real valueWeight , Real gradientWeight ) + { +#if POINT_DATA_RES + Real constraint = 0; + for( int c=0 ; c::SAMPLES ; c++ ) if( p[c].weight ) + { + const Point3D< Real > q = p[c].position; + constraint += (Real)( px( q[0] ) * py( q[1] ) * pz( q[2] ) * p[c].weight * p[c].value ); + } + return constraint * valueWeight; +#else // !POINT_DATA_RES + const Point3D< Real > q = p.position; + return (Real)( px( q[0] ) * py( q[1] ) * pz( q[2] ) * p.weight * p.value ) * valueWeight; +#endif // POINT_DATA_RES + } + static inline Real _CalculateConstraint_( const PointData< Real , false >& p , const Polynomial< Degree >& px , const Polynomial< Degree >& py , const Polynomial< Degree >& pz , const Polynomial< Degree >& dpx , const Polynomial< Degree >& dpy , const Polynomial< Degree >& dpz ) + { +#if POINT_DATA_RES + Real constraint = 0; + for( int c=0 ; c::SAMPLES ; c++ ) if( p[c].weight ) + { + const Point3D< Real > q = p[c].position; + constraint += (Real)( px( q[0] ) * py( q[1] ) * pz( q[2] ) * p[c]._value ); + } + return constraint; +#else // !POINT_DATA_RES + const Point3D< Real > q = p.position; + return (Real)( px( q[0] ) * py( q[1] ) * pz( q[2] ) * p._value ); +#endif // POINT_DATA_RES + } +#if POINT_DATA_RES + static inline void _CalculateCoarser_( int c , PointData< Real , false >& p , Real value , Point3D< Real > gradient , Real valueWeight , Real gradientWeight ){ p[c]._value = value * valueWeight * p[c].weight; } +#else // !POINT_DATA_RES + static inline void _CalculateCoarser_( PointData< Real , false >& p , Real value , Point3D< Real > gradient , Real valueWeight , Real gradientWeight ){ p._value = value * valueWeight * p.weight; } +#endif // POINT_DATA_RES +}; +template< class Real , int Degree > +struct _ConstraintCalculator_< Real , Degree , true > +{ + static inline Real _CalculateConstraint_( const PointData< Real , true >& p , const Polynomial< Degree >& px , const Polynomial< Degree >& py , const Polynomial< Degree >& pz , const Polynomial< Degree >& dpx , const Polynomial< Degree >& dpy , const Polynomial< Degree >& dpz , Real valueWeight , Real gradientWeight ) + { +#if POINT_DATA_RES + Real constraint = 0; + for( int c=0 ; c::SAMPLES ; c++ ) if( p[c].weight ) + { + const Point3D< Real > q = p[c].position; + double _px = px( q[0] ) , _py = py( q[1] ) , _pz = pz( q[2] ); + constraint += + ( + (Real)( _px * _py * _pz * p[c].value ) * valueWeight + + Point3D< Real >::Dot( Point3D< Real >( dpx( q[0] ) * _py * _pz , _px * dpy( q[1] ) * _pz , _px * _py * dpz( q[2] ) ) , p[c].gradient ) * gradientWeight + ) * p[c].weight; + } + return constraint; +#else // !POINT_DATA_RES + const Point3D< Real > q = p.position; + double _px = px( q[0] ) , _py = py( q[1] ) , _pz = pz( q[2] ); + return + ( + (Real)( _px * _py * _pz * p.value ) * valueWeight + + Point3D< Real >::Dot( Point3D< Real >( dpx( q[0] ) * _py * _pz , _px * dpy( q[1] ) * _pz , _px * _py * dpz( q[2] ) ) , p.gradient ) * gradientWeight + ) * p.weight; +#endif // POINT_DATA_RES + } + static inline Real _CalculateConstraint_( const PointData< Real , true >& p , const Polynomial< Degree >& px , const Polynomial< Degree >& py , const Polynomial< Degree >& pz , const Polynomial< Degree >& dpx , const Polynomial< Degree >& dpy , const Polynomial< Degree >& dpz ) + { +#if POINT_DATA_RES + Real constraint = 0; + for( int c=0 ; c::SAMPLES ; c++ ) if( p[c].weight ) + { + const Point3D< Real > q = p[c].position; + double _px = px( q[0] ) , _py = py( q[1] ) , _pz = pz( q[2] ); + constraint += + (Real)( _px * _py * _pz * p[c]._value ) + + Point3D< Real >::Dot( Point3D< Real >( dpx( q[0] ) * _py * _pz , _px * dpy( q[1] ) * _pz , _px * _py * dpz( q[2] ) ) , p[c]._gradient ); + } + return constraint; +#else // !POINT_DATA_RES + const Point3D< Real > q = p.position; + double _px = px( q[0] ) , _py = py( q[1] ) , _pz = pz( q[2] ); + return + (Real)( _px * _py * _pz * p._value ) + + Point3D< Real >::Dot( Point3D< Real >( dpx( q[0] ) * _py * _pz , _px * dpy( q[1] ) * _pz , _px * _py * dpz( q[2] ) ) , p._gradient ); +#endif // POINT_DATA_RES + } +#if POINT_DATA_RES + static inline void _CalculateCoarser_( int c , PointData< Real , true >& p , Real value , Point3D< Real > gradient , Real valueWeight , Real gradientWeight ){ p[c]._value = value * valueWeight * p[c].weight ; p[c]._gradient = gradient * gradientWeight * p[c].weight; } +#else // !POINT_DATA_RES + static inline void _CalculateCoarser_( PointData< Real , true >& p , Real value , Point3D< Real > gradient , Real valueWeight , Real gradientWeight ){ p._value = value * valueWeight * p.weight ; p._gradient = gradient * gradientWeight * p.weight; } +#endif // POINT_DATA_RES +}; + +template< > +template< class I > +double FEMSystemFunctor< 0 , BOUNDARY_FREE >::_integrate( const I& integrator , const int off1[] , const int off2[] ) const +{ +#define D_DOT( D1 , D2 ) { integrator.dot( off1[0] , off2[0] , D1 , D2 ) , integrator.dot( off1[1] , off2[1] , D1 , D2 ) , integrator.dot( off1[2] , off2[2] , D1 , D2 ) } + double d00[] = D_DOT( 0 , 0 ); + return + ( + d00[0] * d00[1] * d00[2] + ) * massWeight; +#undef D_DOT +} +template< > +template< class I > +double FEMSystemFunctor< 0 , BOUNDARY_NEUMANN >::_integrate( const I& integrator , const int off1[] , const int off2[] ) const +{ +#define D_DOT( D1 , D2 ) { integrator.dot( off1[0] , off2[0] , D1 , D2 ) , integrator.dot( off1[1] , off2[1] , D1 , D2 ) , integrator.dot( off1[2] , off2[2] , D1 , D2 ) } + double d00[] = D_DOT( 0 , 0 ); + return + ( + d00[0] * d00[1] * d00[2] + ) * massWeight; +#undef D_DOT +} +template< > +template< class I > +double FEMSystemFunctor< 0 , BOUNDARY_DIRICHLET >::_integrate( const I& integrator , const int off1[] , const int off2[] ) const +{ +#define D_DOT( D1 , D2 ) { integrator.dot( off1[0] , off2[0] , D1 , D2 ) , integrator.dot( off1[1] , off2[1] , D1 , D2 ) , integrator.dot( off1[2] , off2[2] , D1 , D2 ) } + double d00[] = D_DOT( 0 , 0 ); + return + ( + d00[0] * d00[1] * d00[2] + ) * massWeight; +#undef D_DOT +} +template< > +template< class I > +double FEMSystemFunctor< 1 , BOUNDARY_FREE >::_integrate( const I& integrator , const int off1[] , const int off2[] ) const +{ +#define D_DOT( D1 , D2 ) { integrator.dot( off1[0] , off2[0] , D1 , D2 ) , integrator.dot( off1[1] , off2[1] , D1 , D2 ) , integrator.dot( off1[2] , off2[2] , D1 , D2 ) } + double d00[] = D_DOT( 0 , 0 ) , d11[] = D_DOT( 1 , 1 ); + return + ( + d00[0] * d00[1] * d00[2] + ) * massWeight + + + ( + d11[0] * d00[1] * d00[2] + + d11[1] * d00[2] * d00[0] + + d11[2] * d00[0] * d00[1] + ) * lapWeight; +#undef D_DOT +} +template< > +template< class I > +double FEMSystemFunctor< 1 , BOUNDARY_NEUMANN >::_integrate( const I& integrator , const int off1[] , const int off2[] ) const +{ +#define D_DOT( D1 , D2 ) { integrator.dot( off1[0] , off2[0] , D1 , D2 ) , integrator.dot( off1[1] , off2[1] , D1 , D2 ) , integrator.dot( off1[2] , off2[2] , D1 , D2 ) } + double d00[] = D_DOT( 0 , 0 ) , d11[] = D_DOT( 1 , 1 ); + return + ( + d00[0] * d00[1] * d00[2] + ) * massWeight + + + ( + d11[0] * d00[1] * d00[2] + + d11[1] * d00[2] * d00[0] + + d11[2] * d00[0] * d00[1] + ) * lapWeight; +#undef D_DOT +} +template< > +template< class I > +double FEMSystemFunctor< 1 , BOUNDARY_DIRICHLET >::_integrate( const I& integrator , const int off1[] , const int off2[] ) const +{ +#define D_DOT( D1 , D2 ) { integrator.dot( off1[0] , off2[0] , D1 , D2 ) , integrator.dot( off1[1] , off2[1] , D1 , D2 ) , integrator.dot( off1[2] , off2[2] , D1 , D2 ) } + double d00[] = D_DOT( 0 , 0 ) , d11[] = D_DOT( 1 , 1 ); + return + ( + d00[0] * d00[1] * d00[2] + ) * massWeight + + + ( + d11[0] * d00[1] * d00[2] + + d11[1] * d00[2] * d00[0] + + d11[2] * d00[0] * d00[1] + ) * lapWeight; +#undef D_DOT +} + +template< int FEMDegree , BoundaryType BType > +template< class I > +double FEMSystemFunctor< FEMDegree , BType >::_integrate( const I& integrator , const int off1[] , const int off2[] ) const +{ +#define D_DOT( D1 , D2 ) { integrator.dot( off1[0] , off2[0] , D1 , D2 ) , integrator.dot( off1[1] , off2[1] , D1 , D2 ) , integrator.dot( off1[2] , off2[2] , D1 , D2 ) } + double d00[] = D_DOT( 0 , 0 ) , d02[] = D_DOT( 0 , 2 ) , d20[] = D_DOT( 2 , 0 ) , d22[] = D_DOT( 2 , 2 ) , d11[] = D_DOT( 1 , 1 ); + return + ( + d00[0] * d00[1] * d00[2] + ) * massWeight + + + ( + d11[0] * d00[1] * d00[2] + + d11[1] * d00[2] * d00[0] + + d11[2] * d00[0] * d00[1] + ) * lapWeight + + + ( + d22[0] * d00[1] * d00[2] + // Unmixed + d22[1] * d00[2] * d00[0] + // Unmixed + d22[2] * d00[0] * d00[1] + // Unmixed + d00[0] * ( d02[1] * d20[2] + d20[1] * d02[2] ) + // Mixed + d00[1] * ( d02[2] * d20[0] + d20[2] * d02[0] ) + // Mixed + d00[2] * ( d02[0] * d20[1] + d20[0] * d02[1] ) // Mixed + ) * biLapWeight; +#undef D_DOT +} +template< int SFDegree , BoundaryType SFBType , int FEMDegree , BoundaryType FEMBType > +template< bool Reverse , class I > +double FEMSFConstraintFunctor< SFDegree , SFBType , FEMDegree , FEMBType >::_integrate( const I& integrator , const int off1[] , const int off2[] ) const +{ +#define D_DOT( D1 , D2 ) { integrator.dot( off1[0] , off2[0] , Reverse ? D2 : D1 , Reverse ? D1 : D2 ) , integrator.dot( off1[1] , off2[1] , Reverse ? D2 : D1 , Reverse ? D1 : D2 ) , integrator.dot( off1[2] , off2[2] , Reverse ? D2 : D1 , Reverse ? D1 : D2 ) } + double d00[] = D_DOT( 0 , 0 ) , d02[] = D_DOT( 0 , 2 ) , d20[] = D_DOT( 2 , 0 ) , d22[] = D_DOT( 2 , 2 ) , d11[] = D_DOT( 1 , 1 ); + if( SFDegree==0 || FEMDegree==0 ) + return d00[0] * d00[1] * d00[2] * massWeight; + else if( SFDegree<=1 || FEMDegree<=1 ) + return + ( + d00[0] * d00[1] * d00[2] + ) * massWeight + + + ( + d11[0] * d00[1] * d00[2] + + d11[1] * d00[2] * d00[0] + + d11[2] * d00[0] * d00[1] + ) * lapWeight; + else + return + ( + d00[0] * d00[1] * d00[2] + ) * massWeight + + + ( + d11[0] * d00[1] * d00[2] + + d11[1] * d00[2] * d00[0] + + d11[2] * d00[0] * d00[1] + ) * lapWeight + + + ( + d22[0] * d00[1] * d00[2] + // Unmixed + d22[1] * d00[2] * d00[0] + // Unmixed + d22[2] * d00[0] * d00[1] + // Unmixed + d00[0] * ( d02[1] * d20[2] + d20[1] * d02[2] ) + // Mixed + d00[1] * ( d02[2] * d20[0] + d20[2] * d02[0] ) + // Mixed + d00[2] * ( d02[0] * d20[1] + d20[0] * d02[1] ) // Mixed + ) * biLapWeight; +#undef D_DOT +} +template< int VFDegree , BoundaryType VFBType , int FEMDegree , BoundaryType FEMBType > +template< bool Reverse , class I > +Point3D< double > FEMVFConstraintFunctor< VFDegree , VFBType , FEMDegree , FEMBType >::_integrate( const I& integrator , const int off1[] , const int off2[] ) const +{ +#define D_DOT( D1 , D2 ) { integrator.dot( off1[0] , off2[0] , Reverse ? D2 : D1 , Reverse ? D1 : D2 ) , integrator.dot( off1[1] , off2[1] , Reverse ? D2 : D1 , Reverse ? D1 : D2 ) , integrator.dot( off1[2] , off2[2] , Reverse ? D2 : D1 , Reverse ? D1 : D2 ) } + if( FEMDegree==0 ) fprintf( stderr , "[ERROR] FEMDegree does not support differentiation: %d\n" , FEMDegree ) , exit( 0 ); + if( VFDegree==0 || FEMDegree==1 ) + { + double d00[] = D_DOT( 0 , 0 ) , d01[] = D_DOT( 0 , 1 ); + return + Point3D< double > + ( + d01[0] * d00[1] * d00[2] , + d01[1] * d00[2] * d00[0] , + d01[2] * d00[0] * d00[1] + ) * lapWeight; + } + else + { + double d00[] = D_DOT( 0 , 0 ) , d10[] = D_DOT( 1 , 0 ) , d01[] = D_DOT( 0 , 1 ) , d02[] = D_DOT( 0 , 2 ) , d12[] = D_DOT( 1 , 2 ); + return + Point3D< double > + ( + d01[0] * d00[1] * d00[2] , + d01[1] * d00[2] * d00[0] , + d01[2] * d00[0] * d00[1] + ) * lapWeight + + + Point3D< double > + ( + d12[0] * d00[1] * d00[2] + d10[0] * ( d00[1] * d02[2] + d02[1] * d00[2] ) , + d12[1] * d00[2] * d00[0] + d10[1] * ( d00[2] * d02[0] + d02[2] * d00[0] ) , + d12[2] * d00[0] * d00[1] + d10[2] * ( d00[0] * d02[1] + d02[0] * d00[1] ) + ) * biLapWeight; + } +#undef D_DOT +} + +template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 > +template< bool Reverse , class _FEMSystemFunctor > +void SystemCoefficients< Degree1 , BType1 , Degree2 , BType2 >::SetCentralConstraintStencil( const _FEMSystemFunctor& F , const Integrator& integrator , Stencil< double , OverlapSize >& stencil ) +{ + int center = ( 1<>1; + int offset[] = { center , center , center }; + for( int x=0 ; x( integrator , _offset , offset ); + } +} +template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 > +template< bool Reverse , class _FEMSystemFunctor > +void SystemCoefficients< Degree1 , BType1 , Degree2 , BType2 >::SetCentralConstraintStencils( const _FEMSystemFunctor& F , const ChildIntegrator& integrator , Stencil< double , OverlapSize > stencils[2][2][2] ) +{ + int center = ( 1<>1; + // [NOTE] We want the center to be at the first node of the brood + // Which is not the case when childDepth is 1. + center = ( center>>1 )<<1; + for( int i=0 ; i<2 ; i++ ) for( int j=0 ; j<2 ; j++ ) for( int k=0 ; k<2 ; k++ ) + { + int offset[] = { center+i , center+j , center+k }; + for( int x=0 ; x( integrator , _offset , offset ); + } + } +} +template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 > +template< bool Reverse , class _FEMSystemFunctor > +void SystemCoefficients< Degree1 , BType1 , Degree2 , BType2 >::SetCentralConstraintStencil( const _FEMSystemFunctor& F , const Integrator& integrator , Stencil< Point3D< double > , OverlapSize >& stencil ) +{ + int center = ( 1<>1; + int offset[] = { center , center , center }; + for( int x=0 ; x( integrator , _offset , offset ); + } +} +template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 > +template< bool Reverse , class _FEMSystemFunctor > +void SystemCoefficients< Degree1 , BType1 , Degree2 , BType2 >::SetCentralConstraintStencils( const _FEMSystemFunctor& F , const ChildIntegrator& integrator , Stencil< Point3D< double > , OverlapSize > stencils[2][2][2] ) +{ + int center = ( 1<>1; + // [NOTE] We want the center to be at the first node of the brood + // Which is not the case when childDepth is 1. + center = ( center>>1 )<<1; + for( int i=0 ; i<2 ; i++ ) for( int j=0 ; j<2 ; j++ ) for( int k=0 ; k<2 ; k++ ) + { + int offset[] = { center+i , center+j , center+k }; + for( int x=0 ; x( integrator , _offset , offset ); + } + } +} +template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 > +template< class _FEMSystemFunctor > +void SystemCoefficients< Degree1 , BType1 , Degree2 , BType2 >::SetCentralSystemStencil( const _FEMSystemFunctor& F , const Integrator& integrator , Stencil< double , OverlapSize >& stencil ) +{ + int center = ( 1<>1; + int offset[] = { center , center , center }; + for( int x=0 ; x +template< class _FEMSystemFunctor > +void SystemCoefficients< Degree1 , BType1 , Degree2 , BType2 >::SetCentralSystemStencils( const _FEMSystemFunctor& F , const ChildIntegrator& integrator , Stencil< double , OverlapSize > stencils[2][2][2] ) +{ + int center = ( 1<>1; + // [NOTE] We want the center to be at the first node of the brood + // Which is not the case when childDepth is 1. + center = ( center>>1 )<<1; + for( int i=0 ; i<2 ; i++ ) for( int j=0 ; j<2 ; j++ ) for( int k=0 ; k<2 ; k++ ) + { + int offset[] = { center+i , center+j , center+k }; + for( int x=0 ; x +template< int FEMDegree > +void Octree< Real >::_setMultiColorIndices( int start , int end , std::vector< std::vector< int > >& indices ) const +{ + static const int OverlapRadius = - BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + + const int modulus = OverlapRadius+1; + indices.resize( modulus*modulus*modulus ); + int count[modulus*modulus*modulus]; + memset( count , 0 , sizeof(int)*modulus*modulus*modulus ); +#pragma omp parallel for num_threads( threads ) + for( int i=start ; idepthAndOffset( d , off ); + int idx = (modulus*modulus) * ( off[2]%modulus ) + modulus * ( off[1]%modulus ) + ( off[0]%modulus ); +#pragma omp atomic + count[idx]++; + } + + for( int i=0 ; idepthAndOffset( d , off ); + int idx = (modulus*modulus) * ( off[2]%modulus ) + modulus * ( off[1]%modulus ) + ( off[0]%modulus ); + indices[idx].push_back( i - start ); + } +} + +template< class Real > +template< class C , int FEMDegree , BoundaryType BType > +void Octree< Real >::_downSample( LocalDepth highDepth , DenseNodeData< C , FEMDegree >& constraints ) const +{ + typedef typename TreeOctNode::NeighborKey< -BSplineSupportSizes< FEMDegree >::UpSampleStart , BSplineSupportSizes< FEMDegree >::UpSampleEnd > UpSampleKey; + + LocalDepth lowDepth = highDepth-1; + if( lowDepth<0 ) return; + + typename BSplineEvaluationData< FEMDegree , BType >::UpSampleEvaluator upSampleEvaluator; + BSplineEvaluationData< FEMDegree , BType >::SetUpSampleEvaluator( upSampleEvaluator , lowDepth ); + std::vector< UpSampleKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i::UpSampleSize > upSampleStencil; + int lowCenter = ( 1<>1; + for( int i=0 ; i::UpSampleSize ; i++ ) for( int j=0 ; j::UpSampleSize ; j++ ) for( int k=0 ; k::UpSampleSize ; k++ ) + upSampleStencil( i , j , k ) = + upSampleEvaluator.value( lowCenter , 2*lowCenter + i + BSplineSupportSizes< FEMDegree >::UpSampleStart ) * + upSampleEvaluator.value( lowCenter , 2*lowCenter + j + BSplineSupportSizes< FEMDegree >::UpSampleStart ) * + upSampleEvaluator.value( lowCenter , 2*lowCenter + k + BSplineSupportSizes< FEMDegree >::UpSampleStart ); + + // Iterate over all (valid) parent nodes +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodesBegin(lowDepth) ; i<_sNodesEnd(lowDepth) ; i++ ) if( _isValidFEMNode( _sNodes.treeNodes[i] ) ) + { + TreeOctNode* pNode = _sNodes.treeNodes[i]; + + UpSampleKey& neighborKey = neighborKeys[ omp_get_thread_num() ]; + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( pNode , d , off ); + + neighborKey.template getNeighbors< false >( pNode ); + + // Get the child neighbors + typename TreeOctNode::Neighbors< BSplineSupportSizes< FEMDegree >::UpSampleSize > neighbors; + neighborKey.template getChildNeighbors< false >( 0 , _localToGlobal( d ) , neighbors ); + + C& coarseConstraint = constraints[i]; + + // Want to make sure test if contained children are interior. + // This is more conservative because we are test that overlapping children are interior + bool isInterior = _isInteriorlyOverlapped< FEMDegree , FEMDegree >( pNode ); + if( isInterior ) + { + for( int ii=0 ; ii::UpSampleSize ; ii++ ) for( int jj=0 ; jj::UpSampleSize ; jj++ ) for( int kk=0 ; kk::UpSampleSize ; kk++ ) + { + const TreeOctNode* cNode = neighbors.neighbors[ii][jj][kk]; + if( IsActiveNode( cNode ) ) coarseConstraint += (C)( constraints[ cNode->nodeData.nodeIndex ] * upSampleStencil( ii , jj , kk ) ); + } + } + else + { + double upSampleValues[3][ BSplineSupportSizes< FEMDegree >::UpSampleSize ]; + for( int ii=0 ; ii::UpSampleSize ; ii++ ) + { + upSampleValues[0][ii] = upSampleEvaluator.value( off[0] , 2*off[0] + ii + BSplineSupportSizes< FEMDegree >::UpSampleStart ); + upSampleValues[1][ii] = upSampleEvaluator.value( off[1] , 2*off[1] + ii + BSplineSupportSizes< FEMDegree >::UpSampleStart ); + upSampleValues[2][ii] = upSampleEvaluator.value( off[2] , 2*off[2] + ii + BSplineSupportSizes< FEMDegree >::UpSampleStart ); + } + for( int ii=0 ; ii::UpSampleSize ; ii++ ) for( int jj=0 ; jj::UpSampleSize ; jj++ ) + { + double dxy = upSampleValues[0][ii] * upSampleValues[1][jj]; + for( int kk=0 ; kk::UpSampleSize ; kk++ ) + { + const TreeOctNode* cNode = neighbors.neighbors[ii][jj][kk]; + if( _isValidFEMNode( cNode ) ) coarseConstraint += (C)( constraints[ cNode->nodeData.nodeIndex ] * dxy * upSampleValues[2][kk] ); + } + } + } + } +} +template< class Real > +template< class C , int FEMDegree , BoundaryType BType > +void Octree< Real >::_upSample( LocalDepth highDepth , DenseNodeData< C , FEMDegree >& coefficients ) const +{ + static const int LeftDownSampleRadius = -( ( BSplineSupportSizes< FEMDegree >::DownSample0Start < BSplineSupportSizes< FEMDegree >::DownSample1Start ) ? BSplineSupportSizes< FEMDegree >::DownSample0Start : BSplineSupportSizes< FEMDegree >::DownSample1Start ); + static const int RightDownSampleRadius = ( ( BSplineSupportSizes< FEMDegree >::DownSample0End > BSplineSupportSizes< FEMDegree >::DownSample1End ) ? BSplineSupportSizes< FEMDegree >::DownSample0End : BSplineSupportSizes< FEMDegree >::DownSample1End ); + typedef TreeOctNode::NeighborKey< LeftDownSampleRadius , RightDownSampleRadius > DownSampleKey; + + LocalDepth lowDepth = highDepth-1; + if( lowDepth<0 ) return; + + typename BSplineEvaluationData< FEMDegree , BType >::UpSampleEvaluator upSampleEvaluator; + BSplineEvaluationData< FEMDegree , BType >::SetUpSampleEvaluator( upSampleEvaluator , lowDepth ); + std::vector< DownSampleKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i::DownSample0Size > BSplineSupportSizes< FEMDegree >::DownSample1Size ? BSplineSupportSizes< FEMDegree >::DownSample0Size : BSplineSupportSizes< FEMDegree >::DownSample1Size; + Stencil< double , DownSampleSize > downSampleStencils[ Cube::CORNERS ]; + int lowCenter = ( 1<>1; + for( int c=0 ; c::DownSampleSize[cx] ; ii++ ) + for( int jj=0 ; jj::DownSampleSize[cy] ; jj++ ) + for( int kk=0 ; kk::DownSampleSize[cz] ; kk++ ) + downSampleStencils[c]( ii , jj , kk ) = + upSampleEvaluator.value( lowCenter + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] , 2*lowCenter + cx ) * + upSampleEvaluator.value( lowCenter + jj + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] , 2*lowCenter + cy ) * + upSampleEvaluator.value( lowCenter + kk + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] , 2*lowCenter + cz ) ; + } + + // For Dirichlet constraints, can't get to all children from parents because boundary nodes are invalid +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodesBegin(highDepth) ; i<_sNodesEnd(highDepth) ; i++ ) if( _isValidFEMNode( _sNodes.treeNodes[i] ) ) + { + TreeOctNode *cNode = _sNodes.treeNodes[i] , *pNode = cNode->parent; + int c = (int)( cNode-pNode->children ); + + DownSampleKey& neighborKey = neighborKeys[ omp_get_thread_num() ]; + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( pNode , d , off ); + typename TreeOctNode::Neighbors< LeftDownSampleRadius + RightDownSampleRadius + 1 >& neighbors = neighborKey.template getNeighbors< false >( pNode ); + + // Want to make sure test if contained children are interior. + // This is more conservative because we are test that overlapping children are interior + bool isInterior = _isInteriorlyOverlapped< FEMDegree , FEMDegree >( pNode ); + + C& fineCoefficient = coefficients[ cNode->nodeData.nodeIndex ]; + + int cx , cy , cz; + Cube::FactorCornerIndex( c , cx , cy , cz ); + + if( isInterior ) + { + for( int ii=0 ; ii::DownSampleSize[cx] ; ii++ ) for( int jj=0 ; jj::DownSampleSize[cy] ; jj++ ) + { + int _ii = ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] + LeftDownSampleRadius; + int _jj = jj + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] + LeftDownSampleRadius; + for( int kk=0 ; kk::DownSampleSize[cz] ; kk++ ) + { + int _kk = kk + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] + LeftDownSampleRadius; + const TreeOctNode* _pNode = neighbors.neighbors[_ii][_jj][_kk]; + if( _pNode ) fineCoefficient += (C)( coefficients[ _pNode->nodeData.nodeIndex ] * downSampleStencils[c]( ii , jj , kk ) ); + } + } + } + else + { + double downSampleValues[3][ BSplineSupportSizes< FEMDegree >::DownSample0Size > BSplineSupportSizes< FEMDegree >::DownSample1Size ? BSplineSupportSizes< FEMDegree >::DownSample0Size : BSplineSupportSizes< FEMDegree >::DownSample1Size ]; + for( int ii=0 ; ii::DownSampleSize[cx] ; ii++ ) downSampleValues[0][ii] = upSampleEvaluator.value( off[0] + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] , 2*off[0] + cx ); + for( int ii=0 ; ii::DownSampleSize[cy] ; ii++ ) downSampleValues[1][ii] = upSampleEvaluator.value( off[1] + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] , 2*off[1] + cy ); + for( int ii=0 ; ii::DownSampleSize[cz] ; ii++ ) downSampleValues[2][ii] = upSampleEvaluator.value( off[2] + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] , 2*off[2] + cz ); + + for( int ii=0 ; ii::DownSampleSize[cx] ; ii++ ) for( int jj=0 ; jj::DownSampleSize[cy] ; jj++ ) + { + double dxy = downSampleValues[0][ii] * downSampleValues[1][jj]; + int _ii = ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] + LeftDownSampleRadius; + int _jj = jj + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] + LeftDownSampleRadius; + for( int kk=0 ; kk::DownSampleSize[cz] ; kk++ ) + { + int _kk = kk + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] + LeftDownSampleRadius; + const TreeOctNode* _pNode = neighbors.neighbors[_ii][_jj][_kk]; + if( _isValidFEMNode( _pNode ) ) fineCoefficient += (C)( coefficients[ _pNode->nodeData.nodeIndex ] * dxy * downSampleValues[2][kk] ); + } + } + } + } +} + +template< class Real > +template< class C , int FEMDegree , BoundaryType BType > +void Octree< Real >::_UpSample( LocalDepth highDepth , ConstPointer( C ) lowCoefficients , Pointer( C ) highCoefficients , int threads ) +{ + static const int LeftDownSampleRadius = -( ( BSplineSupportSizes< FEMDegree >::DownSample0Start < BSplineSupportSizes< FEMDegree >::DownSample1Start ) ? BSplineSupportSizes< FEMDegree >::DownSample0Start : BSplineSupportSizes< FEMDegree >::DownSample1Start ); + static const int RightDownSampleRadius = ( ( BSplineSupportSizes< FEMDegree >::DownSample0End > BSplineSupportSizes< FEMDegree >::DownSample1End ) ? BSplineSupportSizes< FEMDegree >::DownSample0End : BSplineSupportSizes< FEMDegree >::DownSample1End ); + typedef TreeOctNode::NeighborKey< LeftDownSampleRadius , RightDownSampleRadius > DownSampleKey; + + LocalDepth lowDepth = highDepth - 1; + if( lowDepth<0 ) return; + + typename BSplineEvaluationData< FEMDegree , BType >::UpSampleEvaluator upSampleEvaluator; + BSplineEvaluationData< FEMDegree , BType >::SetUpSampleEvaluator( upSampleEvaluator , lowDepth ); + std::vector< DownSampleKey > neighborKeys( std::max< int >( 1 , threads ) ); + + static const int DownSampleSize = BSplineSupportSizes< FEMDegree >::DownSample0Size > BSplineSupportSizes< FEMDegree >::DownSample1Size ? BSplineSupportSizes< FEMDegree >::DownSample0Size : BSplineSupportSizes< FEMDegree >::DownSample1Size; + Stencil< double , DownSampleSize > downSampleStencils[ Cube::CORNERS ]; + int lowCenter = ( 1<>1; + for( int c=0 ; c::DownSample0Size > BSplineSupportSizes< FEMDegree >::DownSample1Size ? BSplineSupportSizes< FEMDegree >::DownSample0Size : BSplineSupportSizes< FEMDegree >::DownSample1Size; + for( int ii=0 ; ii::DownSampleSize[cx] ; ii++ ) + for( int jj=0 ; jj::DownSampleSize[cy] ; jj++ ) + for( int kk=0 ; kk::DownSampleSize[cz] ; kk++ ) + downSampleStencils[c]( ii , jj , kk ) = + upSampleEvaluator.value( lowCenter + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] , 2*lowCenter + cx ) * + upSampleEvaluator.value( lowCenter + jj + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] , 2*lowCenter + cy ) * + upSampleEvaluator.value( lowCenter + kk + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] , 2*lowCenter + cz ) ; + } + int lowBegin = _BSplineBegin< FEMDegree , BType >( lowDepth ) , lowEnd = _BSplineEnd< FEMDegree , BType >( lowDepth ); + int highBegin = _BSplineBegin< FEMDegree , BType >( highDepth ) , highEnd = _BSplineEnd< FEMDegree , BType >( highDepth ); + int lowDim = lowEnd - lowBegin , highDim = highEnd - highBegin; + // Iterate over all child nodes. (This is required since there can be child nodes whose parent is inactive.) +#pragma omp parallel for num_threads( threads ) + for( int k=0 ; k>1 , _off[1] = off[1]>>1 , _off[2] = off[2]>>1; + + // Want to make sure test if contained children are interior. + // This is more conservative because we are test that overlapping children are interior + bool isInterior = _IsInteriorlyOverlapped< FEMDegree , FEMDegree >( lowDepth , _off ); + int cx = off[0]&1 , cy = off[1]&1 , cz = off[2]&1; + int c = Cube::CornerIndex( cx , cy , cz ); + + C& highCoefficient = highCoefficients[ highIdx ]; + + if( isInterior ) + { + for( int ii=0 ; ii::DownSampleSize[cx] ; ii++ ) for( int jj=0 ; jj::DownSampleSize[cy] ; jj++ ) + { + int _i = _off[0] + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] - lowBegin; + int _j = _off[1] + jj + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] - lowBegin; + for( int kk=0 ; kk::DownSampleSize[cz] ; kk++ ) + { + int _k = _off[2] + kk + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] - lowBegin; + highCoefficient += (C)( lowCoefficients[ _i + _j*lowDim + _k*lowDim*lowDim ] * downSampleStencils[c]( ii , jj , kk ) ); + } + } + } + else + { + double downSampleValues[3][ BSplineSupportSizes< FEMDegree >::DownSample0Size > BSplineSupportSizes< FEMDegree >::DownSample1Size ? BSplineSupportSizes< FEMDegree >::DownSample0Size : BSplineSupportSizes< FEMDegree >::DownSample1Size ]; + + for( int ii=0 ; ii::DownSampleSize[cx] ; ii++ ) downSampleValues[0][ii] = upSampleEvaluator.value( _off[0] + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] , off[0] ); + for( int ii=0 ; ii::DownSampleSize[cy] ; ii++ ) downSampleValues[1][ii] = upSampleEvaluator.value( _off[1] + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] , off[1] ); + for( int ii=0 ; ii::DownSampleSize[cz] ; ii++ ) downSampleValues[2][ii] = upSampleEvaluator.value( _off[2] + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] , off[2] ); + + for( int ii=0 ; ii::DownSampleSize[cx] ; ii++ ) for( int jj=0 ; jj::DownSampleSize[cy] ; jj++ ) + { + double dxy = downSampleValues[0][ii] * downSampleValues[1][jj]; + int _i = _off[0] + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] - lowBegin; + int _j = _off[1] + jj + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] - lowBegin; + if( _i>=0 && _i=0 && _j::DownSampleSize[cz] ; kk++ ) + { + int _k = _off[2] + kk + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] - lowBegin; + if( _k>=0 && _k +template< class C , int FEMDegree , BoundaryType BType > +DenseNodeData< C , FEMDegree > Octree< Real >::coarseCoefficients( const DenseNodeData< C , FEMDegree >& coefficients ) const +{ + DenseNodeData< Real , FEMDegree > coarseCoefficients( _sNodesEnd(_maxDepth-1) ); + memset( &coarseCoefficients[0] , 0 , sizeof(Real)*_sNodesEnd(_maxDepth-1) ); +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodesBegin(0) ; i<_sNodesEnd(_maxDepth-1) ; i++ ) coarseCoefficients[i] = coefficients[i]; + for( LocalDepth d=1 ; d<_maxDepth ; d++ ) _upSample< C , FEMDegree , BType >( d , coarseCoefficients ); + return coarseCoefficients; +} +template< class Real > +template< class C , int FEMDegree , BoundaryType BType > +DenseNodeData< C , FEMDegree > Octree< Real >::coarseCoefficients( const SparseNodeData< C , FEMDegree >& coefficients ) const +{ + DenseNodeData< Real , FEMDegree > coarseCoefficients( _sNodesEnd(_maxDepth-1) ); + memset( &coarseCoefficients[0] , 0 , sizeof(Real)*_sNodesEnd(_maxDepth-1) ); +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodesBegin(0) ; i<_sNodesEnd(_maxDepth-1) ; i++ ) + { + const C* c = coefficients( _sNodes.treeNodes[i] ); + if( c ) coarseCoefficients[i] = *c; + } + for( LocalDepth d=1 ; d<_maxDepth ; d++ ) _upSample< C , FEMDegree , BType >( d , coarseCoefficients ); + return coarseCoefficients; +} + +template< class Real > +template< int FEMDegree , BoundaryType BType > +Real Octree< Real >::_coarserFunctionValue( Point3D< Real > p , const PointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* pointNode , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& upSampledCoefficients ) const +{ + static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize; + static const int LeftSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart; + + double pointValue = 0; + LocalDepth depth = _localDepth( pointNode ); + if( depth<0 ) return (Real)0.; + + // Iterate over all basis functions that overlap the point at the coarser resolution + { + const typename TreeOctNode::Neighbors< SupportSize >& neighbors = neighborKey.neighbors[ _localToGlobal( depth-1 ) ]; + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( pointNode->parent , _d , _off ); + int fStart , fEnd; + BSplineData< FEMDegree , BType >::FunctionSpan( _d , fStart , fEnd ); + + double pointValues[ DIMENSION ][SupportSize]; + memset( pointValues , 0 , sizeof(double) * DIMENSION * SupportSize ); + + for( int dd=0 ; dd::FunctionIndex( _d , _off[dd]+i ); + if( fIdx>=fStart && fIdxnodeData.nodeIndex] ); + } + pointValue += _pointValue * xyValue; + } + } + return Real( pointValue ); +} +template< class Real > +template< int FEMDegree , BoundaryType BType > +Point3D< Real > Octree< Real >::_coarserFunctionGradient( Point3D< Real > p , const PointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* pointNode , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& upSampledCoefficients ) const +{ + static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize; + static const int LeftSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart; + + Point3D< double > pointGradient; + LocalDepth depth = _localDepth( pointNode ); + if( depth<=0 ) return Real(0.); + + // Iterate over all basis functions that overlap the point at the coarser resolution + { + const typename TreeOctNode::Neighbors< SupportSize >& neighbors = neighborKey.neighbors[ _localToGlobal( depth-1 ) ]; + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( pointNode->parent , _d , _off ); + int fStart , fEnd; + BSplineData< FEMDegree , BType >::FunctionSpan( _d , fStart , fEnd ); + + double _pointValues[ DIMENSION ][SupportSize] , dPointValues[ DIMENSION ][SupportSize]; + memset( _pointValues , 0 , sizeof(double) * DIMENSION * SupportSize ); + memset( dPointValues , 0 , sizeof(double) * DIMENSION * SupportSize ); + + for( int dd=0 ; dd::FunctionIndex( _d , _off[dd]+i ); + if( fIdx>=fStart && fIdxnodeData.nodeIndex] ); + _dPointValue += dPointValues[2][l] * double( upSampledCoefficients[_node->nodeData.nodeIndex] ); + } + } + + pointGradient += Point3D< double >( __pointValue * dx_yValue , __pointValue * _xdyValue , _dPointValue * _x_yValue ); + } + } + return Point3D< Real >( pointGradient ); +} + +template< class Real > +template< int FEMDegree , BoundaryType BType > +Real Octree< Real >::_finerFunctionValue( Point3D< Real > p , const PointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* pointNode , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& finerCoefficients ) const +{ + typename TreeOctNode::Neighbors< BSplineSupportSizes< FEMDegree >::SupportSize > childNeighbors; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + + double pointValue = 0; + LocalDepth depth = _localDepth( pointNode ); + neighborKey.template getChildNeighbors< false >( _childIndex( pointNode , p ) , _localToGlobal( depth ) , childNeighbors ); + for( int j=-LeftPointSupportRadius ; j<=RightPointSupportRadius ; j++ ) + for( int k=-LeftPointSupportRadius ; k<=RightPointSupportRadius ; k++ ) + for( int l=-LeftPointSupportRadius ; l<=RightPointSupportRadius ; l++ ) + { + const TreeOctNode* _node = childNeighbors.neighbors[j+LeftPointSupportRadius][k+LeftPointSupportRadius][l+LeftPointSupportRadius]; + if( _isValidFEMNode( _node ) ) + { + int fIdx[3]; + functionIndex< FEMDegree , BType >( _node , fIdx ); + pointValue += + bsData.baseBSplines[ fIdx[0] ][LeftSupportRadius-j]( p[0] ) * + bsData.baseBSplines[ fIdx[1] ][LeftSupportRadius-k]( p[1] ) * + bsData.baseBSplines[ fIdx[2] ][LeftSupportRadius-l]( p[2] ) * + double( finerCoefficients[ _node->nodeData.nodeIndex ] ); + } + } + return Real( pointValue ); +} +template< class Real > +template< int FEMDegree , BoundaryType BType > +Point3D< Real > Octree< Real >::_finerFunctionGradient( Point3D< Real > p , const PointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* pointNode , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& finerCoefficients ) const +{ + typename TreeOctNode::Neighbors< BSplineSupportSizes< FEMDegree >::SupportSize > childNeighbors; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + + Point3D< double > pointGradient = 0; + LocalDepth depth = _localDepth( pointNode ); + neighborKey.template getChildNeighbors< false >( _childIndex( pointNode , p ) , _localToGlobal( depth ) , childNeighbors ); + for( int j=-LeftPointSupportRadius ; j<=RightPointSupportRadius ; j++ ) + for( int k=-LeftPointSupportRadius ; k<=RightPointSupportRadius ; k++ ) + for( int l=-LeftPointSupportRadius ; l<=RightPointSupportRadius ; l++ ) + { + const TreeOctNode* _node = childNeighbors.neighbors[j+LeftPointSupportRadius][k+LeftPointSupportRadius][l+LeftPointSupportRadius]; + if( _isValidFEMNode( _node ) ) + { + int fIdx[3]; + functionIndex< FEMDegree , BType >( _node , fIdx ); + double x = bsData. baseBSplines[ fIdx[0] ][LeftSupportRadius-j]( p[0] ) , y = bsData. baseBSplines[ fIdx[1] ][LeftSupportRadius-k]( p[1] ) , z = bsData. baseBSplines[ fIdx[2] ][LeftSupportRadius-l]( p[2] ); + double dx = bsData.dBaseBSplines[ fIdx[0] ][LeftSupportRadius-j]( p[0] ) , dy = bsData.dBaseBSplines[ fIdx[1] ][LeftSupportRadius-k]( p[1] ) , dz = bsData.dBaseBSplines[ fIdx[2] ][LeftSupportRadius-l]( p[2] ); + pointGradient += Point3D< double >( dx * y * z , x * dy * z , x * y * dz ) * (double)( finerCoefficients[ _node->nodeData.nodeIndex ] ); + } + } + return Point3D< Real >( pointGradient ); +} + +template< class Real > +template< int FEMDegree , BoundaryType BType , bool HasGradients > +void Octree< Real >::_setPointValuesFromCoarser( InterpolationInfo< HasGradients >& interpolationInfo , LocalDepth highDepth , const BSplineData< FEMDegree , BType >& bsData , const DenseNodeData< Real , FEMDegree >& upSampledCoefficients ) +{ + LocalDepth lowDepth = highDepth-1; + if( lowDepth<0 ) return; + std::vector< PointSupportKey< FEMDegree > > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i& neighborKey = neighborKeys[ omp_get_thread_num() ]; + PointData< Real , HasGradients >* pData = interpolationInfo( _sNodes.treeNodes[i] ); + if( pData ) + { + neighborKey.template getNeighbors< false >( _sNodes.treeNodes[i]->parent ); +#if POINT_DATA_RES + for( int c=0 ; c::SAMPLES ; c++ ) if( (*pData)[c].weight ) + _ConstraintCalculator_< Real , FEMDegree , HasGradients >::_CalculateCoarser_ + ( + c , *pData , + _coarserFunctionValue( (*pData)[c].position , neighborKey , _sNodes.treeNodes[i] , bsData , upSampledCoefficients ) , + HasGradients ? _coarserFunctionGradient( (*pData)[c].position , neighborKey , _sNodes.treeNodes[i] , bsData , upSampledCoefficients ) : Point3D< Real >() , + interpolationInfo.valueWeight , interpolationInfo.gradientWeight + ); +#else // !POINT_DATA_RES + _ConstraintCalculator_< Real , FEMDegree , HasGradients >::_CalculateCoarser_ + ( + *pData , + _coarserFunctionValue( pData->position , neighborKey , _sNodes.treeNodes[i] , bsData , upSampledCoefficients ) , + HasGradients ? _coarserFunctionGradient( pData->position , neighborKey , _sNodes.treeNodes[i] , bsData , upSampledCoefficients ) : Point3D< Real >() , + interpolationInfo.valueWeight , interpolationInfo.gradientWeight + ); +#endif // POINT_DATA_RES + } + } +} + +template< class Real > +template< int FEMDegree , BoundaryType BType , bool HasGradients > +void Octree< Real >::_updateCumulativeInterpolationConstraintsFromFiner( const InterpolationInfo< HasGradients >& interpolationInfo , const BSplineData< FEMDegree , BType >& bsData , LocalDepth highDepth , const DenseNodeData< Real , FEMDegree >& finerCoefficients , DenseNodeData< Real , FEMDegree >& coarserConstraints ) const +{ + static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + + // Note: We can't iterate over the finer point nodes as the point weights might be + // scaled incorrectly, due to the adaptive exponent. So instead, we will iterate + // over the coarser nodes and evaluate the finer solution at the associated points. + LocalDepth lowDepth = highDepth-1; + if( lowDepth<0 ) return; + size_t start = _sNodesBegin(lowDepth) , end = _sNodesEnd(lowDepth); + std::vector< PointSupportKey< FEMDegree > > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i& neighborKey = neighborKeys[ omp_get_thread_num() ]; + const PointData< Real , HasGradients >* pData = interpolationInfo( _sNodes.treeNodes[i] ); + if( pData ) + { + typename TreeOctNode::Neighbors< SupportSize >& neighbors = neighborKey.template getNeighbors< false >( _sNodes.treeNodes[i] ); + // evaluate the solution @( depth ) at the current point @( depth-1 ) +#if POINT_DATA_RES + for( int c=0 ; c::SAMPLES ; c++ ) if( (*pData)[c].weight ) +#endif // POINT_DATA_RES + { +#if POINT_DATA_RES + Real finerPointDValue = _finerFunctionValue( (*pData)[c].position , neighborKey , _sNodes.treeNodes[i] , bsData , finerCoefficients ) * interpolationInfo.valueWeight * (*pData)[c].weight; + Point3D< Real > finerPointDGradient = HasGradients ? _finerFunctionGradient( (*pData)[c].position , neighborKey , _sNodes.treeNodes[i] , bsData , finerCoefficients ) * interpolationInfo.gradientWeight * (*pData)[c].weight : Point3D< Real >(); + Point3D< Real > p = (*pData)[c].position; +#else // !POINT_DATA_RES + Real finerPointDValue = _finerFunctionValue( pData->position , neighborKey , _sNodes.treeNodes[i] , bsData , finerCoefficients ) * interpolationInfo.valueWeight * pData->weight; + Point3D< Real > finerPointDGradient = HasGradients ? _finerFunctionGradient( pData->position , neighborKey , _sNodes.treeNodes[i] , bsData , finerCoefficients ) * interpolationInfo.gradientWeight * pData->weight : Point3D< Real >(); + Point3D< Real > p = pData->position; +#endif // POINT_DATA_RES + // Update constraints for all nodes @( depth-1 ) that overlap the point + int idx[3]; + functionIndex< FEMDegree , BType >( _sNodes.treeNodes[i] , idx ); + for( int x=-LeftPointSupportRadius ; x<=RightPointSupportRadius ; x++ ) for( int y=-LeftPointSupportRadius ; y<=RightPointSupportRadius ; y++ ) for( int z=-LeftPointSupportRadius ; z<=RightPointSupportRadius ; z++ ) + { + const TreeOctNode* _node = neighbors.neighbors[x+LeftPointSupportRadius][y+LeftPointSupportRadius][z+LeftPointSupportRadius]; + if( _isValidFEMNode( _node ) ) + { + double px = bsData.baseBSplines[idx[0]+x][LeftSupportRadius-x]( p[0] ) , py = bsData.baseBSplines[idx[1]+y][LeftSupportRadius-y]( p[1] ) , pz = bsData.baseBSplines[idx[2]+z][LeftSupportRadius-z]( p[2] ); +#pragma omp atomic + coarserConstraints[ _node->nodeData.nodeIndex ] += (Real)( px * py * pz * finerPointDValue ); + if( HasGradients ) + { + double dpx = bsData.dBaseBSplines[idx[0]+x][LeftSupportRadius-x]( p[0] ) , dpy = bsData.dBaseBSplines[idx[1]+y][LeftSupportRadius-y]( p[1] ) , dpz = bsData.dBaseBSplines[idx[2]+z][LeftSupportRadius-z]( p[2] ); +#pragma omp atomic + coarserConstraints[ _node->nodeData.nodeIndex ] += Point3D< Real >::Dot( finerPointDGradient , Point3D< Real >( dpx * py * pz , px * dpy * pz , px * py * dpz ) ); + } + } + } + } + } + } +} + +template< class Real > +template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients > +int Octree< Real >::_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 +{ + static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize; + static const int OverlapRadius = - BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + static const int OverlapSize = BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize; + static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int RightPointSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + + bool hasYZPoints[SupportSize] , hasZPoints[SupportSize][SupportSize]; + Real diagonal = 0; + // Given a node: + // -- for each node in its support: + // ---- if the supporting node contains a point: + // ------ evaluate the x, y, and z B-splines of the nodes supporting the point + // splineValues \in [-LeftSupportRadius,RightSupportRadius] x [-LeftSupportRadius,RightSupportRadius] x [-LeftSupportRadius,RightSupportRadius] x [0,Dimension) x [-LeftPointSupportRadius,RightPointSupportRadius] +#if POINT_DATA_RES + Real _splineValues[PointData< Real , HasGradients >::SAMPLES][SupportSize][SupportSize][SupportSize][DIMENSION][SupportSize]; + Real wSplineValues[PointData< Real , HasGradients >::SAMPLES][SupportSize][SupportSize][SupportSize][DIMENSION][SupportSize]; + Real dSplineValues[PointData< Real , HasGradients >::SAMPLES][SupportSize][SupportSize][SupportSize][DIMENSION][SupportSize]; + memset( _splineValues , 0 , sizeof( Real ) * PointData< Real , HasGradients >::SAMPLES * SupportSize * SupportSize * SupportSize * DIMENSION *SupportSize ); + memset( wSplineValues , 0 , sizeof( Real ) * PointData< Real , HasGradients >::SAMPLES * SupportSize * SupportSize * SupportSize * DIMENSION *SupportSize ); + memset( dSplineValues , 0 , sizeof( Real ) * PointData< Real , HasGradients >::SAMPLES * SupportSize * SupportSize * SupportSize * DIMENSION *SupportSize ); +#else // !POINT_DATA_RES + Real _splineValues[SupportSize][SupportSize][SupportSize][DIMENSION][SupportSize]; + Real wSplineValues[SupportSize][SupportSize][SupportSize][DIMENSION][SupportSize]; + Real dSplineValues[SupportSize][SupportSize][SupportSize][DIMENSION][SupportSize]; + memset( _splineValues , 0 , sizeof( Real ) * SupportSize * SupportSize * SupportSize * DIMENSION *SupportSize ); + memset( wSplineValues , 0 , sizeof( Real ) * SupportSize * SupportSize * SupportSize * DIMENSION *SupportSize ); + memset( dSplineValues , 0 , sizeof( Real ) * SupportSize * SupportSize * SupportSize * DIMENSION *SupportSize ); +#endif // NEW_POINT_DATA + + int count = 0; + const TreeOctNode* node = neighbors.neighbors[OverlapRadius][OverlapRadius][OverlapRadius]; + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( node , d , off ); + int fStart , fEnd; + BSplineData< FEMDegree , BType >::FunctionSpan( d , fStart , fEnd ); + bool isInterior = _isInteriorlyOverlapped< FEMDegree , FEMDegree >( node ); + + if( interpolationInfo ) + { + // Iterate over all neighboring nodes that may have a constraining point + // -- For each one, compute the values of the spline functions supported on the point + for( int j=0 ; j& pData = *( (*interpolationInfo)( _node ) ); + +#if POINT_DATA_RES + for( int c=0 ; c::SAMPLES ; c++ ) if( pData[c].weight ) +#endif // POINT_DATA_RES + { +#if POINT_DATA_RES + Real (*__splineValues)[SupportSize] = _splineValues[c][jj][kk][ll]; + Real (*_wSplineValues)[SupportSize] = wSplineValues[c][jj][kk][ll]; + Real (*_dSplineValues)[SupportSize] = dSplineValues[c][jj][kk][ll]; + Real weight = pData[c].weight; + Point3D< Real > p = pData[c].position; +#else // !POINT_DATA_RES + Real (*__splineValues)[SupportSize] = _splineValues[jj][kk][ll]; + Real (*_wSplineValues)[SupportSize] = wSplineValues[jj][kk][ll]; + Real (*_dSplineValues)[SupportSize] = dSplineValues[jj][kk][ll]; + Real weight = pData.weight; + Point3D< Real > p = pData.position; +#endif // POINT_DATA_RES + + // evaluate the point p at all the nodes whose functions have it in their support + for( int s=-LeftPointSupportRadius ; s<=RightPointSupportRadius ; s++ ) for( int dd=0 ; dd::FunctionIndex( d , pOff[dd]+s ); + if( fIdx>=fStart && fIdxvalueWeight * weight; + Point3D< Real > weightedGradient; + if( HasGradients ) + { + Point3D< Real > gradient + ( + _dSplineValues[0][-j+LeftPointSupportRadius] * __splineValues[1][-k+LeftPointSupportRadius] * __splineValues[2][-l+LeftPointSupportRadius] , + __splineValues[0][-j+LeftPointSupportRadius] * _dSplineValues[1][-k+LeftPointSupportRadius] * __splineValues[2][-l+LeftPointSupportRadius] , + __splineValues[0][-j+LeftPointSupportRadius] * __splineValues[1][-k+LeftPointSupportRadius] * _dSplineValues[2][-l+LeftPointSupportRadius] + ); + weightedGradient = gradient * interpolationInfo->gradientWeight * weight; + diagonal += value * weightedValue + Point3D< Real >::Dot( gradient , weightedGradient ); + } + else diagonal += value * weightedValue; + + // Pre-multiply the x-coordinate values so that when we evaluate at one of the neighboring basis functions + // we get the product of the values of the center base function and the base function of the neighboring node + if( HasGradients ) for( int s=0 ; s& pData = *( (*interpolationInfo)( _node ) ); +#if POINT_DATA_RES + for( int c=0 ; c::SAMPLES ; c++ ) if( pData[c].weight ) +#endif // POINT_DATA_RES + { +#if POINT_DATA_RES + Real (*__splineValues)[SupportSize] = _splineValues[c][i+LeftSupportRadius][j+LeftSupportRadius][k+LeftSupportRadius]; + Real (*_wSplineValues)[SupportSize] = wSplineValues[c][i+LeftSupportRadius][j+LeftSupportRadius][k+LeftSupportRadius]; + Real (*_dSplineValues)[SupportSize] = dSplineValues[c][i+LeftSupportRadius][j+LeftSupportRadius][k+LeftSupportRadius]; +#else // !POINT_DATA_RES + Real (*__splineValues)[SupportSize] = _splineValues[i+LeftSupportRadius][j+LeftSupportRadius][k+LeftSupportRadius]; + Real (*_wSplineValues)[SupportSize] = wSplineValues[i+LeftSupportRadius][j+LeftSupportRadius][k+LeftSupportRadius]; + Real (*_dSplineValues)[SupportSize] = dSplineValues[i+LeftSupportRadius][j+LeftSupportRadius][k+LeftSupportRadius]; +#endif // POINT_DATA_RES + // Iterate over all neighbors whose support contains the point and accumulate the mutual integral + for( int ii=-LeftPointSupportRadius ; ii<=RightPointSupportRadius ; ii++ ) + for( int jj=-LeftPointSupportRadius ; jj<=RightPointSupportRadius ; jj++ ) + if( HasGradients ) + { + Real partialW_SplineValue = _wSplineValues[0][ii+LeftPointSupportRadius ] * __splineValues[1][jj+LeftPointSupportRadius ]; + Real partial__SplineValue = __splineValues[0][ii+LeftPointSupportRadius ] * __splineValues[1][jj+LeftPointSupportRadius ]; + Real partialD0SplineValue = _dSplineValues[0][ii+LeftPointSupportRadius ] * __splineValues[1][jj+LeftPointSupportRadius ]; + Real partialD1SplineValue = __splineValues[0][ii+LeftPointSupportRadius ] * _dSplineValues[1][jj+LeftPointSupportRadius ]; + Real* _pointValues = pointValues[i+ii+OverlapRadius][j+jj+OverlapRadius] + k + OverlapRadius; + Real* ___splineValues = __splineValues[2] + LeftPointSupportRadius; + Real* __dSplineValues = _dSplineValues[2] + LeftPointSupportRadius; + TreeOctNode* const * _neighbors = neighbors.neighbors[i+ii+OverlapRadius][j+jj+OverlapRadius] + k + OverlapRadius; + for( int kk=-LeftPointSupportRadius ; kk<=RightPointSupportRadius ; kk++ ) if( _isValidFEMNode( _neighbors[kk] ) ) + _pointValues[kk] += + partialW_SplineValue * ___splineValues[kk] + partialD0SplineValue * ___splineValues[kk] + partialD1SplineValue * ___splineValues[kk] + partial__SplineValue * __dSplineValues[kk]; + } + else + { + Real partialWSplineValue = _wSplineValues[0][ii+LeftPointSupportRadius ] * __splineValues[1][jj+LeftPointSupportRadius ]; + Real* _pointValues = pointValues[i+ii+OverlapRadius][j+jj+OverlapRadius] + k + OverlapRadius; + Real* ___splineValues = __splineValues[2] + LeftPointSupportRadius; + TreeOctNode* const * _neighbors = neighbors.neighbors[i+ii+OverlapRadius][j+jj+OverlapRadius] + k + OverlapRadius; + for( int kk=-LeftPointSupportRadius ; kk<=RightPointSupportRadius ; kk++ ) if( _isValidFEMNode( _neighbors[kk] ) ) + _pointValues[kk] += partialWSplineValue * ___splineValues[kk]; + } + } + } + } + } + pointValues[OverlapRadius][OverlapRadius][OverlapRadius] = diagonal; + int nodeIndex = neighbors.neighbors[OverlapRadius][OverlapRadius][OverlapRadius]->nodeData.nodeIndex; + if( isInterior ) // General case, so try to make fast + { + const TreeOctNode* const * _nodes = &neighbors.neighbors[0][0][0]; + const double* _stencil = &stencil( 0 , 0 , 0 ); + Real* _values = &pointValues[0][0][0]; + const static int CenterIndex = OverlapSize*OverlapSize*OverlapRadius + OverlapSize*OverlapRadius + OverlapRadius; + if( interpolationInfo ) for( int i=0 ; i( nodeIndex-offset , _values[CenterIndex] ); + for( int i=0 ; i( _nodes[i]->nodeData.nodeIndex-offset , _values[i] ); + } + else + { + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( node , d , off ); + Real temp = (Real)F.integrate( integrator , off , off ); + if( interpolationInfo ) temp += pointValues[OverlapRadius][OverlapRadius][OverlapRadius]; + row[count++] = MatrixEntry< Real >( nodeIndex-offset , temp ); + for( int x=0 ; x( _node->nodeData.nodeIndex-offset , temp ); + } + } + return count; +} + +template< class Real > +template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients > +int Octree< Real >::_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 ) +{ + static const int OverlapRadius = - BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + static const int OverlapSize = BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize; + + size_t start = _sNodesBegin(depth) , end = _sNodesEnd(depth) , range = end-start; + Stencil< double , OverlapSize > stencil , stencils[2][2][2]; + SystemCoefficients< FEMDegree , BType , FEMDegree , BType >::SetCentralSystemStencil ( F , integrator , stencil ); + SystemCoefficients< FEMDegree , BType , FEMDegree , BType >::SetCentralSystemStencils( F , childIntegrator , stencils ); + matrix.Resize( (int)range ); + std::vector< AdjacenctNodeKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i neighbors; + neighborKey.template getNeighbors< false , OverlapRadius , OverlapRadius >( node , neighbors ); + int count = _getMatrixRowSize< FEMDegree , BType >( neighbors ); + // Allocate memory for the row + matrix.SetRowSize( i , count ); + + // Set the row entries + matrix.rowSizes[i] = _setMatrixRow( F , interpolationInfo , neighbors , matrix[i] , (int)start , integrator , stencil , bsData ); + if( coarseToFine && depth>0 ) + { + // Offset the constraints using the solution from lower resolutions. + int x , y , z; + Cube::FactorCornerIndex( (int)( node - node->parent->children ) , x , y , z ); + typename TreeOctNode::Neighbors< OverlapSize > pNeighbors; + neighborKey.template getNeighbors< false , OverlapRadius , OverlapRadius >( node->parent , pNeighbors ); + _updateConstraintsFromCoarser( F , interpolationInfo , neighbors , pNeighbors , node , constraints , metSolution , childIntegrator , stencils[x][y][z] , bsData ); + } + } + memoryUsage(); + return 1; +} + +template< class Real > +template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients > +int Octree< Real >::_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 ) +{ + static const int OverlapSize = BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize; + static const int OverlapRadius = -BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + + int nStart = _sNodesBegin( depth , slice ) , nEnd = _sNodesEnd( depth , slice ); + size_t range = nEnd - nStart; + Stencil< double , OverlapSize > stencil , stencils[2][2][2]; + SystemCoefficients< FEMDegree , BType , FEMDegree , BType >::SetCentralSystemStencil ( F , integrator , stencil ); + SystemCoefficients< FEMDegree , BType , FEMDegree , BType >::SetCentralSystemStencils( F , childIntegrator , stencils ); + + matrix.Resize( (int)range ); + std::vector< AdjacenctNodeKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i neighbors; + neighborKey.template getNeighbors< false , OverlapRadius , OverlapRadius >( node , neighbors ); + int count = _getMatrixRowSize< FEMDegree , BType >( neighbors ); + + // Allocate memory for the row + matrix.SetRowSize( i , count ); + + // Set the row entries + matrix.rowSizes[i] = _setMatrixRow( F , interpolationInfo , neighbors , matrix[i] , _sNodesBegin( depth , slice ) , integrator , stencil , bsData ); + + if( coarseToFine && depth>0 ) + { + // Offset the constraints using the solution from lower resolutions. + int x , y , z; + Cube::FactorCornerIndex( (int)( node - node->parent->children ) , x , y , z ); + typename TreeOctNode::Neighbors< OverlapSize > pNeighbors; + neighborKey.template getNeighbors< false, OverlapRadius , OverlapRadius >( node->parent , pNeighbors ); + _updateConstraintsFromCoarser( F , interpolationInfo , neighbors , pNeighbors , node , constraints , metSolution , childIntegrator , stencils[x][y][z] , bsData ); + } + } +#if !defined( _WIN32 ) && !defined( _WIN64 ) +#pragma message( "[WARNING] I'm not sure how expensive this system call is on non-Windows system. (You may want to comment this out.)" ) +#endif // !_WIN32 && !_WIN64 + memoryUsage(); + return 1; +} + +#ifndef MOD +#define MOD( a , b ) ( (a)>0 ? (a) % (b) : ( (b) - ( -(a) % (b) ) ) % (b) ) +#endif // MOD +template< class Real > +template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients > +int Octree< Real >::_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 ) +{ + const int OverlapRadius = -BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + 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; + BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::SetIntegrator( integrator , depth ); + if( depth>0 ) BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::SetChildIntegrator( childIntegrator , depth-1 ); + + DenseNodeData< Real , FEMDegree >& metSolution = metSolutionConstraints; // This stores the up-sampled solution up to depth-2 + DenseNodeData< Real , FEMDegree >& metConstraints = metSolutionConstraints; // This stores the down-sampled constraints up to depth + + int sliceBegin = _BSplineBegin< FEMDegree , BType >( depth ) , sliceEnd = _BSplineEnd< FEMDegree , BType >( depth ); + double& systemTime = stats. systemTime; + double& solveTime = stats. solveTime; + double& evaluateTime = stats.evaluateTime; + systemTime = solveTime = evaluateTime = 0.; + + if( coarseToFine ) + { + if( depth>0 ) + { + // Up-sample the cumulative change in solution @(depth-2) into the cumulative change in solution @(depth-1) + if( depth-2>=0 ) _upSample< Real , FEMDegree , BType >( depth-1 , metSolution ); + // Add in the change in solution @(depth-1) +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodesBegin(depth-1) ; i<_sNodesEnd(depth-1) ; i++ ) metSolution[i] += solution[i]; + // evaluate the points @(depth) using the cumulative change in solution @(depth-1) + if( interpolationInfo ) + { + evaluateTime = Time(); + _setPointValuesFromCoarser( *interpolationInfo , depth , bsData , metSolution ); + evaluateTime = Time() - evaluateTime; + } + } + } + else if( depth<_maxDepth ) for( int i=_sNodesBegin(depth) ; i<_sNodesEnd(depth) ; i++ ) constraints[i] -= metConstraints[i]; + double bNorm = 0 , inRNorm = 0 , outRNorm = 0; + if( depth>=0 ) + { + // Add padding space if we are computing residuals + int frontOffset = computeNorms ? OverlapRadius : 0 , backOffset = computeNorms ? OverlapRadius : 0; + // Set the number of in-memory slices required for a temporally blocked solver + int solveSlices = std::max< int >( 0 , std::min< int >( OverlapRadius*iters - (OverlapRadius-1) , sliceEnd-sliceBegin ) ) , matrixSlices = std::max< int >( 1 , std::min< int >( solveSlices+frontOffset+backOffset , sliceEnd-sliceBegin ) ); + // The list of matrices for each in-memory slices + std::vector< SparseMatrix< Real > > _M( matrixSlices ); + // The list of multi-colored indices for each in-memory slice + std::vector< std::vector< std::vector< int > > > __mcIndices( solveSlices ); + + int dir = coarseToFine ? -1 : 1 , start = coarseToFine ? sliceEnd-1 : sliceBegin , end = coarseToFine ? sliceBegin-1 : sliceEnd; + for( int frontSlice=start-frontOffset*dir , backSlice = frontSlice-OverlapRadius*(iters-1)*dir ; backSlice!=end+backOffset*dir ; frontSlice+=dir , backSlice+=dir ) + { + double t; + if( frontSlice+frontOffset*dir>=sliceBegin && frontSlice+frontOffset*dir ) start = _M[_s][j]; + ConstPointer( MatrixEntry< Real > ) end = start + _M[_s].rowSizes[j]; + ConstPointer( MatrixEntry< Real > ) e; + for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value; + bNorm += B[j]*B[j]; + inRNorm += (temp-B[j]) * (temp-B[j]); + } + } + } + t = Time(); + // Compute the multicolor indices + if( iters && frontSlice>=sliceBegin && frontSlice( _sNodesBegin( depth , s ) , _sNodesEnd( depth , s ) , __mcIndices[__s] ); + } + // Advance through the in-memory slices, taking an appropriately sized stride + for( int slice=frontSlice ; slice*dir>=backSlice*dir ; slice-=OverlapRadius*dir ) + if( slice>=sliceBegin && slice::SolveGS( __mcIndices[__s] , _M[_s] , B , X , !coarseToFine , threads ); + } + solveTime += Time() - t; + // Compute residuals + if( computeNorms && backSlice-backOffset*dir>=sliceBegin && backSlice-backOffset*dir ) start = _M[_s][j]; + ConstPointer( MatrixEntry< Real > ) end = start + _M[_s].rowSizes[j]; + ConstPointer( MatrixEntry< Real > ) e; + for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value; + outRNorm += (temp-B[j]) * (temp-B[j]); + } + } + } + } + if( computeNorms ) stats.bNorm2 = bNorm , stats.inRNorm2 = inRNorm , stats.outRNorm2 = outRNorm; + + if( !coarseToFine && depth>0 ) + { + // Explicitly compute the restriction of the met solution onto the coarser nodes + // and down-sample the previous accumulation + { + _updateCumulativeIntegralConstraintsFromFiner( F , bsData , depth , solution , metConstraints ); + if( interpolationInfo ) _updateCumulativeInterpolationConstraintsFromFiner( *interpolationInfo , bsData , depth , solution , metConstraints ); + if( depth<_maxDepth ) _downSample< Real , FEMDegree , BType >( depth , metConstraints ); + } + } + memoryUsage(); + + return iters; +} +#undef MOD + +template< class Real > +template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients > +int Octree< Real >::_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 ) +{ + 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; + BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::SetIntegrator( integrator , depth ); + if( depth>0 ) BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::SetChildIntegrator( childIntegrator , depth-1 ); + + DenseNodeData< Real , FEMDegree >& metSolution = metSolutionConstraints; // This stores the up-sampled solution up to depth-2 + DenseNodeData< Real , FEMDegree >& metConstraints = metSolutionConstraints; // This stores the down-sampled constraints up to depth + + int iter = 0; + Pointer( Real ) X = GetPointer( & solution[0] + _sNodesBegin(depth) , _sNodesSize(depth) ); + Pointer( Real ) B = GetPointer( &constraints[0] + _sNodesBegin(depth) , _sNodesSize(depth) ); + SparseMatrix< Real > M; + double& systemTime = stats. systemTime; + double& solveTime = stats. solveTime; + double& evaluateTime = stats.evaluateTime; + systemTime = solveTime = evaluateTime = 0.; + + if( coarseToFine ) + { + if( depth>0 ) + { + // Up-sample the cumulative change in solution @(depth-2) into the cumulative change in solution @(depth-1) + if( depth-2>=0 ) _upSample< Real , FEMDegree , BType >( depth-1 , metSolution ); + // Add in the change in solution @(depth-1) +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodesBegin(depth-1) ; i<_sNodesEnd(depth-1) ; i++ ) metSolution[i] += solution[i]; + // evaluate the points @(depth) using the cumulative change in solution @(depth-1) + if( interpolationInfo ) + { + evaluateTime = Time(); + _setPointValuesFromCoarser( *interpolationInfo , depth , bsData , metSolution ); + evaluateTime = Time() - evaluateTime; + } + } + } + else if( depth<_maxDepth ) for( int i=_sNodesBegin(depth) ; i<_sNodesEnd(depth) ; i++ ) constraints[i] -= metConstraints[i]; + + // Get the system matrix (and adjust the right-hand-side based on the coarser solution if prolonging) + systemTime = Time(); + _getMatrixAndUpdateConstraints( F , interpolationInfo , M , constraints , integrator , childIntegrator , bsData , depth , metSolution , coarseToFine ); + systemTime = Time()-systemTime; + + solveTime = Time(); + // Solve the linear system + accuracy = Real( accuracy / 100000 ) * M.rows; + int dim = _BSplineEnd< FEMDegree , BType >( depth ) - _BSplineBegin< FEMDegree , BType >( depth ); + int nonZeroRows = 0; + for( int i=0 ; ivalueWeight ) && HasPartitionOfUnity< BType >() && F.vanishesOnConstants() ); + double bNorm = 0 , inRNorm = 0 , outRNorm = 0; + if( computeNorms ) + { +#pragma omp parallel for num_threads( threads ) reduction( + : bNorm , inRNorm ) + for( int j=0 ; j ) start = M[j]; + ConstPointer( MatrixEntry< Real > ) end = start + M.rowSizes[j]; + ConstPointer( MatrixEntry< Real > ) e; + for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value; + bNorm += B[j] * B[j]; + inRNorm += ( temp-B[j] ) * ( temp-B[j] ); + } + } + + iters = std::min< int >( nonZeroRows , iters ); + if( iters ) iter += SparseMatrix< Real >::SolveCG( M , ( ConstPointer( Real ) )B , iters , X , Real( accuracy ) , 0 , addDCTerm , false , threads ); + + solveTime = Time()-solveTime; + if( computeNorms ) + { +#pragma omp parallel for num_threads( threads ) reduction( + : outRNorm ) + for( int j=0 ; j ) start = M[j]; + ConstPointer( MatrixEntry< Real > ) end = start + M.rowSizes[j]; + ConstPointer( MatrixEntry< Real > ) e; + for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value; + outRNorm += ( temp-B[j] ) * ( temp-B[j] ); + } + stats.bNorm2 = bNorm , stats.inRNorm2 = inRNorm , stats.outRNorm2 = outRNorm; + } + + // Copy the old solution into the buffer, write in the new solution, compute the change, and update the met solution + if( !coarseToFine && depth>0 ) + { + // Explicitly compute the restriction of the met solution onto the coarser nodes + // and down-sample the previous accumulation + { + _updateCumulativeIntegralConstraintsFromFiner( F , bsData , depth , solution , metConstraints ); + if( interpolationInfo ) _updateCumulativeInterpolationConstraintsFromFiner( *interpolationInfo , bsData , depth , solution , metConstraints ); + if( depth>_maxDepth ) _downSample< Real , FEMDegree , BType >( depth , metConstraints ); + } + } + memoryUsage(); + return iter; +} + +template< class Real > +template< int FEMDegree , BoundaryType BType > +int Octree< Real >::_getMatrixRowSize( const typename TreeOctNode::Neighbors< BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize >& neighbors ) const +{ + static const int OverlapSize = BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize; + static const int OverlapRadius = - BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + + int count = 0; + int nodeIndex = neighbors.neighbors[OverlapRadius][OverlapRadius][OverlapRadius]->nodeData.nodeIndex; + const TreeOctNode* const * _nodes = &neighbors.neighbors[0][0][0]; + for( int i=0 ; i +template< int FEMDegree1 , int FEMDegree2 > +void Octree< Real >::_SetParentOverlapBounds( const TreeOctNode* node , int& startX , int& endX , int& startY , int& endY , int& startZ , int& endZ ) +{ + const int OverlapStart = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::OverlapStart; + + if( node->parent ) + { + int x , y , z , c = (int)( node - node->parent->children ); + Cube::FactorCornerIndex( c , x , y , z ); + startX = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::ParentOverlapStart[x]-OverlapStart , endX = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::ParentOverlapEnd[x]-OverlapStart+1; + startY = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::ParentOverlapStart[y]-OverlapStart , endY = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::ParentOverlapEnd[y]-OverlapStart+1; + startZ = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::ParentOverlapStart[z]-OverlapStart , endZ = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::ParentOverlapEnd[z]-OverlapStart+1; + } +} + +// It is assumed that at this point, the evaluationg of the current depth's points, using the coarser resolution solution +// has already happened +template< class Real > +template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients > +void Octree< Real >::_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 >& lapStencil , const BSplineData< FEMDegree , BType >& bsData ) const +{ + static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int OverlapRadius = - BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + + if( _localDepth( node )<=0 ) return; + // This is a conservative estimate as we only need to make sure that the parent nodes don't overlap the child (not the parent itself) + bool isInterior = _isInteriorlyOverlapped< FEMDegree , FEMDegree >( node->parent ); + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( node , d , off ); + + // Offset the constraints using the solution from lower resolutions. + int startX , endX , startY , endY , startZ , endZ; + _SetParentOverlapBounds< FEMDegree , FEMDegree >( node , startX , endX , startY , endY , startZ , endZ ); + + for( int x=startX ; xnodeData.nodeIndex ]; + { + if( isInterior ) constraints[ node->nodeData.nodeIndex ] -= Real( lapStencil( x , y , z ) * _solution ); + else + { + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( _node , _d , _off ); + constraints[ node->nodeData.nodeIndex ] -= (Real)F.integrate( childIntegrator , _off , off ) * _solution; + } + } + } + + if( interpolationInfo ) + { + double constraint = 0; + int fIdx[3]; + functionIndex< FEMDegree , BType >( node , fIdx ); + // evaluate the current node's basis function at adjacent points + for( int x=-LeftSupportRadius ; x<=RightSupportRadius ; x++ ) for( int y=-LeftSupportRadius ; y<=RightSupportRadius ; y++ ) for( int z=-LeftSupportRadius ; z<=RightSupportRadius ; z++ ) + { + const TreeOctNode* _node = neighbors.neighbors[x+OverlapRadius][y+OverlapRadius][z+OverlapRadius]; + if( _isValidSpaceNode( _node ) && (*interpolationInfo)( _node ) ) + { + const PointData< Real , HasGradients >& pData = *( (*interpolationInfo)( _node ) ); + constraint += _ConstraintCalculator_< Real , FEMDegree , HasGradients >::_CalculateConstraint_ + ( + pData , + bsData. baseBSplines[ fIdx[0] ][x+LeftSupportRadius] , + bsData. baseBSplines[ fIdx[1] ][y+LeftSupportRadius] , + bsData. baseBSplines[ fIdx[2] ][z+LeftSupportRadius] , + bsData.dBaseBSplines[ fIdx[0] ][x+LeftSupportRadius] , + bsData.dBaseBSplines[ fIdx[1] ][y+LeftSupportRadius] , + bsData.dBaseBSplines[ fIdx[2] ][z+LeftSupportRadius] + ); + } + } + constraints[ node->nodeData.nodeIndex ] -= Real( constraint ); + } +} + +// Given the solution @( depth ) add to the met constraints @( depth-1 ) +template< class Real > +template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor > +void Octree< Real >::_updateCumulativeIntegralConstraintsFromFiner( const FEMSystemFunctor& F , const BSplineData< FEMDegree , BType >& bsData , LocalDepth highDepth , const DenseNodeData< Real , FEMDegree >& fineSolution , DenseNodeData< Real , FEMDegree >& coarseConstraints ) const +{ + typename BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::FunctionIntegrator::template ChildIntegrator< DERIVATIVES( FEMDegree ) , DERIVATIVES( FEMDegree ) > childIntegrator; + BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::SetChildIntegrator( childIntegrator , highDepth-1 ); + + static const int OverlapSize = BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize; + static const int OverlapRadius = - BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + typedef typename TreeOctNode::NeighborKey< -BSplineSupportSizes< FEMDegree >::SupportStart , BSplineSupportSizes< FEMDegree >::SupportEnd >SupportKey; + + if( highDepth<=0 ) return; + // Get the stencil describing the Laplacian relating coefficients @(depth) with coefficients @(depth-1) + Stencil< double , OverlapSize > stencils[2][2][2]; + SystemCoefficients< FEMDegree , BType , FEMDegree , BType >::SetCentralSystemStencils( F , childIntegrator , stencils ); + size_t start = _sNodesBegin( highDepth) , end = _sNodesEnd(highDepth) , range = end-start; + int lStart = _sNodesBegin(highDepth-1); + + // Iterate over the nodes @( depth ) + std::vector< SupportKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; iparent->children ); + Cube::FactorCornerIndex( c , x , y , z ); + { + typename TreeOctNode::Neighbors< OverlapSize > pNeighbors; + neighborKey.template getNeighbors< false , OverlapRadius , OverlapRadius >( node->parent , pNeighbors ); + const Stencil< double , OverlapSize >& stencil = stencils[x][y][z]; + + bool isInterior = _isInteriorlyOverlapped< FEMDegree , FEMDegree >( node->parent ); + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( node , d , off ); + + // Offset the constraints using the solution from finer resolutions. + int startX , endX , startY , endY , startZ , endZ; + _SetParentOverlapBounds< FEMDegree , FEMDegree >( node , startX , endX , startY , endY , startZ , endZ ); + + Real solution = fineSolution[ node->nodeData.nodeIndex ]; + for( int x=startX ; xnodeData.nodeIndex ] += Real( stencil( x , y , z ) * solution ); + else + { + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( _node , _d , _off ); +#pragma omp atomic + coarseConstraints[ _node->nodeData.nodeIndex ] += Real( F.integrate( childIntegrator , _off , off ) * solution ); + } + } + } + } +} + +template< class Real > +template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients > +void Octree< Real >::setSystemMatrix( const FEMSystemFunctor& F , const InterpolationInfo< HasGradients >* interpolationInfo , LocalDepth depth , SparseMatrix< Real >& matrix ) const +{ + if( depth<0 || depth>_maxDepth ) fprintf( stderr , "[ERROR] System depth out of bounds: %d <= %d <= %d\n" , 0 , depth , _maxDepth ) , exit( 0 ); + typename BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::FunctionIntegrator::template Integrator< DERIVATIVES( FEMDegree ) , DERIVATIVES( FEMDegree ) > integrator; + BSplineIntegrationData< FEMDegree , BType , FEMDegree , BType >::SetIntegrator( integrator , depth ); + BSplineData< FEMDegree , BType > bsData( depth ); + + static const int OverlapRadius = - BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + static const int OverlapSize = BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize; + + Stencil< double , OverlapSize > stencil; + SystemCoefficients< FEMDegree , BType , FEMDegree , BType >::SetCentralSystemStencil ( F , integrator , stencil ); + + matrix.Resize( _sNodesSize(depth) ); + std::vector< AdjacenctNodeKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i neighbors; + neighborKey.template getNeighbors< false , OverlapRadius , OverlapRadius >( _sNodes.treeNodes[i] , neighbors ); + + matrix.SetRowSize( ii , _getMatrixRowSize< FEMDegree , BType >( neighbors ) ); + matrix.rowSizes[ii] = _setMatrixRow( F , interpolationInfo , neighbors , matrix[ii] , _sNodesBegin(depth) , integrator , stencil , bsData ); + } +} + +template< class Real > +template< int FEMDegree , BoundaryType BType , class FEMSystemFunctor , bool HasGradients > +DenseNodeData< Real , FEMDegree > Octree< Real >::solveSystem( const FEMSystemFunctor& F , InterpolationInfo< HasGradients >* interpolationInfo , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxSolveDepth , const typename Octree< Real >::SolverInfo& solverInfo ) +{ + BSplineData< FEMDegree , BType > bsData( maxSolveDepth ); + + maxSolveDepth = std::min< LocalDepth >( maxSolveDepth , _maxDepth ); + int iter = 0; + const int _iters = std::max< int >( 0 , solverInfo.iters ); + + DenseNodeData< Real , FEMDegree > solution( _sNodesEnd( _maxDepth ) ); + memset( &solution[0] , 0 , sizeof(Real) * _sNodesEnd( _maxDepth ) ); + + DenseNodeData< Real , FEMDegree > metSolution( _sNodesEnd( _maxDepth-1 ) ); + memset( &metSolution[0] , 0 , sizeof(Real)*_sNodesEnd( _maxDepth-1 ) ); + for( LocalDepth d=0 ; d<=maxSolveDepth ; d++ ) + { + int iters = (int)ceil( _iters * pow( solverInfo.lowResIterMultiplier , maxSolveDepth-d ) ); + _SolverStats sStats; + if( !d ) iter = _solveSystemCG( F , bsData , interpolationInfo , d , solution , constraints , metSolution , _sNodesSize(d) , true , sStats , solverInfo.showResidual , 0 ); + else + { + if( d>solverInfo.cgDepth ) iter = _solveSystemGS( F , bsData , interpolationInfo , d , solution , constraints , metSolution , iters , true , sStats , solverInfo.showResidual ); + else iter = _solveSystemCG( F , bsData , interpolationInfo , d , solution , constraints , metSolution , iters , true , sStats , solverInfo.showResidual , solverInfo.cgAccuracy ); + } + int femNodes = 0; +#pragma omp parallel for reduction( + : femNodes ) + for( int i=_sNodesBegin(d) ; i<_sNodesEnd(d) ; i++ ) if( _isValidFEMNode( _sNodes.treeNodes[i] ) ) femNodes++; + if( solverInfo.verbose ) + { + if( maxSolveDepth<10 ) printf( "Depth[%d/%d]:\t" , d , maxSolveDepth ); + else printf( "Depth[%2d/%d]:\t" , d , maxSolveDepth ); + printf( "Evaluated / Got / Solved in: %6.3f / %6.3f / %6.3f\t(%.3f MB)\tNodes: %d\n" , sStats.evaluateTime , sStats.systemTime , sStats.solveTime , _localMemoryUsage , femNodes ); + } + if( solverInfo.showResidual && iters ) + { + for( LocalDepth dd=0 ; dd %.4e -> %.4e (%.2e) [%d]\n" , d<=solverInfo.cgDepth ? "CG" : "GS" , sqrt( sStats.bNorm2 ) , sqrt( sStats.inRNorm2 ) , sqrt( sStats.outRNorm2 ) , sqrt( sStats.outRNorm2 / sStats.bNorm2 ) , iters ); + } + } + memoryUsage(); + return solution; +} + +template< class Real > +template< int FEMDegree > +DenseNodeData< Real , FEMDegree > Octree< Real >::initDenseNodeData( void ) +{ + DenseNodeData< Real , FEMDegree > constraints( _sNodes.size() ); + memset( &constraints[0] , 0 , sizeof(Real)*_sNodes.size() ); + return constraints; +} +template< > template< > float Octree< float >::_Dot( const float & r1 , const float & r2 ){ return r1*r2; } +template< > template< > double Octree< double >::_Dot( const double& r1 , const double& r2 ){ return r1*r2; } +template< > template< > float Octree< float >::_Dot( const Point3D< float >& p1 , const Point3D< float >& p2 ){ return Point3D< float >::Dot( p1 , p2 ); } +template< > template< > double Octree< double >::_Dot( const Point3D< double >& p1 , const Point3D< double >& p2 ){ return Point3D< double >::Dot( p1 , p2 ); } +template< > template< > bool Octree< float >::_IsZero( const float & r ){ return r==0; } +template< > template< > bool Octree< double >::_IsZero( const double& r ){ return r==0; } +template< > template< > bool Octree< float >::_IsZero( const Point3D< float >& p ){ return p[0]==0 && p[1]==0 && p[2]==0; } +template< > template< > bool Octree< double >::_IsZero( const Point3D< double >& p ){ return p[0]==0 && p[1]==0 && p[2]==0; } +template< class Real > +template< int FEMDegree , BoundaryType FEMBType , int CDegree , BoundaryType CBType , class FEMConstraintFunctor , class Coefficients , class D , class _D > +void Octree< Real >::_addFEMConstraints( const FEMConstraintFunctor& F , const Coefficients& coefficients , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxDepth ) +{ + typedef typename TreeOctNode::NeighborKey< -BSplineSupportSizes< FEMDegree >::SupportStart , BSplineSupportSizes< FEMDegree >::SupportEnd > SupportKey; + const int CFEMOverlapSize = BSplineOverlapSizes< CDegree , FEMDegree >::OverlapSize; + const int LeftCFEMOverlapRadius = -BSplineOverlapSizes< CDegree , FEMDegree >::OverlapStart; + const int RightCFEMOverlapRadius = BSplineOverlapSizes< CDegree , FEMDegree >::OverlapEnd; + const int LeftFEMCOverlapRadius = -BSplineOverlapSizes< FEMDegree , CDegree >::OverlapStart; + const int RightFEMCOverlapRadius = BSplineOverlapSizes< FEMDegree , CDegree >::OverlapEnd; + + // To set the constraints, we iterate over the + // splatted normals and compute the dot-product of the + // divergence of the normal field with all the basis functions. + // Within the same depth: set directly as a gather + // Coarser depths + maxDepth = std::min< LocalDepth >( maxDepth , _maxDepth ); + DenseNodeData< Real , FEMDegree >* __constraints = new DenseNodeData< Real , FEMDegree >( _sNodesEnd(maxDepth-1) ); + DenseNodeData< Real , FEMDegree >& _constraints = *__constraints; + memset( &_constraints[0] , 0 , sizeof(Real)*( _sNodesEnd(maxDepth-1) ) ); + memoryUsage(); + + for( LocalDepth d=maxDepth ; d>=0 ; d-- ) + { + Stencil< _D , CFEMOverlapSize > stencil , stencils[2][2][2]; + typename SystemCoefficients< CDegree , CBType , FEMDegree , FEMBType >:: Integrator integrator; + typename SystemCoefficients< FEMDegree , FEMBType , CDegree , CBType >::ChildIntegrator childIntegrator; + BSplineIntegrationData< CDegree , CBType , FEMDegree , FEMBType >::SetIntegrator( integrator , d ); + if( d>0 ) BSplineIntegrationData< FEMDegree , FEMBType , CDegree , CBType >::SetChildIntegrator( childIntegrator , d-1 ); + SystemCoefficients< CDegree , CBType , FEMDegree , FEMBType >::template SetCentralConstraintStencil < false >( F, integrator , stencil ); + SystemCoefficients< FEMDegree , FEMBType , CDegree , CBType >::template SetCentralConstraintStencils< true >( F, childIntegrator , stencils ); + + std::vector< SupportKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i neighbors; + neighborKey.template getNeighbors< false , LeftFEMCOverlapRadius , RightFEMCOverlapRadius >( node , neighbors ); + bool isInterior = _isInteriorlyOverlapped< FEMDegree , CDegree >( node ) , isInterior2 = _isInteriorlyOverlapped< CDegree , FEMDegree >( node->parent ); + + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( node , d , off ); + // Set constraints from current depth + // Gather the constraints from the vector-field at _node into the constraint stored with node + if( _isValidFEMNode( node ) ) + { + for( int x=startX ; x( _node ) ) + { + const D* d = coefficients( _node ); + if( d ) + if( isInterior ) constraints[i] += _Dot( (D)stencil( x , y , z ) , *d ); + else + { + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( _node , _d , _off ); + constraints[i] += _Dot( *d , (D)F.template integrate< false >( integrator , _off , off ) ); + } + } + } + _SetParentOverlapBounds< CDegree , FEMDegree >( node , startX , endX , startY , endY , startZ , endZ ); + } + if( !isValidFEMNode< CDegree , CBType >( node ) ) continue; + const D* _data = coefficients( node ); + if( !_data ) continue; + const D& data = *_data; + if( _IsZero( data ) ) continue; + + // Set the _constraints for the parents + if( d>0 ) + { + int cx , cy , cz; + Cube::FactorCornerIndex( (int)( node - node->parent->children ) , cx , cy ,cz ); + const Stencil< _D , CFEMOverlapSize >& _stencil = stencils[cx][cy][cz]; + + neighborKey.template getNeighbors< false , LeftCFEMOverlapRadius , RightCFEMOverlapRadius >( node->parent , neighbors ); + + for( int x=startX ; x( childIntegrator , _off , off ) ); + } +#pragma omp atomic + _constraints[ _node->nodeData.nodeIndex ] += c; + } + } + } + } + memoryUsage(); + } + + // Fine-to-coarse down-sampling of constraints + for( LocalDepth d=maxDepth-1 ; d>0 ; d-- ) _downSample< Real , FEMDegree , FEMBType >( d , _constraints ); + + // Add the accumulated constraints from all finer depths +#pragma omp parallel for num_threads( threads ) + for( int i=0 ; i<_sNodesEnd(maxDepth-1) ; i++ ) constraints[i] += _constraints[i]; + + delete __constraints; + + DenseNodeData< D , CDegree > _coefficients( _sNodesEnd(maxDepth-1) ); + memset( &_coefficients[0] , 0 , sizeof(D) * _sNodesEnd(maxDepth-1) ); + for( LocalDepth d=maxDepth-1 ; d>=0 ; d-- ) + { +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodesBegin(d) ; i<_sNodesEnd(d) ; i++ ) if( isValidFEMNode< CDegree , CBType >( _sNodes.treeNodes[i] ) ) + { + const D* d = coefficients( _sNodes.treeNodes[i] ); + if( d ) _coefficients[i] += *d; + } + } + + // Coarse-to-fine up-sampling of coefficients + for( LocalDepth d=1 ; d( d , _coefficients ); + + // Compute the contribution from all coarser depths + for( LocalDepth d=1 ; d<=maxDepth ; d++ ) + { + size_t start = _sNodesBegin( d ) , end = _sNodesEnd( d ) , range = end - start; + Stencil< _D , CFEMOverlapSize > stencils[2][2][2]; + typename SystemCoefficients< CDegree , CBType , FEMDegree , FEMBType >::ChildIntegrator childIntegrator; + BSplineIntegrationData< CDegree , CBType , FEMDegree , FEMBType >::SetChildIntegrator( childIntegrator , d-1 ); + SystemCoefficients< CDegree , CBType , FEMDegree , FEMBType >::template SetCentralConstraintStencils< false >( F , childIntegrator , stencils ); + std::vector< SupportKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i( node , startX , endX , startY , endY , startZ , endZ ); + typename TreeOctNode::Neighbors< CFEMOverlapSize > pNeighbors; + neighborKey.template getNeighbors< false , LeftFEMCOverlapRadius , RightFEMCOverlapRadius >( node->parent , pNeighbors ); + + bool isInterior = _isInteriorlyOverlapped< FEMDegree , CDegree >( node->parent ); + int cx , cy , cz; + if( d>0 ) + { + int c = int( node - node->parent->children ); + Cube::FactorCornerIndex( c , cx , cy , cz ); + } + else cx = cy = cz = 0; + Stencil< _D , CFEMOverlapSize >& _stencil = stencils[cx][cy][cz]; + + Real constraint = Real(0); + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( node , d , off ); + for( int x=startX ; x( _node ) ) + { + if( isInterior ) constraint += _Dot( _coefficients[ _node->nodeData.nodeIndex ] , (D)_stencil( x , y , z ) ); + else + { + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset ( _node , _d , _off ); + constraint += _Dot( _coefficients[ _node->nodeData.nodeIndex ] , (D)F.template integrate< false >( childIntegrator , _off , off ) ); + } + } + } + constraints[i] += constraint; + } + } + memoryUsage(); +} + +template< class Real > +template< int FEMDegree , BoundaryType BType , bool HasGradients > +void Octree< Real >::addInterpolationConstraints( const InterpolationInfo< HasGradients >& interpolationInfo , DenseNodeData< Real , FEMDegree >& constraints , LocalDepth maxDepth ) +{ + typedef typename TreeOctNode::NeighborKey< -BSplineSupportSizes< FEMDegree >::SupportStart , BSplineSupportSizes< FEMDegree >::SupportEnd > SupportKey; + maxDepth = std::min< LocalDepth >( maxDepth , _maxDepth ); + { + static const int OverlapSize = BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapSize; + static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart; + static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd; + static const int OverlapRadius = - BSplineOverlapSizes< FEMDegree , FEMDegree >::OverlapStart; + BSplineData< FEMDegree , BType > bsData( _maxDepth ); + for( int d=0 ; d<=maxDepth ; d++ ) + { + std::vector< SupportKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i neighbors; + neighborKey.template getNeighbors< false , OverlapRadius , OverlapRadius >( node , neighbors ); + + double constraint = 0; + int fIdx[3]; + functionIndex< FEMDegree , BType >( node , fIdx ); + // evaluate the current node's basis function at adjacent points + for( int x=-LeftSupportRadius ; x<=RightSupportRadius ; x++ ) for( int y=-LeftSupportRadius ; y<=RightSupportRadius ; y++ ) for( int z=-LeftSupportRadius ; z<=RightSupportRadius ; z++ ) + { + const TreeOctNode* _node = neighbors.neighbors[x+OverlapRadius][y+OverlapRadius][z+OverlapRadius]; + if( _isValidSpaceNode( _node ) && interpolationInfo( _node ) ) + { + const PointData< Real , HasGradients >& pData = *( interpolationInfo( _node ) ); + constraint += _ConstraintCalculator_< Real , FEMDegree , HasGradients >::_CalculateConstraint_ + ( + pData , + bsData. baseBSplines[ fIdx[0] ][x+LeftSupportRadius] , + bsData. baseBSplines[ fIdx[1] ][y+LeftSupportRadius] , + bsData. baseBSplines[ fIdx[2] ][z+LeftSupportRadius] , + bsData.dBaseBSplines[ fIdx[0] ][x+LeftSupportRadius] , + bsData.dBaseBSplines[ fIdx[1] ][y+LeftSupportRadius] , + bsData.dBaseBSplines[ fIdx[2] ][z+LeftSupportRadius] , + interpolationInfo.valueWeight , interpolationInfo.gradientWeight + ); + } + } + constraints[ node->nodeData.nodeIndex ] += (Real)constraint; + } + } + memoryUsage(); + } +} +template< class Real > +template< int FEMDegree1 , BoundaryType FEMBType1 , int FEMDegree2 , BoundaryType FEMBType2 , class DotFunctor , bool HasGradients , class Coefficients1 , class Coefficients2 > +double Octree< Real >::_dot( const DotFunctor& F , const InterpolationInfo< HasGradients >* iInfo , const Coefficients1& coefficients1 , const Coefficients2& coefficients2 ) const +{ + double dot = 0; + + // Calculate the contribution from @(depth,depth) + { + typedef typename TreeOctNode::ConstNeighborKey< -BSplineSupportSizes< FEMDegree1 >::SupportStart , BSplineSupportSizes< FEMDegree1 >::SupportEnd > SupportKey; + const int OverlapSize = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::OverlapSize; + const int LeftOverlapRadius = -BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::OverlapStart; + const int RightOverlapRadius = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::OverlapEnd; + + for( LocalDepth d=0 ; d<=_maxDepth ; d++ ) + { + Stencil< double , OverlapSize > stencil; + typename SystemCoefficients< FEMDegree1 , FEMBType1 , FEMDegree2 , FEMBType2 >::Integrator integrator; + BSplineIntegrationData< FEMDegree1 , FEMBType1 , FEMDegree2 , FEMBType2 >::SetIntegrator( integrator , d ); + SystemCoefficients< FEMDegree1 , FEMBType1 , FEMDegree2 , FEMBType2 >::template SetCentralConstraintStencil< false , DotFunctor >( F , integrator , stencil ); + + std::vector< SupportKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i( node ) && ( _data1=coefficients1(node) ) ) + { + SupportKey& neighborKey = neighborKeys[ omp_get_thread_num() ]; + typename TreeOctNode::ConstNeighbors< OverlapSize > neighbors; + neighborKey.template getNeighbors< LeftOverlapRadius , RightOverlapRadius >( node , neighbors ); + bool isInterior = _isInteriorlyOverlapped< FEMDegree1 , FEMDegree2 >( node ); + + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( node , d , off ); + + for( int x=0 ; x( _node ) && ( _data2=coefficients2( _node ) ) ) + if( isInterior ) dot += (*_data1) * (*_data2 ) * stencil( x , y , z ); + else + { + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( _node , _d , _off ); + dot += (*_data1) * (*_data2) * F.template integrate< false >( integrator , off , _off ); + } + } + } + } + } + } + // Calculate the contribution from @(::SupportStart , BSplineSupportSizes< FEMDegree1 >::SupportEnd > SupportKey; + const int OverlapSize = BSplineOverlapSizes< FEMDegree2 , FEMDegree1 >::OverlapSize; + const int LeftOverlapRadius = -BSplineOverlapSizes< FEMDegree2 , FEMDegree1 >::OverlapStart; + const int RightOverlapRadius = BSplineOverlapSizes< FEMDegree2 , FEMDegree1 >::OverlapEnd; + + DenseNodeData< Real , FEMDegree1 > cumulative1( _sNodesEnd( _maxDepth-1 ) ); + if( _maxDepth>0 ) memset( &cumulative1[0] , 0 , sizeof(Real) * _sNodesEnd( _maxDepth-1 ) ); + + for( LocalDepth d=1 ; d<=_maxDepth ; d++ ) + { + // Update the cumulative coefficients with the coefficients @(depth-1) +#pragma omp parallel for + for( int i=_sNodesBegin(d-1) ; i<_sNodesEnd(d-1) ; i++ ) + { + const Real* _data1 = coefficients1( _sNodes.treeNodes[i] ); + if( _data1 ) cumulative1[i] += *_data1; + } + + Stencil< double , OverlapSize > stencils[2][2][2]; + typename SystemCoefficients< FEMDegree1 , FEMBType1 , FEMDegree2 , FEMBType2 >::ChildIntegrator childIntegrator; + BSplineIntegrationData< FEMDegree1 , FEMBType1 , FEMDegree2 , FEMBType2 >::SetChildIntegrator( childIntegrator , d-1 ); + SystemCoefficients< FEMDegree1 , FEMBType1 , FEMDegree2 , FEMBType2 >::template SetCentralConstraintStencils< false >( F, childIntegrator , stencils ); + + std::vector< SupportKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i( node ) && ( _data2=coefficients2( node ) ) ) + { + SupportKey& neighborKey = neighborKeys[ omp_get_thread_num() ]; + bool isInterior = _isInteriorlyOverlapped< FEMDegree1 , FEMDegree2 >( node->parent ); + + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( node , d , off ); + + int cx , cy , cz; + Cube::FactorCornerIndex( (int)( node - node->parent->children ) , cx , cy ,cz ); + const Stencil< double , OverlapSize >& _stencil = stencils[cx][cy][cz]; + typename TreeOctNode::ConstNeighbors< OverlapSize > neighbors; + neighborKey.template getNeighbors< LeftOverlapRadius , RightOverlapRadius >( node->parent , neighbors ); + + int startX , endX , startY , endY , startZ , endZ; + _SetParentOverlapBounds< FEMDegree2 , FEMDegree1 >( node , startX , endX , startY , endY , startZ , endZ ); + for( int x=startX ; x( _node ) && ( _data1=cumulative1(_node) ) ) + { + if( isInterior ) dot += (*_data1) * (*_data2) * _stencil( x , y , z ); + else + { + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( _node , _d , _off ); + dot += (*_data1) * (*_data2) * F.template integrate< false >( childIntegrator , _off , off ); + } + } + } + } + } + // Up sample the cumulative coefficients for the next level + if( d<_maxDepth ) _upSample< Real , FEMDegree1 , FEMBType1 >( d , cumulative1 ); + } + } + + // Calculate the contribution from @(>depth,depth) + { + typedef typename TreeOctNode::ConstNeighborKey< -BSplineSupportSizes< FEMDegree2 >::SupportStart , BSplineSupportSizes< FEMDegree2 >::SupportEnd > SupportKey; + const int OverlapSize = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::OverlapSize; + const int LeftOverlapRadius = -BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::OverlapStart; + const int RightOverlapRadius = BSplineOverlapSizes< FEMDegree1 , FEMDegree2 >::OverlapEnd; + + DenseNodeData< Real , FEMDegree2 > cumulative2( _sNodesEnd( _maxDepth-1 ) ); + if( _maxDepth>0 ) memset( &cumulative2[0] , 0 , sizeof(Real) * _sNodesEnd( _maxDepth-1 ) ); + + for( LocalDepth d=_maxDepth ; d>0 ; d-- ) + { + Stencil< double , OverlapSize > stencils[2][2][2]; + typename SystemCoefficients< FEMDegree2 , FEMBType2 , FEMDegree1 , FEMBType1 >::ChildIntegrator childIntegrator; + BSplineIntegrationData< FEMDegree2 , FEMBType2 , FEMDegree1 , FEMBType1 >::SetChildIntegrator( childIntegrator , d-1 ); + SystemCoefficients< FEMDegree2 , FEMBType2 , FEMDegree1 , FEMBType1 >::template SetCentralConstraintStencils< true >( F , childIntegrator , stencils ); + + std::vector< SupportKey > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i( node ) && ( _data1=coefficients1( node ) ) ) + { + SupportKey& neighborKey = neighborKeys[ omp_get_thread_num() ]; + bool isInterior = _isInteriorlyOverlapped< FEMDegree2 , FEMDegree1 >( node->parent ); + + LocalDepth d ; LocalOffset off; + _localDepthAndOffset( node , d , off ); + + int cx , cy , cz; + Cube::FactorCornerIndex( (int)( node - node->parent->children ) , cx , cy ,cz ); + const Stencil< double , OverlapSize >& _stencil = stencils[cx][cy][cz]; + typename TreeOctNode::ConstNeighbors< OverlapSize > neighbors; + neighborKey.template getNeighbors< LeftOverlapRadius , RightOverlapRadius >( node->parent , neighbors ); + + int startX , endX , startY , endY , startZ , endZ; + _SetParentOverlapBounds< FEMDegree1 , FEMDegree2 >( node , startX , endX , startY , endY , startZ , endZ ); + + for( int x=startX ; x( _node ) ) + { + Real _dot; + if( isInterior ) _dot = (*_data1) * _stencil( x , y , z ); + else + { + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( _node , _d , _off ); + _dot = (*_data1) * F.template integrate< true >( childIntegrator , _off , off ); + } +#pragma omp atomic + cumulative2[ _node->nodeData.nodeIndex ] += _dot; + } + } + } + } + // Update the dot-product using the cumulative constraints @(depth-1) +#pragma omp parallel for num_threads( threads ) reduction( + : dot ) + for( int i=_sNodesBegin(d-1) ; i<_sNodesEnd(d-1) ; i++ ) + { + const TreeOctNode* node = _sNodes.treeNodes[i]; + const Real* _data2; + if( isValidFEMNode< FEMDegree2 , FEMBType2 >( node ) && ( _data2=coefficients2( node ) ) ) dot += cumulative2[ node->nodeData.nodeIndex ] * (*_data2); + } + + // Down-sample the cumulative constraints from @(depth-1) to @(depth-2) for the next pass + if( d-1>0 ) _downSample< Real , FEMDegree2 , FEMBType2 >( d-1 , cumulative2 ); + } + } + + if( iInfo ) + { + MultiThreadedEvaluator< FEMDegree1 , FEMBType1 > mt1( this , coefficients1 , threads ); + MultiThreadedEvaluator< FEMDegree2 , FEMBType2 > mt2( this , coefficients2 , threads ); + +#pragma omp parallel for num_threads( threads ) reduction( + : dot ) + for( int i=_sNodesBegin(0) ; i<_sNodesEnd(_maxDepth) ; i++ ) + { + if( _isValidSpaceNode( _sNodes.treeNodes[i] ) && !_isValidSpaceNode( _sNodes.treeNodes[i]->children ) && (*iInfo)( _sNodes.treeNodes[i] ) ) + { + + const PointData< Real , HasGradients >& pData = *( (*iInfo)( _sNodes.treeNodes[i] ) ); +#if POINT_DATA_RES + for( int c=0 ; c::SAMPLES ; c++ ) if( pData[c].weight ) + { + Point3D< Real > p = pData[c].position; + Real w = pData[c].weight; + if( HasGradients ) + { + std::pair< Real , Point3D< Real > > v1 = mt1.valueAndGradient( p , omp_get_thread_num() ); + std::pair< Real , Point3D< Real > > v2 = mt2.valueAndGradient( p , omp_get_thread_num() ); + dot += v1.first * v2.first * w * iInfo->valueWeight + Point3D< Real >::Dot( v1.second , v2.second ) * w * iInfo->gradientWeight; + } + else dot += mt1.value( p , omp_get_thread_num() ) * mt2.value( p , omp_get_thread_num() ) * w * iInfo->valueWeight; + } +#else // !POINT_DATA_RES + Point3D< Real > p = pData.position; + Real w = pData.weight; + if( HasGradients ) + { + std::pair< Real , Point3D< Real > > v1 = mt1.valueAndGradient( p , omp_get_thread_num() ); + std::pair< Real , Point3D< Real > > v2 = mt2.valueAndGradient( p , omp_get_thread_num() ); + dot += v1.first * v2.first * w * iInfo->valueWeight + Point3D< Real >::Dot( v1.second , v2.second ) * w * iInfo->gradientWeight; + } + else dot += mt1.value( p , omp_get_thread_num() ) * mt2.value( p , omp_get_thread_num() ) * w * iInfo->valueWeight; +#endif // POINT_DATA_RES + } + } + } + + return dot; +} diff --git a/src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.WeightedSamples.inl b/src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.WeightedSamples.inl new file mode 100644 index 000000000..1ab0fc4aa --- /dev/null +++ b/src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.WeightedSamples.inl @@ -0,0 +1,450 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. Redistributions in binary form must reproduce +the above copyright notice, this list of conditions and the following disclaimer +in the documentation and/or other materials provided with the distribution. + +Neither the name of the Johns Hopkins University nor the names of its contributors +may be used to endorse or promote products derived from this software without specific +prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +// evaluate the result of splatting along a plane and then evaluating at a point on the plane. +template< int Degree > double GetScaleValue( void ) +{ + double centerValues[Degree+1]; + Polynomial< Degree >::BSplineComponentValues( 0.5 , centerValues ); + double scaleValue = 0; + for( int i=0 ; i<=Degree ; i++ ) scaleValue += centerValues[i] * centerValues[i]; + return 1./ scaleValue; +} +template< class Real > +template< int WeightDegree > +void Octree< Real >::_addWeightContribution( SparseNodeData< Real , WeightDegree >& densityWeights , TreeOctNode* node , Point3D< Real > position , PointSupportKey< WeightDegree >& weightKey , Real weight ) +{ + static const double ScaleValue = GetScaleValue< WeightDegree >(); + double dx[ DIMENSION ][ PointSupportKey< WeightDegree >::Size ]; + typename TreeOctNode::Neighbors< PointSupportKey< WeightDegree >::Size >& neighbors = weightKey.template getNeighbors< true >( node , _NodeInitializer ); + densityWeights.reserve( NodeCount() ); + Point3D< Real > start; + Real w; + _startAndWidth( node , start , w ); + for( int dim=0 ; dim::BSplineComponentValues( ( position[dim]-start[dim] ) / w , dx[dim] ); + + weight *= (Real)ScaleValue; + + for( int i=0 ; i::Size ; i++ ) for( int j=0 ; j::Size ; j++ ) + { + double dxdy = dx[0][i] * dx[1][j] * weight; + TreeOctNode** _neighbors = neighbors.neighbors[i][j]; + for( int k=0 ; k::Size ; k++ ) if( _neighbors[k] ) densityWeights[ _neighbors[k] ] += Real( dxdy * dx[2][k] ); + } +} + +template< class Real > +template< int WeightDegree , class PointSupportKey > +Real Octree< Real >::_getSamplesPerNode( const SparseNodeData< Real , WeightDegree >& densityWeights , const TreeOctNode* node , Point3D< Real > position , PointSupportKey& weightKey ) const +{ + Real weight = 0; + double dx[ DIMENSION ][ PointSupportKey::Size ]; + const typename PointSupportKey::template Neighbors< PointSupportKey::Size >& neighbors = weightKey.getNeighbors( node ); + + Point3D< Real > start; + Real w; + _startAndWidth( node , start , w ); + + for( int dim=0 ; dim::BSplineComponentValues( ( position[dim]-start[dim] ) / w , dx[dim] ); + + for( int i=0 ; i +template< int WeightDegree , class PointSupportKey > +void Octree< Real >::_getSampleDepthAndWeight( const SparseNodeData< Real , WeightDegree >& densityWeights , const TreeOctNode* node , Point3D< Real > position , PointSupportKey& weightKey , Real& depth , Real& weight ) const +{ + const TreeOctNode* temp = node; + weight = _getSamplesPerNode( densityWeights , temp , position , weightKey ); + if( weight>=(Real)1. ) depth = Real( _localDepth( temp ) + log( weight ) / log(double(1<<(DIMENSION-1))) ); + else + { + Real oldWeight , newWeight; + oldWeight = newWeight = weight; + while( newWeight<(Real)1. && temp->parent ) + { + temp=temp->parent; + oldWeight = newWeight; + newWeight = _getSamplesPerNode( densityWeights , temp , position , weightKey ); + } + if( newWeight<=0 || oldWeight<=0 ) + { + fprintf( stderr , "[WARNING] Octree::_getSampleDepthAndWeight: Weights should be positive: %g %g\n" , oldWeight , newWeight ); + depth = weight = 0; + return; + } + else depth = Real( _localDepth( temp ) + log( newWeight ) / log( newWeight / oldWeight ) ); + } + weight = Real( pow( double(1<<(DIMENSION-1)) , -double(depth) ) ); +} +template< class Real > +template< int WeightDegree , class PointSupportKey > +void Octree< Real >::_getSampleDepthAndWeight( const SparseNodeData< Real , WeightDegree >& densityWeights , Point3D< Real > position , PointSupportKey& weightKey , Real& depth , Real& weight ) const +{ + TreeOctNode* temp; + Point3D< Real > myCenter( (Real)0.5 , (Real)0.5 , (Real)0.5 ); + Real myWidth = Real( 1. ); + + // Get the finest node with depth less than or equal to the splat depth that contains the point + temp = _spaceRoot; + while( true ) + { + if( !IsActiveNode( temp->children ) ) break;// fprintf( stderr , "[ERROR] Octree::GetSampleDepthAndWeight\n" ) , exit( 0 ); + int cIndex = TreeOctNode::CornerIndex( myCenter , position ); + if( densityWeights( temp->children + cIndex ) ) temp = temp->children + cIndex; + else break; + myWidth /= 2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + } + return _getSampleDepthAndWeight( densityWeights , temp , position , weightKey , depth , weight ); +} + +template< class Real > +template< bool CreateNodes , int DataDegree , class V > +void Octree< Real >::_splatPointData( TreeOctNode* node , Point3D< Real > position , V v , SparseNodeData< V , DataDegree >& dataInfo , PointSupportKey< DataDegree >& dataKey ) +{ + double dx[ DIMENSION ][ PointSupportKey< DataDegree >::Size ]; + typename TreeOctNode::Neighbors< PointSupportKey< DataDegree >::Size >& neighbors = dataKey.template getNeighbors< CreateNodes >( node , _NodeInitializer ); + Point3D< Real > start; + Real w; + _startAndWidth( node , start , w ); + + for( int dd=0 ; dd::BSplineComponentValues( ( position[dd]-start[dd] ) / w , dx[dd] ); + + for( int i=0 ; i::Size ; i++ ) for( int j=0 ; j::Size ; j++ ) + { + double dxdy = dx[0][i] * dx[1][j]; + for( int k=0 ; k::Size ; k++ ) + if( IsActiveNode( neighbors.neighbors[i][j][k] ) ) + { + TreeOctNode* _node = neighbors.neighbors[i][j][k]; + + double dxdydz = dxdy * dx[2][k]; + dataInfo[ _node ] += v * (Real)dxdydz; + } + } +} +template< class Real > +template< bool CreateNodes , int WeightDegree , int DataDegree , class V > +Real Octree< Real >::_splatPointData( const SparseNodeData< Real , WeightDegree >& densityWeights , Point3D< Real > position , V v , SparseNodeData< V , DataDegree >& dataInfo , PointSupportKey< WeightDegree >& weightKey , PointSupportKey< DataDegree >& dataKey , LocalDepth minDepth , LocalDepth maxDepth , int dim ) +{ + double dx; + V _v; + TreeOctNode* temp; + int cnt=0; + double width; + Point3D< Real > myCenter( (Real)0.5 , (Real)0.5 , (Real)0.5 ); + Real myWidth = (Real)1.; + + temp = _spaceRoot; + while( true ) + { + if( !IsActiveNode( temp->children ) ) break; + int cIndex = TreeOctNode::CornerIndex( myCenter , position ); + if( densityWeights( temp->children + cIndex ) ) temp = temp->children + cIndex; + else break; + myWidth /= 2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + } + Real weight , depth; + _getSampleDepthAndWeight( densityWeights , temp , position , weightKey , depth , weight ); + + if( depthmaxDepth ) depth = Real(maxDepth); + int topDepth = int(ceil(depth)); + + dx = 1.0-(topDepth-depth); + if ( topDepth<=minDepth ) topDepth = minDepth , dx = 1; + else if( topDepth> maxDepth ) topDepth = maxDepth , dx = 1; + + while( _localDepth( temp )>topDepth ) temp=temp->parent; + while( _localDepth( temp )children ) temp->initChildren( _NodeInitializer ); + int cIndex = TreeOctNode::CornerIndex( myCenter , position ); + temp = &temp->children[cIndex]; + myWidth/=2; + if( cIndex&1 ) myCenter[0] += myWidth/2; + else myCenter[0] -= myWidth/2; + if( cIndex&2 ) myCenter[1] += myWidth/2; + else myCenter[1] -= myWidth/2; + if( cIndex&4 ) myCenter[2] += myWidth/2; + else myCenter[2] -= myWidth/2; + } + width = 1.0 / ( 1<<_localDepth( temp ) ); + _v = v * weight / Real( pow( width , dim ) ) * Real( dx ); + _splatPointData< CreateNodes >( temp , position , _v , dataInfo , dataKey ); + if( fabs(1.0-dx) > EPSILON ) + { + dx = Real(1.0-dx); + temp = temp->parent; + width = 1.0 / ( 1<<_localDepth( temp ) ); + + _v = v * weight / Real( pow( width , dim ) ) * Real( dx ); + _splatPointData< CreateNodes >( temp , position , _v , dataInfo , dataKey ); + } + return weight; +} +template< class Real > +template< bool CreateNodes , int WeightDegree , int DataDegree , class V > +Real Octree< Real >::_multiSplatPointData( const SparseNodeData< Real , WeightDegree >* densityWeights , TreeOctNode* node , Point3D< Real > position , V v , SparseNodeData< V , DataDegree >& dataInfo , PointSupportKey< WeightDegree >& weightKey , PointSupportKey< DataDegree >& dataKey , int dim ) +{ + Real _depth , weight; + if( densityWeights ) _getSampleDepthAndWeight( *densityWeights , position , weightKey , _depth , weight ); + else weight = (Real)1.; + V _v = v * weight; + + double dx[ DIMENSION ][ PointSupportKey< DataDegree >::Size ]; + dataKey.template getNeighbors< CreateNodes >( node , _NodeInitializer ); + + for( TreeOctNode* _node=node ; _localDepth( _node )>=0 ; _node=_node->parent ) + { + V __v = _v * (Real)pow( 1<<_localDepth( _node ) , dim ); + Point3D< Real > start; + Real w; + _startAndWidth( _node , start , w ); + for( int dd=0 ; dd::BSplineComponentValues( ( position[dd]-start[dd] ) / w , dx[dd] ); + typename TreeOctNode::Neighbors< PointSupportKey< DataDegree >::Size >& neighbors = dataKey.neighbors[ _localToGlobal( _localDepth( _node ) ) ]; + for( int i=0 ; i::Size ; i++ ) for( int j=0 ; j::Size ; j++ ) + { + double dxdy = dx[0][i] * dx[1][j]; + for( int k=0 ; k::Size ; k++ ) + if( IsActiveNode( neighbors.neighbors[i][j][k] ) ) + { + TreeOctNode* _node = neighbors.neighbors[i][j][k]; + double dxdydz = dxdy * dx[2][k]; + dataInfo[ _node ] += __v * (Real)dxdydz; + } + } + } + return weight; +} + +template< class Real > +template< class V , int DataDegree , BoundaryType BType , class Coefficients > +V Octree< Real >::_evaluate( const Coefficients& coefficients , Point3D< Real > p , const BSplineData< DataDegree , BType >& bsData , const ConstPointSupportKey< DataDegree >& dataKey ) const +{ + V value = V(0); + + for( int d=_localToGlobal( 0 ) ; d<=dataKey.depth() ; d++ ) + { + double dx[ DIMENSION ][ PointSupportKey< DataDegree >::Size ]; + memset( dx , 0 , sizeof( double ) * DIMENSION * PointSupportKey< DataDegree >::Size ); + { + const TreeOctNode* n = dataKey.neighbors[d].neighbors[ PointSupportKey< DataDegree >::LeftRadius ][ PointSupportKey< DataDegree >::LeftRadius ][ PointSupportKey< DataDegree >::LeftRadius ]; + if( !n ) fprintf( stderr , "[ERROR] Point is not centered on a node\n" ) , exit( 0 ); + int fIdx[3]; + functionIndex< DataDegree , BType >( n , fIdx ); + int fStart , fEnd; + BSplineData< DataDegree , BType >::FunctionSpan( _localDepth( n ) , fStart , fEnd ); + for( int dd=0 ; dd::LeftRadius ; i<=PointSupportKey< DataDegree >::RightRadius ; i++ ) + if( fIdx[dd]+i>=fStart && fIdx[dd]+i::RightRadius ]( p[dd] ); + } + for( int i=0 ; i::Size ; i++ ) for( int j=0 ; j::Size ; j++ ) for( int k=0 ; k::Size ; k++ ) + { + const TreeOctNode* n = dataKey.neighbors[d].neighbors[i][j][k]; + if( isValidFEMNode< DataDegree , BType >( n ) ) + { + const V* v = coefficients( n ); + if( v ) value += (*v) * (Real) ( dx[0][i] * dx[1][j] * dx[2][k] ); + } + } + } + + return value; +} + +template< class Real > +template< class V , int DataDegree , BoundaryType BType > +Pointer( V ) Octree< Real >::voxelEvaluate( const DenseNodeData< V , DataDegree >& coefficients , int& res , Real isoValue , LocalDepth depth , bool primal ) +{ + int begin , end , dim; + if( depth<=0 || depth>_maxDepth ) depth = _maxDepth; + + // Initialize the coefficients at the coarsest level + Pointer( V ) _coefficients = NullPointer( V ); + { + LocalDepth d = 0; + begin = _BSplineBegin< DataDegree , BType >( d ) , end = _BSplineEnd< DataDegree , BType >( d ) , dim = end - begin; + _coefficients = NewPointer< V >( dim * dim * dim ); + memset( _coefficients , 0 , sizeof( V ) * dim * dim * dim ); +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodesBegin(d) ; i<_sNodesEnd(d) ; i++ ) if( !_outOfBounds< DataDegree , BType >( _sNodes.treeNodes[i] ) ) + { + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( _sNodes.treeNodes[i] , _d , _off ); + _off[0] -= begin , _off[1] -= begin , _off[2] -= begin; + _coefficients[ _off[0] + _off[1]*dim + _off[2]*dim*dim ] = coefficients[i]; + } + } + + // Up-sample and add in the existing coefficients + for( LocalDepth d=1 ; d<=depth ; d++ ) + { + begin = _BSplineBegin< DataDegree , BType >( d ) , end = _BSplineEnd< DataDegree , BType >( d ) , dim = end - begin; + Pointer( V ) __coefficients = NewPointer< V >( dim * dim *dim ); + memset( __coefficients , 0 , sizeof( V ) * dim * dim * dim ); +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodesBegin(d) ; i<_sNodesEnd(d) ; i++ ) if( !_outOfBounds< DataDegree , BType >( _sNodes.treeNodes[i] ) ) + { + LocalDepth _d ; LocalOffset _off; + _localDepthAndOffset( _sNodes.treeNodes[i] , _d , _off ); + _off[0] -= begin , _off[1] -= begin , _off[2] -= begin; + __coefficients[ _off[0] + _off[1]*dim + _off[2]*dim*dim ] = coefficients[i]; + } + _UpSample< V , DataDegree , BType >( d , ( ConstPointer(V) )_coefficients , __coefficients , threads ); + DeletePointer( _coefficients ); + _coefficients = __coefficients; + } + + res = 1<( res*res*res ); + memset( values , 0 , sizeof(V)*res*res*res ); + + if( primal ) + { + // evaluate at the cell corners + typename BSplineEvaluationData< DataDegree , BType >::CornerEvaluator::Evaluator evaluator; + BSplineEvaluationData< DataDegree , BType >::SetCornerEvaluator( evaluator , depth ); +#pragma omp parallel for num_threads( threads ) + for( int k=0 ; k::CornerEnd ; kk<=-BSplineSupportSizes< DataDegree >::CornerStart ; kk++ ) if( k+kk>=begin && k+kk::CornerEnd ; jj<=-BSplineSupportSizes< DataDegree >::CornerStart ; jj++ ) if( j+jj>=begin && j+jj::CornerEnd ; ii<=-BSplineSupportSizes< DataDegree >::CornerStart ; ii++ ) if( i+ii>=begin && i+ii::CenterEvaluator::Evaluator evaluator; + BSplineEvaluationData< DataDegree , BType >::SetCenterEvaluator( evaluator , depth ); +#pragma omp parallel for num_threads( threads ) + for( int k=0 ; k::SupportEnd ; kk<=-BSplineSupportSizes< DataDegree >::SupportStart ; kk++ ) if( k+kk>=begin && k+kk::SupportEnd ; jj<=-BSplineSupportSizes< DataDegree >::SupportStart ; jj++ ) if( j+jj>=begin && j+jj::SupportEnd ; ii<=-BSplineSupportSizes< DataDegree >::SupportStart ; ii++ ) if( i+ii>=begin && i+ii +template< int FEMDegree , BoundaryType BType > +SparseNodeData< Real , 0 > Octree< Real >::leafValues( const DenseNodeData< Real , FEMDegree >& coefficients ) const +{ + SparseNodeData< Real , 0 > values; + DenseNodeData< Real , FEMDegree > _coefficients( _sNodesEnd(_maxDepth-1) ); + memset( &_coefficients[0] , 0 , sizeof(Real)*_sNodesEnd(_maxDepth-1) ); + for( int i=_sNodes.begin( _localToGlobal( 0 ) ) ; i<_sNodesEnd(_maxDepth-1) ; i++ ) _coefficients[i] = coefficients[i]; + for( LocalDepth d=1 ; d<_maxDepth ; d++ ) _upSample( d , _coefficients ); + for( LocalDepth d=_maxDepth ; d>=0 ; d-- ) + { + _Evaluator< FEMDegree , BType > evaluator; + evaluator.set( d ); + std::vector< ConstPointSupportKey< FEMDegree > > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i& neighborKey = neighborKeys[ omp_get_thread_num() ]; + TreeOctNode* node = _sNodes.treeNodes[i]; + if( !IsActiveNode( node->children ) ) + { + neighborKey.getNeighbors( node ); + bool isInterior = _IsInteriorlySupported< FEMDegree >( node->parent ); + values[ node ] = _getCenterValue( neighborKey , node , coefficients , _coefficients , evaluator , isInterior ); + } + } + } + return values; +} +template< class Real > +template< int FEMDegree , BoundaryType BType > +SparseNodeData< Point3D< Real > , 0 > Octree< Real >::leafGradients( const DenseNodeData< Real , FEMDegree >& coefficients ) const +{ + SparseNodeData< Point3D< Real > , 0 > gradients; + DenseNodeData< Real , FEMDegree > _coefficients( _sNodesEnd(_maxDepth-1 ) ); + memset( &_coefficients[0] , 0 , sizeof(Real)*_sNodesEnd(_maxDepth-1) ); + for( int i=_sNodesBegin(0) ; i<_sNodesEnd(_maxDepth-1) ; i++ ) _coefficients[i] = coefficients[i]; + for( LocalDepth d=1 ; d<_maxDepth ; d++ ) _upSample( d , _coefficients ); + for( LocalDepth d=_maxDepth ; d>=0 ; d-- ) + { + _Evaluator< FEMDegree , BType > evaluator; + evaluator.set( d ); + std::vector< ConstPointSupportKey< FEMDegree > > neighborKeys( std::max< int >( 1 , threads ) ); + for( size_t i=0 ; i& neighborKey = neighborKeys[ omp_get_thread_num() ]; + TreeOctNode* node = _sNodes.treeNodes[i]; + if( !IsActiveNode( node->children ) ) + { + neighborKey.getNeighbors( node ); + bool isInterior = _IsInteriorlySupported< FEMDegree >( node->parent ); + gradients[ node ] = _getCenterValueAndGradient( neighborKey , node , coefficients , _coefficients , evaluator , isInterior ).second; + } + } + } + return gradients; +} diff --git a/src/meshlabplugins/filter_screened_poisson/Src/MyTime.h b/src/meshlabplugins/filter_screened_poisson/Src/MyTime.h new file mode 100644 index 000000000..09c084e00 --- /dev/null +++ b/src/meshlabplugins/filter_screened_poisson/Src/MyTime.h @@ -0,0 +1,51 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. Redistributions in binary form must reproduce +the above copyright notice, this list of conditions and the following disclaimer +in the documentation and/or other materials provided with the distribution. + +Neither the name of the Johns Hopkins University nor the names of its contributors +may be used to endorse or promote products derived from this software without specific +prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#ifndef MY_TIME_INCLUDED +#define MY_TIME_INCLUDED + +#include +#include +#ifndef WIN32 +#include +#endif // WIN32 + +inline double Time( void ) +{ +#ifdef WIN32 + struct _timeb t; + _ftime( &t ); + return double( t.time ) + double( t.millitm ) / 1000.0; +#else // WIN32 + struct timeval t; + gettimeofday( &t , NULL ); + return t.tv_sec + double( t.tv_usec ) / 1000000; +#endif // WIN32 +} + +#endif // MY_TIME_INCLUDED diff --git a/src/meshlabplugins/filter_screened_poisson/Src/PPolynomial.h b/src/meshlabplugins/filter_screened_poisson/Src/PPolynomial.h index fc0651460..1b6b1c206 100755 --- a/src/meshlabplugins/filter_screened_poisson/Src/PPolynomial.h +++ b/src/meshlabplugins/filter_screened_poisson/Src/PPolynomial.h @@ -30,57 +30,56 @@ DAMAGE. #define P_POLYNOMIAL_INCLUDED #include #include "Polynomial.h" +#include "Array.h" -template -class StartingPolynomial{ +template< int Degree > +class StartingPolynomial +{ public: - Polynomial p; + Polynomial< Degree > p; double start; - template - StartingPolynomial operator * (const StartingPolynomial& p) const; - StartingPolynomial scale(double s) const; - StartingPolynomial shift(double t) const; - int operator < (const StartingPolynomial& sp) const; - static int Compare(const void* v1,const void* v2); + template< int Degree2 > + StartingPolynomial< Degree+Degree2 > operator * ( const StartingPolynomial< Degree2 >& p ) const; + StartingPolynomial scale( double s ) const; + StartingPolynomial shift( double t ) const; + int operator < ( const StartingPolynomial& sp ) const; + static int Compare( const void* v1 , const void* v2 ); }; -template +template< int Degree > class PPolynomial { public: size_t polyCount; - StartingPolynomial* polys; + Pointer( StartingPolynomial< Degree > ) polys; - PPolynomial(void); - PPolynomial(const PPolynomial& p); - ~PPolynomial(void); + PPolynomial( void ); + PPolynomial( const PPolynomial& p ); + ~PPolynomial( void ); - PPolynomial& operator = (const PPolynomial& p); + PPolynomial& operator = ( const PPolynomial& p ); - int size(void) const; + int size( void ) const; void set( size_t size ); // Note: this method will sort the elements in sps - void set( StartingPolynomial* sps , int count ); + void set( Pointer( StartingPolynomial ) sps , int count ); void reset( size_t newSize ); + PPolynomial& compress( double delta=0. ); double operator()( double t ) const; double integral( double tMin , double tMax ) const; double Integral( void ) const; - template - PPolynomial& operator = (const PPolynomial& p); + template< int Degree2 > PPolynomial< Degree >& operator = ( const PPolynomial< Degree2 >& p ); - PPolynomial operator + (const PPolynomial& p) const; - PPolynomial operator - (const PPolynomial& p) const; + PPolynomial operator + ( const PPolynomial& p ) const; + PPolynomial operator - ( const PPolynomial& p ) const; - template - PPolynomial operator * (const Polynomial& p) const; - - template - PPolynomial operator * (const PPolynomial& p) const; + template< int Degree2 > PPolynomial< Degree+Degree2 > operator * ( const Polynomial< Degree2 >& p ) const; + template< int Degree2 > PPolynomial< Degree+Degree2 > operator * ( const PPolynomial< Degree2 >& p) const; PPolynomial& operator += ( double s ); @@ -92,15 +91,16 @@ public: PPolynomial operator * ( double s ) const; PPolynomial operator / ( double s ) const; - PPolynomial& addScaled(const PPolynomial& poly,double scale); + PPolynomial& addScaled( const PPolynomial& poly , double scale ); PPolynomial scale( double s ) const; PPolynomial shift( double t ) const; + PPolynomial reflect( double r=0. ) const; PPolynomial< Degree-1 > derivative(void) const; PPolynomial< Degree+1 > integral(void) const; - void getSolutions(double c,std::vector& roots,double EPS,double min=-DBL_MAX,double max=DBL_MAX) const; + void getSolutions( double c , std::vector< double >& roots , double EPS , double min=-DBL_MAX , double max=DBL_MAX ) const; void printnl( void ) const; diff --git a/src/meshlabplugins/filter_screened_poisson/Src/Ply.h b/src/meshlabplugins/filter_screened_poisson/Src/Ply.h index 934cb3615..699381ffb 100755 --- a/src/meshlabplugins/filter_screened_poisson/Src/Ply.h +++ b/src/meshlabplugins/filter_screened_poisson/Src/Ply.h @@ -33,6 +33,8 @@ #ifndef __PLY_H__ #define __PLY_H__ +#define USE_PLY_WRAPPER 1 + #ifndef WIN32 #define _strdup strdup #endif @@ -239,12 +241,22 @@ static PlyProperty face_props[] = { _strdup( "vertex_indices" ) , PLY_INT , PLY_INT , offsetof( PlyFace , vertices ) , 1 , PLY_UCHAR, PLY_UCHAR , offsetof(PlyFace,nr_vertices) }, }; + +/////////////////// +// PlyVertexType // +/////////////////// + +// The "Wrapper" class indicates the class to cast to/from in order to support linear operations. template< class Real > class PlyVertex { public: - const static int Components=3; - static PlyProperty Properties[]; + typedef PlyVertex Wrapper; + + const static int ReadComponents=3; + const static int WriteComponents=3; + static PlyProperty ReadProperties[]; + static PlyProperty WriteProperties[]; Point3D< Real > point; @@ -260,7 +272,13 @@ public: template< class _Real > PlyVertex& operator /= ( _Real s ) { point /= s ; return *this; } }; template< class Real , class _Real > PlyVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyVertex< Real > v ) { return PlyVertex< Real >( xForm * v.point ); } -template< class Real > PlyProperty PlyVertex< Real >::Properties[]= +template< class Real > PlyProperty PlyVertex< Real >::ReadProperties[]= +{ + { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 } +}; +template< class Real > PlyProperty PlyVertex< Real >::WriteProperties[]= { { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 }, { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 }, @@ -270,8 +288,12 @@ template< class Real > class PlyValueVertex { public: - const static int Components=4; - static PlyProperty Properties[]; + typedef PlyValueVertex Wrapper; + + const static int ReadComponents=4; + const static int WriteComponents=4; + static PlyProperty ReadProperties[]; + static PlyProperty WriteProperties[]; Point3D point; Real value; @@ -288,8 +310,14 @@ public: template< class _Real > PlyValueVertex& operator /= ( _Real s ) { point /= s , value /= Real(s) ; return *this; } }; template< class Real , class _Real > PlyValueVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyValueVertex< Real > v ) { return PlyValueVertex< Real >( xForm * v.point , v.value ); } -template< class Real > -PlyProperty PlyValueVertex< Real >::Properties[]= +template< class Real > PlyProperty PlyValueVertex< Real >::ReadProperties[]= +{ + { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "value" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , value ) ) , 0 , 0 , 0 , 0 } +}; +template< class Real > PlyProperty PlyValueVertex< Real >::WriteProperties[]= { { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 }, { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 }, @@ -300,8 +328,12 @@ template< class Real > class PlyOrientedVertex { public: - const static int Components=6; - static PlyProperty Properties[]; + typedef PlyOrientedVertex Wrapper; + + const static int ReadComponents=6; + const static int WriteComponents=6; + static PlyProperty ReadProperties[]; + static PlyProperty WriteProperties[]; Point3D point , normal; @@ -317,8 +349,16 @@ public: template< class _Real > PlyOrientedVertex& operator /= ( _Real s ) { point /= s , normal /= s ; return *this; } }; template< class Real , class _Real > PlyOrientedVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyOrientedVertex< Real > v ) { return PlyOrientedVertex< Real >( xForm * v.point , xForm.inverse().transpose() * v.normal ); } -template< class Real > -PlyProperty PlyOrientedVertex< Real >::Properties[]= +template< class Real > PlyProperty PlyOrientedVertex< Real >::ReadProperties[]= +{ + { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "nx" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , normal.coords[0] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "ny" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , normal.coords[1] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "nz" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , normal.coords[2] ) ) , 0 , 0 , 0 , 0 } +}; +template< class Real > PlyProperty PlyOrientedVertex< Real >::WriteProperties[]= { { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 }, { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 }, @@ -331,19 +371,61 @@ template< class Real > class PlyColorVertex { public: - const static int Components=6; - static PlyProperty Properties[]; + struct _PlyColorVertex + { + Point3D< Real > point , color; + _PlyColorVertex( void ) { ; } + _PlyColorVertex( Point3D< Real > p , Point3D< Real > c ) : point(p) , color(c) { ; } + _PlyColorVertex( PlyColorVertex< Real > p ){ point = p.point ; for( int c=0 ; c<3 ; c++ ) color[c] = (Real) p.color[c]; } + operator PlyColorVertex< Real > () + { + PlyColorVertex< Real > p; + p.point = point; + for( int c=0 ; c<3 ; c++ ) p.color[c] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( color[c]+0.5 ) ) ); + return p; + } - Point3D point; + _PlyColorVertex operator + ( _PlyColorVertex p ) const { return _PlyColorVertex( point+p.point , color+p.color ); } + _PlyColorVertex operator - ( _PlyColorVertex p ) const { return _PlyColorVertex( point-p.value , color-p.color ); } + template< class _Real > _PlyColorVertex operator * ( _Real s ) const { return _PlyColorVertex( point*s , color*s ); } + template< class _Real > _PlyColorVertex operator / ( _Real s ) const { return _PlyColorVertex( point/s , color/s ); } + _PlyColorVertex& operator += ( _PlyColorVertex p ) { point += p.point , color += p.color ; return *this; } + _PlyColorVertex& operator -= ( _PlyColorVertex p ) { point -= p.point , color -= p.color ; return *this; } + template< class _Real > _PlyColorVertex& operator *= ( _Real s ) { point *= s , color *= s ; return *this; } + template< class _Real > _PlyColorVertex& operator /= ( _Real s ) { point /= s , color /= s ; return *this; } + }; + + typedef _PlyColorVertex Wrapper; + + const static int ReadComponents=9; + const static int WriteComponents=6; + static PlyProperty ReadProperties[]; + static PlyProperty WriteProperties[]; + + Point3D< Real > point; unsigned char color[3]; - operator Point3D& () {return point;} - operator const Point3D& () const {return point;} - PlyColorVertex(void) {point.coords[0]=point.coords[1]=point.coords[2]=0,color[0]=color[1]=color[2]=0;} - PlyColorVertex(const Point3D& p) {point=p;} + operator Point3D< Real >& (){ return point; } + operator const Point3D< Real >& () const { return point; } + PlyColorVertex( void ) { point.coords[0] = point.coords[1] = point.coords[2] = 0 , color[0] = color[1] = color[2] = 0; } + PlyColorVertex( const Point3D& p ) { point=p; } + PlyColorVertex( const Point3D< Real >& p , const unsigned char c[3] ) { point = p , color[0] = c[0] , color[1] = c[1] , color[2] = c[2]; } }; -template< class Real > -PlyProperty PlyColorVertex< Real >::Properties[]= +template< class Real , class _Real > PlyColorVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyColorVertex< Real > v ) { return PlyColorVertex< Real >( xForm * v.point , v.color ); } + +template< class Real > PlyProperty PlyColorVertex< Real >::ReadProperties[]= +{ + { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "z" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "red" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[0] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[1] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "blue" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[2] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "r" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[0] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "g" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[1] ) ) , 0 , 0 , 0 , 0 }, + { _strdup( "b" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[2] ) ) , 0 , 0 , 0 , 0 } +}; +template< class Real > PlyProperty PlyColorVertex< Real >::WriteProperties[]= { { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 }, { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 }, @@ -352,6 +434,77 @@ PlyProperty PlyColorVertex< Real >::Properties[]= { _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[1] ) ) , 0 , 0 , 0 , 0 }, { _strdup( "blue" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[2] ) ) , 0 , 0 , 0 , 0 } }; +template< class Real > +class PlyColorAndValueVertex +{ +public: + struct _PlyColorAndValueVertex + { + Point3D< Real > point , color; + Real value; + _PlyColorAndValueVertex( void ) : value(0) { ; } + _PlyColorAndValueVertex( Point3D< Real > p , Point3D< Real > c , Real v ) : point(p) , color(c) , value(v) { ; } + _PlyColorAndValueVertex( PlyColorAndValueVertex< Real > p ){ point = p.point ; for( int c=0 ; c<3 ; c++ ) color[c] = (Real) p.color[c] ; value = p.value; } + operator PlyColorAndValueVertex< Real > () + { + PlyColorAndValueVertex< Real > p; + p.point = point; + for( int c=0 ; c<3 ; c++ ) p.color[c] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( color[c]+0.5 ) ) ); + p.value = value; + return p; + } + + _PlyColorAndValueVertex operator + ( _PlyColorAndValueVertex p ) const { return _PlyColorAndValueVertex( point+p.point , color+p.color , value+p.value ); } + _PlyColorAndValueVertex operator - ( _PlyColorAndValueVertex p ) const { return _PlyColorAndValueVertex( point-p.value , color-p.color , value+p.value ); } + template< class _Real > _PlyColorAndValueVertex operator * ( _Real s ) const { return _PlyColorAndValueVertex( point*s , color*s , value*s ); } + template< class _Real > _PlyColorAndValueVertex operator / ( _Real s ) const { return _PlyColorAndValueVertex( point/s , color/s , value/s ); } + _PlyColorAndValueVertex& operator += ( _PlyColorAndValueVertex p ) { point += p.point , color += p.color , value += p.value ; return *this; } + _PlyColorAndValueVertex& operator -= ( _PlyColorAndValueVertex p ) { point -= p.point , color -= p.color , value -= p.value ; return *this; } + template< class _Real > _PlyColorAndValueVertex& operator *= ( _Real s ) { point *= s , color *= s , value *= (Real)s ; return *this; } + template< class _Real > _PlyColorAndValueVertex& operator /= ( _Real s ) { point /= s , color /= s , value /= (Real)s ; return *this; } + }; + + typedef _PlyColorAndValueVertex Wrapper; + + const static int ReadComponents=10; + const static int WriteComponents=7; + static PlyProperty ReadProperties[]; + static PlyProperty WriteProperties[]; + + Point3D< Real > point; + unsigned char color[3]; + Real value; + + operator Point3D< Real >& (){ return point; } + operator const Point3D< Real >& () const { return point; } + PlyColorAndValueVertex( void ) { point.coords[0] = point.coords[1] = point.coords[2] = (Real)0 , color[0] = color[1] = color[2] = 0 , value = (Real)0; } + PlyColorAndValueVertex( const Point3D< Real >& p ) { point=p; } + PlyColorAndValueVertex( const Point3D< Real >& p , const unsigned char c[3] , Real v) { point = p , color[0] = c[0] , color[1] = c[1] , color[2] = c[2] , value = v; } +}; +template< class Real , class _Real > PlyColorAndValueVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyColorAndValueVertex< Real > v ) { return PlyColorAndValueVertex< Real >( xForm * v.point , v.color , v.value ); } +template< class Real > PlyProperty PlyColorAndValueVertex< Real >::ReadProperties[]= +{ + { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "value" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , value ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "red" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[0] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[1] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "blue" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[2] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "r" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[0] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "g" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[1] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "b" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[2] ) ) , 0 , 0 , 0 , 0 } +}; +template< class Real > PlyProperty PlyColorAndValueVertex< Real >::WriteProperties[]= +{ + { _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "value" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , value ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "red" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[0] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[1] ) ) , 0 , 0 , 0 , 0 } , + { _strdup( "blue" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[2] ) ) , 0 , 0 , 0 , 0 } +}; template< class Vertex , class Real > int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_type , const Point3D< float >& translate , float scale , char** comments=NULL , int commentNum=0 , XForm4x4< Real > xForm=XForm4x4< Real >::Identity() ); @@ -359,6 +512,94 @@ int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file template< class Vertex , class Real > int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_type , char** comments=NULL , int commentNum=0 , XForm4x4< Real > xForm=XForm4x4< Real >::Identity() ); +inline bool PlyReadHeader( char* fileName , PlyProperty* properties , int propertyNum , bool* readFlags , int& file_type ) +{ + int nr_elems; + char **elist; + float version; + PlyFile* ply; + char* elem_name; + int num_elems; + int nr_props; + PlyProperty** plist; + + ply = ply_open_for_reading( fileName , &nr_elems , &elist , &file_type , &version ); + if( !ply ) return false; + + for( int i=0 ; ielems[i]->name ); + free( ply->elems[i]->store_prop ); + for( int j=0 ; jelems[i]->nprops ; j++ ) + { + free( ply->elems[i]->props[j]->name ); + free( ply->elems[i]->props[j] ); + } + free( ply->elems[i]->props ); + } + for( int i=0 ; ielems[i] ); + free( ply->elems ); + for( int i=0 ; inum_comments ; i++ ) free( ply->comments[i] ); + free( ply->comments ); + for( int i=0 ; inum_obj_info ; i++ ) free( ply->obj_info[i] ); + free( ply->obj_info ); + ply_free_other_elements( ply->other_elems ); + + for( int i=0 ; iname ); + free( plist[j] ); + } + free( plist ); + } // for each type of element + + for( int i=0 ; ielems[i]->name ); + free( ply->elems[i]->store_prop ); + for( int j=0 ; jelems[i]->nprops ; j++ ) + { + free( ply->elems[i]->props[j]->name ); + free( ply->elems[i]->props[j] ); + } + if( ply->elems[i]->props && ply->elems[i]->nprops ) free(ply->elems[i]->props); + } + for( int i=0 ; ielems[i]); + free( ply->elems) ; + for( int i=0 ; inum_comments ; i++ ) free( ply->comments[i] ); + free( ply->comments ); + for( int i=0 ; inum_obj_info ; i++ ) free( ply->obj_info[i] ); + free( ply->obj_info ); + ply_free_other_elements(ply->other_elems); + + + for( int i=0 ; i int PlyReadPolygons(char* fileName, std::vector& vertices,std::vector >& polygons, @@ -626,7 +867,7 @@ int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_ // describe vertex and face properties // ply_element_count( ply , "vertex" , nr_vertices ); - for( int i=0 ; i template< int Degree > -class Polynomial{ +class Polynomial +{ public: double coefficients[Degree+1]; @@ -86,7 +89,11 @@ public: void getSolutions(double c,std::vector& roots,double EPS) const; int getSolutions( double c , double* roots , double EPS ) const; + // [NOTE] Both of these methods define the indexing according to DeBoor's algorithm, so that + // Polynomial< Degree >BSplineComponent( 0 )( 1.0 )=0 for all Degree>0. static Polynomial BSplineComponent( int i ); + static void BSplineComponentValues( double x , double* values ); + static void BinomialCoefficients( int bCoefficients[Degree+1] ); }; #include "Polynomial.inl" diff --git a/src/meshlabplugins/filter_screened_poisson/Src/Polynomial.inl b/src/meshlabplugins/filter_screened_poisson/Src/Polynomial.inl index 6d7ea1d88..ea7ae5044 100755 --- a/src/meshlabplugins/filter_screened_poisson/Src/Polynomial.inl +++ b/src/meshlabplugins/filter_screened_poisson/Src/Polynomial.inl @@ -306,6 +306,7 @@ int Polynomial::getSolutions( double c , double* roots , double EPS ) co for( int i=0 ; i<_rCount ; i++ ) if( fabs(_roots[i][1])<=EPS ) roots[rCount++] = _roots[i][0]; return rCount; } +// The 0-th order B-spline template< > Polynomial< 0 > Polynomial< 0 >::BSplineComponent( int i ) { @@ -313,20 +314,56 @@ Polynomial< 0 > Polynomial< 0 >::BSplineComponent( int i ) p.coefficients[0] = 1.; return p; } + +// The Degree-th order B-spline template< int Degree > Polynomial< Degree > Polynomial< Degree >::BSplineComponent( int i ) { + // B_d^i(x) = \int_x^1 B_{d-1}^{i}(y) dy + \int_0^x B_{d-1}^{i-1} y dy + // = \int_0^1 B_{d-1}^{i}(y) dy - \int_0^x B_{d-1}^{i}(y) dy + \int_0^x B_{d-1}^{i-1} y dy Polynomial p; - if( i>0 ) - { - Polynomial< Degree > _p = Polynomial< Degree-1 >::BSplineComponent( i-1 ).integral(); - p -= _p; - p.coefficients[0] += _p(1); - } if( i _p = Polynomial< Degree-1 >::BSplineComponent( i ).integral(); + p -= _p; + p.coefficients[0] += _p(1); + } + if( i>0 ) + { + Polynomial< Degree > _p = Polynomial< Degree-1 >::BSplineComponent( i-1 ).integral(); p += _p; } return p; -} \ No newline at end of file +} + + +// The 0-th order B-spline values +template< > void Polynomial< 0 >::BSplineComponentValues( double x , double* values ){ values[0] = 1.; } +// The Degree-th order B-spline +template< int Degree > void Polynomial< Degree >::BSplineComponentValues( double x , double* values ) +{ + const double Scale = 1./Degree; + Polynomial< Degree-1 >::BSplineComponentValues( x , values+1 ); + values[0] = values[1] * (1.-x) * Scale; + for( int i=1 ; i void Polynomial< 0 >::BinomialCoefficients( int bCoefficients[1] ){ bCoefficients[0] = 1; } +template< int Degree > void Polynomial< Degree >::BinomialCoefficients( int bCoefficients[Degree+1] ) +{ + Polynomial< Degree-1 >::BinomialCoefficients( bCoefficients ); + int leftValue = 0; + for( int i=0 ; i +#include +#include +#include +#if defined( _WIN32 ) || defined( _WIN64 ) +#include +#include +#endif // _WIN32 || _WIN64 +#include "MyTime.h" +#include "MarchingCubes.h" +#include "Octree.h" +#include "SparseMatrix.h" +#include "CmdLineParser.h" +#include "PPolynomial.h" +#include "Ply.h" +#include "MemoryUsage.h" +#ifdef _OPENMP +#include "omp.h" +#endif // _OPENMP +void DumpOutput( const char* format , ... ); +void DumpOutput2( std::vector< char* >& comments , const char* format , ... ); +#include "MultiGridOctreeData.h" + +#define DEFAULT_FULL_DEPTH 5 + +#define XSTR(x) STR(x) +#define STR(x) #x +#if DEFAULT_FULL_DEPTH +#pragma message ( "[WARNING] Setting default full depth to " XSTR(DEFAULT_FULL_DEPTH) ) +#endif // DEFAULT_FULL_DEPTH + +#include +char* outputFile=NULL; +int echoStdout=0; +void DumpOutput( const char* format , ... ) +{ + if( outputFile ) + { + FILE* fp = fopen( outputFile , "a" ); + va_list args; + va_start( args , format ); + vfprintf( fp , format , args ); + fclose( fp ); + va_end( args ); + } + if( echoStdout ) + { + va_list args; + va_start( args , format ); + vprintf( format , args ); + va_end( args ); + } +} +void DumpOutput2( std::vector< char* >& comments , const char* format , ... ) +{ + if( outputFile ) + { + FILE* fp = fopen( outputFile , "a" ); + va_list args; + va_start( args , format ); + vfprintf( fp , format , args ); + fclose( fp ); + va_end( args ); + } + if( echoStdout ) + { + va_list args; + va_start( args , format ); + vprintf( format , args ); + va_end( args ); + } + comments.push_back( new char[1024] ); + char* str = comments.back(); + va_list args; + va_start( args , format ); + vsprintf( str , format , args ); + va_end( args ); + if( str[strlen(str)-1]=='\n' ) str[strlen(str)-1] = 0; +} + + +cmdLineString + In( "in" ) , + Out( "out" ) , + VoxelGrid( "voxel" ) , + XForm( "xForm" ); + +cmdLineReadable +#if defined( _WIN32 ) || defined( _WIN64 ) + Performance( "performance" ) , +#endif // _WIN32 || _WIN64 + ShowResidual( "showResidual" ) , + NoComments( "noComments" ) , + PolygonMesh( "polygonMesh" ) , + Confidence( "confidence" ) , + NormalWeights( "nWeights" ) , + NonManifold( "nonManifold" ) , + ASCII( "ascii" ) , + Density( "density" ) , + NonLinearFit( "nonLinearFit" ) , + PrimalVoxel( "primalVoxel" ) , +#ifndef FAST_COMPILE + FreeBoundary( "freeBoundary" ) , + Double( "double" ) , +#endif // !FAST_COMPILE + Verbose( "verbose" ); + +cmdLineInt +#ifndef FAST_COMPILE + Degree( "degree" , 2 ) , +#endif // !FAST_COMPILE + Depth( "depth" , 8 ) , + CGDepth( "cgDepth" , 0 ) , + KernelDepth( "kernelDepth" ) , + AdaptiveExponent( "adaptiveExp" , 1 ) , + Iters( "iters" , 8 ) , + VoxelDepth( "voxelDepth" , -1 ) , + FullDepth( "fullDepth" , DEFAULT_FULL_DEPTH ) , + MaxSolveDepth( "maxSolveDepth" ) , + Threads( "threads" , omp_get_num_procs() ); + +cmdLineFloat + Color( "color" , 16.f ) , + SamplesPerNode( "samplesPerNode" , 1.5f ) , + Scale( "scale" , 1.1f ) , + CGSolverAccuracy( "cgAccuracy" , 1e-3f ) , + LowResIterMultiplier( "iterMultiplier" , 1.5f ) , + ValueWeight ( "valueWeight" , 4e-0f ) , + GradientWeight( "gradientWeight" , 1e-3f ) , + BiLapWeight ( "biLapWeight" , 1e-5f ); + + +cmdLineReadable* params[] = +{ +#ifndef FAST_COMPILE + &Degree , &Double , &FreeBoundary , +#endif // !FAST_COMPILE + &In , &Depth , &Out , &XForm , + &Scale , &Verbose , &CGSolverAccuracy , &NoComments , &LowResIterMultiplier , + &KernelDepth , &SamplesPerNode , &Confidence , &NormalWeights , &NonManifold , &PolygonMesh , &ASCII , &ShowResidual , &VoxelDepth , + &BiLapWeight , + &ValueWeight , &GradientWeight , &VoxelGrid , &Threads , &MaxSolveDepth , + &AdaptiveExponent , + &Density , + &FullDepth , + &CGDepth , &Iters , + &Color , + &NonLinearFit , + &PrimalVoxel , +#if defined( _WIN32 ) || defined( _WIN64 ) + &Performance , +#endif // _WIN32 || _WIN64 +}; + + +void ShowUsage( char* ex ) +{ + printf( "Usage: %s\n" , ex ); + printf( "\t --%s \n" , In.name ); + + printf( "\t[--%s ]\n" , Out.name ); + + printf( "\t[--%s ]\n" , VoxelGrid.name ); + +#ifndef FAST_COMPILE + printf( "\t[--%s =%d]\n" , Degree.name , Degree.value ); + +#ifndef FOR_RELEASE + printf( "\t[--%s]\n" , FreeBoundary.name ); +#endif // !FOR_RELEASE +#endif // !FAST_COMPILE + + printf( "\t[--%s =%d]\n" , Depth.name , Depth.value ); + + printf( "\t[--%s =%f]\n" , Scale.name , Scale.value ); + + printf( "\t[--%s =%f]\n" , SamplesPerNode.name, SamplesPerNode.value ); + + printf( "\t[--%s =%.3e]\n" , ValueWeight.name , ValueWeight.value ); + + printf( "\t[--%s =%.3e]\n" , GradientWeight.name , GradientWeight.value ); + + printf( "\t[--%s =%.3e]\n" , BiLapWeight.name , BiLapWeight.value ); + + printf( "\t[--%s]\n" , Confidence.name ); + + printf( "\t[--%s]\n" , NormalWeights.name ); + +#ifndef FOR_RELEASE + printf( "\t[--%s =%d]\n", AdaptiveExponent.name , AdaptiveExponent.value ); +#endif // !FOR_RELEASE + + printf( "\t[--%s =%d]\n" , Iters.name , Iters.value ); + +#ifndef FOR_RELEASE + printf( "\t[--%s =%f]\n" , LowResIterMultiplier.name , LowResIterMultiplier.value ); +#endif // FOR_RELEASE + + printf( "\t[--%s =%d]\n" , CGDepth.name , CGDepth.value ); + +#ifndef FOR_RELEASE + printf( "\t[--%s =%g]\n" , CGSolverAccuracy.name , CGSolverAccuracy.value ); +#endif // !FOR_RELEASE + + printf( "\t[--%s =%d]\n" , FullDepth.name , FullDepth.value ); + + printf( "\t[--%s =<%s>]\n" , VoxelDepth.name , Depth.name ); + + printf( "\t[--%s]\n" , PrimalVoxel.name ); + + printf( "\t[--%s ]\n" , Color.name ); + + printf( "\t[--%s]\n" , Density.name ); + + printf( "\t[--%s]\n" , NonLinearFit.name ); + + printf( "\t[--%s]\n" , PolygonMesh.name); + +#ifndef FOR_RELEASE + printf( "\t[--%s]\n" , NonManifold.name ); +#endif // !FOR_RELEASE + +#ifdef _OPENMP + printf( "\t[--%s =%d]\n" , Threads.name , Threads.value ); +#endif // _OPENMP + + printf( "\t[--%s]\n" , Verbose.name ); + +#ifndef FOR_RELEASE +#if defined( _WIN32 ) || defined( _WIN64 ) + printf( "\t[--%s]\n" , Performance.name ); +#endif // _WIN32 || _WIN64 + +#endif // !FOR_RELEASE +#ifndef FOR_RELEASE + printf( "\t[--%s]\n" , ASCII.name ); + + printf( "\t[--%s]\n" , NoComments.name ); +#endif // !FOR_RELEASE + +#ifndef FAST_COMPILE + printf( "\t[--%s]\n" , Double.name ); +#endif // !FAST_COMPILE +} + +template< class Real > +struct ColorInfo +{ + static Point3D< Real > ReadASCII( FILE* fp ) + { + Point3D< unsigned char > c; + if( fscanf( fp , " %c %c %c " , &c[0] , &c[1] , &c[2] )!=3 ) fprintf( stderr , "[ERROR] Failed to read color\n" ) , exit( 0 ); + return Point3D< Real >( (Real)c[0] , (Real)c[1] , (Real)c[2] ); + }; + static bool ValidPlyProperties( const bool* props ){ return ( props[0] || props[3] ) && ( props[1] || props[4] ) && ( props[2] || props[5] ); } + const static PlyProperty PlyProperties[]; +}; +template<> +const PlyProperty ColorInfo< float >::PlyProperties[] = +{ + { "r" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[0] ) ) , 0 , 0 , 0 , 0 } , + { "g" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[1] ) ) , 0 , 0 , 0 , 0 } , + { "b" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[2] ) ) , 0 , 0 , 0 , 0 } , + { "red" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[0] ) ) , 0 , 0 , 0 , 0 } , + { "green" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[1] ) ) , 0 , 0 , 0 , 0 } , + { "blue" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[2] ) ) , 0 , 0 , 0 , 0 } +}; +template<> +const PlyProperty ColorInfo< double >::PlyProperties[] = +{ + { "r" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[0] ) ) , 0 , 0 , 0 , 0 } , + { "g" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[1] ) ) , 0 , 0 , 0 , 0 } , + { "b" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[2] ) ) , 0 , 0 , 0 , 0 } , + { "red" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[0] ) ) , 0 , 0 , 0 , 0 } , + { "green" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[1] ) ) , 0 , 0 , 0 , 0 } , + { "blue" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[2] ) ) , 0 , 0 , 0 , 0 } +}; + +bool ValidPlyColorProperties( const bool* props ){ return ( props[0] || props[3] ) && ( props[1] || props[4] ) && ( props[2] || props[5] ); } + +double Weight( double v , double start , double end ) +{ + v = ( v - start ) / ( end - start ); + if ( v<0 ) return 1.; + else if( v>1 ) return 0.; + else + { + // P(x) = a x^3 + b x^2 + c x + d + // P (0) = 1 , P (1) = 0 , P'(0) = 0 , P'(1) = 0 + // => d = 1 , a + b + c + d = 0 , c = 0 , 3a + 2b + c = 0 + // => c = 0 , d = 1 , a + b = -1 , 3a + 2b = 0 + // => a = 2 , b = -3 , c = 0 , d = 1 + // => P(x) = 2 x^3 - 3 x^2 + 1 + return 2. * v * v * v - 3. * v * v + 1.; + } +} + +#if defined( _WIN32 ) || defined( _WIN64 ) +double PeakMemoryUsageMB( void ) +{ + HANDLE h = GetCurrentProcess(); + PROCESS_MEMORY_COUNTERS pmc; + return GetProcessMemoryInfo( h , &pmc , sizeof(pmc) ) ? ( (double)pmc.PeakWorkingSetSize )/(1<<20) : 0; +} +#endif // _WIN32 || _WIN64 + + +template< class Real > +struct OctreeProfiler +{ + Octree< Real >& tree; + double t; + + OctreeProfiler( Octree< Real >& t ) : tree(t) { ; } + void start( void ){ t = Time() , tree.resetLocalMemoryUsage(); } + void print( const char* header ) const + { + tree.memoryUsage(); +#if defined( _WIN32 ) || defined( _WIN64 ) + if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() ); + else printf( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() ); +#else // !_WIN32 && !_WIN64 + if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() ); + else printf( "%9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() ); +#endif // _WIN32 || _WIN64 + } + void dumpOutput( const char* header ) const + { + tree.memoryUsage(); +#if defined( _WIN32 ) || defined( _WIN64 ) + if( header ) DumpOutput( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() ); + else DumpOutput( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() ); +#else // !_WIN32 && !_WIN64 + if( header ) DumpOutput( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() ); + else DumpOutput( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() ); +#endif // _WIN32 || _WIN64 + } + void dumpOutput2( std::vector< char* >& comments , const char* header ) const + { + tree.memoryUsage(); +#if defined( _WIN32 ) || defined( _WIN64 ) + if( header ) DumpOutput2( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() ); + else DumpOutput2( comments , "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() ); +#else // !_WIN32 && !_WIN64 + if( header ) DumpOutput2( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() ); + else DumpOutput2( comments , "%9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() ); +#endif // _WIN32 || _WIN64 + } +}; + +template< class Real > +XForm4x4< Real > GetPointXForm( OrientedPointStream< Real >& stream , Real scaleFactor ) +{ + Point3D< Real > min , max; + stream.boundingBox( min , max ); + Point3D< Real > center = ( max + min ) / 2; + Real scale = std::max< Real >( max[0]-min[0] , std::max< Real >( max[1]-min[1] , max[2]-min[2] ) ); + scale *= scaleFactor; + for( int i=0 ; i<3 ; i++ ) center[i] -= scale/2; + XForm4x4< Real > tXForm = XForm4x4< Real >::Identity() , sXForm = XForm4x4< Real >::Identity(); + for( int i=0 ; i<3 ; i++ ) sXForm(i,i) = (Real)(1./scale ) , tXForm(3,i) = -center[i]; + return sXForm * tXForm; +} + +template< class Real , int Degree , BoundaryType BType , class Vertex > +int _Execute( int argc , char* argv[] ) +{ + typedef typename Octree< Real >::template InterpolationInfo< true > InterpolationInfo; + typedef OrientedPointStream< Real > PointStream; + typedef OrientedPointStreamWithData< Real , Point3D< Real > > PointStreamWithData; + typedef TransformedOrientedPointStream< Real > XPointStream; + typedef TransformedOrientedPointStreamWithData< Real , Point3D< Real > > XPointStreamWithData; + Reset< Real >(); + int paramNum = sizeof(params)/sizeof(cmdLineReadable*); + std::vector< char* > comments; + + if( Verbose.set ) echoStdout=1; + + XForm4x4< Real > xForm , iXForm; + if( XForm.set ) + { + FILE* fp = fopen( XForm.value , "r" ); + if( !fp ) + { + fprintf( stderr , "[WARNING] Could not read x-form from: %s\n" , XForm.value ); + xForm = XForm4x4< Real >::Identity(); + } + else + { + for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ ) + { + float f; + if( fscanf( fp , " %f " , &f )!=1 ) fprintf( stderr , "[ERROR] Execute: Failed to read xform\n" ) , exit( 0 ); + xForm(i,j) = (Real)f; + } + fclose( fp ); + } + } + else xForm = XForm4x4< Real >::Identity(); + + DumpOutput2( comments , "Running SSD Reconstruction (Version 9.0)\n" ); + char str[1024]; + for( int i=0 ; iset ) + { + params[i]->writeValue( str ); + if( strlen( str ) ) DumpOutput2( comments , "\t--%s %s\n" , params[i]->name , str ); + else DumpOutput2( comments , "\t--%s\n" , params[i]->name ); + } + + double startTime = Time(); + Real isoValue = 0; + + Octree< Real > tree; + OctreeProfiler< Real > profiler( tree ); + tree.threads = Threads.value; + if( !In.set ) + { + ShowUsage( argv[0] ); + return 0; + } + if( !MaxSolveDepth.set ) MaxSolveDepth.value = Depth.value; + + OctNode< TreeNodeData >::SetAllocator( MEMORY_ALLOCATOR_BLOCK_SIZE ); + + int kernelDepth = KernelDepth.set ? KernelDepth.value : Depth.value-2; + if( kernelDepth>Depth.value ) + { + fprintf( stderr,"[WARNING] %s can't be greater than %s: %d <= %d\n" , KernelDepth.name , Depth.name , KernelDepth.value , Depth.value ); + kernelDepth = Depth.value; + } + + int pointCount; + + Real pointWeightSum; + std::vector< typename Octree< Real >::PointSample >* samples = new std::vector< typename Octree< Real >::PointSample >(); + std::vector< ProjectiveData< Point3D< Real > , Real > >* sampleData = NULL; + SparseNodeData< Real , WEIGHT_DEGREE >* density = NULL; + SparseNodeData< Point3D< Real > , NORMAL_DEGREE >* normalInfo = NULL; + Real targetValue = (Real)0.; + + // Read in the samples (and color data) + { + profiler.start(); + PointStream* pointStream; + char* ext = GetFileExtension( In.value ); + if( Color.set && Color.value>0 ) + { + sampleData = new std::vector< ProjectiveData< Point3D< Real > , Real > >(); + if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStreamWithData< Real , Point3D< Real > , float , Point3D< unsigned char > >( In.value ); + else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYOrientedPointStreamWithData< Real , Point3D< Real > >( In.value , ColorInfo< Real >::PlyProperties , 6 , ColorInfo< Real >::ValidPlyProperties ); + else pointStream = new ASCIIOrientedPointStreamWithData< Real , Point3D< Real > >( In.value , ColorInfo< Real >::ReadASCII ); + } + else + { + if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStream< Real , float >( In.value ); + else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYOrientedPointStream< Real >( In.value ); + else pointStream = new ASCIIOrientedPointStream< Real >( In.value ); + } + delete[] ext; + XPointStream _pointStream( xForm , *pointStream ); + xForm = GetPointXForm( _pointStream , (Real)Scale.value ) * xForm; + if( sampleData ) + { + XPointStreamWithData _pointStream( xForm , ( PointStreamWithData& )*pointStream ); + pointCount = tree.template init< Point3D< Real > >( _pointStream , Depth.value , Confidence.set , *samples , sampleData ); + } + else + { + XPointStream _pointStream( xForm , *pointStream ); + pointCount = tree.template init< Point3D< Real > >( _pointStream , Depth.value , Confidence.set , *samples , sampleData ); + } + iXForm = xForm.inverse(); + delete pointStream; +#pragma omp parallel for num_threads( Threads.value ) + for( int i=0 ; i<(int)samples->size() ; i++ ) (*samples)[i].sample.data.n *= (Real)-1; + + DumpOutput( "Input Points / Samples: %d / %d\n" , pointCount , samples->size() ); + profiler.dumpOutput2( comments , "# Read input into tree:" ); + } + + DenseNodeData< Real , Degree > solution; + // Solve + { + DenseNodeData< Real , Degree > constraints; + InterpolationInfo* iInfo = NULL; + int solveDepth = MaxSolveDepth.value; + + tree.resetNodeIndices(); + + // Get the kernel density estimator + { + profiler.start(); + density = new SparseNodeData< Real , WEIGHT_DEGREE >(); + *density = tree.template setDensityEstimator< WEIGHT_DEGREE >( *samples , kernelDepth , SamplesPerNode.value ); + profiler.dumpOutput2( comments , "# Got kernel density:" ); + } + + // Transform the Hermite samples into a vector field + { + profiler.start(); + normalInfo = new SparseNodeData< Point3D< Real > , NORMAL_DEGREE >(); + *normalInfo = tree.template setNormalField< NORMAL_DEGREE >( *samples , *density , pointWeightSum , BType==BOUNDARY_NEUMANN ); + profiler.dumpOutput2( comments , "# Got normal field:" ); + } + + if( !Density.set ) delete density , density = NULL; + + // Trim the tree and prepare for multigrid + { + profiler.start(); + std::vector< int > indexMap; + + constexpr int MAX_DEGREE = NORMAL_DEGREE > Degree ? NORMAL_DEGREE : Degree; + tree.template inalizeForBroodedMultigrid< MAX_DEGREE , Degree , BType >( FullDepth.value , typename Octree< Real >::template HasNormalDataFunctor< NORMAL_DEGREE >( *normalInfo ) , &indexMap ); + + if( normalInfo ) normalInfo->remapIndices( indexMap ); + if( density ) density->remapIndices( indexMap ); + profiler.dumpOutput2( comments , "# Finalized tree:" ); + } + + // Free up the normal info + if( normalInfo ) delete normalInfo , normalInfo = NULL; + + // Add the interpolation constraints + if( ValueWeight.value>0 || GradientWeight.value>0 ) + { + profiler.start(); + iInfo = new InterpolationInfo( tree , *samples , targetValue , AdaptiveExponent.value , (Real)ValueWeight.value * pointWeightSum , (Real)GradientWeight.value * pointWeightSum ); + constraints = tree.template initDenseNodeData< Degree >( ); + tree.template addInterpolationConstraints< Degree , BType >( *iInfo , constraints , solveDepth ); + profiler.dumpOutput2( comments , "#Set point constraints:" ); + } + + DumpOutput( "Leaf Nodes / Active Nodes / Ghost Nodes: %d / %d / %d\n" , (int)tree.leaves() , (int)tree.nodes() , (int)tree.ghostNodes() ); + DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage())/(1<<20) ); + + // Solve the linear system + { + profiler.start(); + typename Octree< Real >::SolverInfo solverInfo; + solverInfo.cgDepth = CGDepth.value , solverInfo.iters = Iters.value , solverInfo.cgAccuracy = CGSolverAccuracy.value , solverInfo.verbose = Verbose.set , solverInfo.showResidual = ShowResidual.set , solverInfo.lowResIterMultiplier = std::max< double >( 1. , LowResIterMultiplier.value ); + solution = tree.template solveSystem< Degree , BType >( FEMSystemFunctor< Degree , BType >( 0 , 0 , BiLapWeight.value ) , iInfo , constraints , solveDepth , solverInfo ); + profiler.dumpOutput2( comments , "# Linear system solved:" ); + DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) ); + if( iInfo ) delete iInfo , iInfo = NULL; + } + } + + CoredFileMeshData< Vertex > mesh; + + { + profiler.start(); + double valueSum = 0 , weightSum = 0; + typename Octree< Real >::template MultiThreadedEvaluator< Degree , BType > evaluator( &tree , solution , Threads.value ); +#pragma omp parallel for num_threads( Threads.value ) reduction( + : valueSum , weightSum ) + for( int j=0 ; jsize() ; j++ ) + { + ProjectiveData< OrientedPoint3D< Real > , Real >& sample = (*samples)[j].sample; + Real w = sample.weight; + if( w>0 ) weightSum += w , valueSum += evaluator.value( sample.data.p / sample.weight , omp_get_thread_num() , (*samples)[j].node ) * w; + } + isoValue = (Real)( valueSum / weightSum ); + if( !( Color.set && Color.value>0 ) && samples ) delete samples , samples = NULL; + profiler.dumpOutput( "Got average:" ); + DumpOutput( "Iso-Value: %e\n" , isoValue ); + } + + if( VoxelGrid.set ) + { + profiler.start(); + FILE* fp = fopen( VoxelGrid.value , "wb" ); + if( !fp ) fprintf( stderr , "Failed to open voxel file for writing: %s\n" , VoxelGrid.value ); + else + { + int res = 0; + Pointer( Real ) values = tree.template voxelEvaluate< Real , Degree , BType >( solution , res , isoValue , VoxelDepth.value , PrimalVoxel.set ); + fwrite( &res , sizeof(int) , 1 , fp ); + if( sizeof(Real)==sizeof(float) ) fwrite( values , sizeof(float) , res*res*res , fp ); + else + { + float *fValues = new float[res*res*res]; + for( int i=0 ; i , Real > , DATA_DEGREE >* colorData = NULL; + if( sampleData ) + { + colorData = new SparseNodeData< ProjectiveData< Point3D< Real > , Real > , DATA_DEGREE >(); + *colorData = tree.template setDataField< DATA_DEGREE , false >( *samples , *sampleData , (SparseNodeData< Real , WEIGHT_DEGREE >*)NULL ); + delete sampleData , sampleData = NULL; + for( const OctNode< TreeNodeData >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) ) + { + ProjectiveData< Point3D< Real > , Real >* clr = (*colorData)( n ); + if( clr ) (*clr) *= (Real)pow( Color.value , tree.depth( n ) ); + } + } + tree.template getMCIsoSurface< Degree , BType , WEIGHT_DEGREE , DATA_DEGREE >( density , colorData , solution , isoValue , mesh , NonLinearFit.set , !NonManifold.set , PolygonMesh.set ); + DumpOutput( "Vertices / Polygons: %d / %d\n" , mesh.outOfCorePointCount()+mesh.inCorePoints.size() , mesh.polygonCount() ); + if( PolygonMesh.set ) profiler.dumpOutput2( comments , "# Got polygons:" ); + else profiler.dumpOutput2( comments , "# Got triangles:" ); + + if( colorData ) delete colorData , colorData = NULL; + + if( NoComments.set ) + { + if( ASCII.set ) PlyWritePolygons( Out.value , &mesh , PLY_ASCII , NULL , 0 , iXForm ); + else PlyWritePolygons( Out.value , &mesh , PLY_BINARY_NATIVE , NULL , 0 , iXForm ); + } + else + { + if( ASCII.set ) PlyWritePolygons( Out.value , &mesh , PLY_ASCII , &comments[0] , (int)comments.size() , iXForm ); + else PlyWritePolygons( Out.value , &mesh , PLY_BINARY_NATIVE , &comments[0] , (int)comments.size() , iXForm ); + } + } + if( density ) delete density , density = NULL; + DumpOutput2( comments , "# Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-startTime , tree.maxMemoryUsage() ); + + return 1; +} + +#if defined( _WIN32 ) || defined( _WIN64 ) +inline double to_seconds( const FILETIME& ft ) +{ + const double low_to_sec=100e-9; // 100 nanoseconds + const double high_to_sec=low_to_sec*4294967296.0; + return ft.dwLowDateTime*low_to_sec+ft.dwHighDateTime*high_to_sec; +} +#endif // _WIN32 || _WIN64 + +#ifndef FAST_COMPILE +template< class Real , class Vertex > +int Execute( int argc , char* argv[] ) +{ + if( FreeBoundary.set ) + switch( Degree.value ) + { + case 2: return _Execute< Real , 2 , BOUNDARY_FREE , Vertex >( argc , argv ); + case 3: return _Execute< Real , 3 , BOUNDARY_FREE , Vertex >( argc , argv ); + case 4: return _Execute< Real , 4 , BOUNDARY_FREE , Vertex >( argc , argv ); + default: fprintf( stderr , "[ERROR] Only B-Splines of degree 2 - 4 are supported" ) ; return EXIT_FAILURE; + } + else + switch( Degree.value ) + { + case 2: return _Execute< Real , 2 , BOUNDARY_NEUMANN , Vertex >( argc , argv ); + case 3: return _Execute< Real , 3 , BOUNDARY_NEUMANN , Vertex >( argc , argv ); + case 4: return _Execute< Real , 4 , BOUNDARY_NEUMANN , Vertex >( argc , argv ); + default: fprintf( stderr , "[ERROR] Only B-Splines of degree 2 - 4 are supported" ) ; return EXIT_FAILURE; + } +} +#endif // !FAST_COMPILE +int main( int argc , char* argv[] ) +{ +#if defined(WIN32) && defined(MAX_MEMORY_GB) + if( MAX_MEMORY_GB>0 ) + { + SIZE_T peakMemory = 1; + peakMemory <<= 30; + peakMemory *= MAX_MEMORY_GB; + printf( "Limiting memory usage to %.2f GB\n" , float( peakMemory>>30 ) ); + HANDLE h = CreateJobObject( NULL , NULL ); + AssignProcessToJobObject( h , GetCurrentProcess() ); + + JOBOBJECT_EXTENDED_LIMIT_INFORMATION jeli = { 0 }; + jeli.BasicLimitInformation.LimitFlags = JOB_OBJECT_LIMIT_JOB_MEMORY; + jeli.JobMemoryLimit = peakMemory; + if( !SetInformationJobObject( h , JobObjectExtendedLimitInformation , &jeli , sizeof( jeli ) ) ) + fprintf( stderr , "Failed to set memory limit\n" ); + } +#endif // defined(WIN32) && defined(MAX_MEMORY_GB) + double t = Time(); + + cmdLineParse( argc-1 , &argv[1] , sizeof(params)/sizeof(cmdLineReadable*) , params , 1 ); + if( GradientWeight.value<=0 ) fprintf( stderr , "[ERROR] Gradient weight must be positive: %g>=0\n" , GradientWeight.value ) , exit( 0 ); + if( BiLapWeight.value<=0 ) fprintf( stderr , "[ERROR] Bi-Laplacian weight must be positive: %g>=0\n" , BiLapWeight.value ) , exit( 0 ); +#ifdef FAST_COMPILE + static const int Degree = 2; + static const BoundaryType BType = BOUNDARY_NEUMANN; + fprintf( stderr , "[WARNING] Compiling for degree-%d, boundary-%s, single-precision _only_\n" , Degree , BoundaryNames[ BType ] ); + if( Density.set ) + if( Color.set && Color.value>0 ) _Execute< float , Degree , BType , PlyColorAndValueVertex< float > >( argc , argv ); + else _Execute< float , Degree , BType , PlyValueVertex< float > >( argc , argv ); + else + if( Color.set && Color.value>0 ) _Execute< float , Degree , BType , PlyColorVertex< float > >( argc , argv ); + else _Execute< float , Degree , BType , PlyVertex< float > >( argc , argv ); +#else // !FAST_COMPILE + if( Density.set ) + if( Color.set && Color.value>0 ) + if( Double.set ) Execute< double , PlyColorAndValueVertex< float > >( argc , argv ); + else Execute< float , PlyColorAndValueVertex< float > >( argc , argv ); + else + if( Double.set ) Execute< double , PlyValueVertex< float > >( argc , argv ); + else Execute< float , PlyValueVertex< float > >( argc , argv ); + else + if( Color.set && Color.value>0 ) + if( Double.set ) Execute< double , PlyColorVertex< float > >( argc , argv ); + else Execute< float , PlyColorVertex< float > >( argc , argv ); + else + if( Double.set ) Execute< double , PlyVertex< float > >( argc , argv ); + else Execute< float , PlyVertex< float > >( argc , argv ); +#endif // FAST_COMPILE +#if defined( _WIN32 ) || defined( _WIN64 ) + if( Performance.set ) + { + HANDLE cur_thread=GetCurrentThread(); + FILETIME tcreat, texit, tkernel, tuser; + if( GetThreadTimes( cur_thread , &tcreat , &texit , &tkernel , &tuser ) ) + printf( "Time (Wall/User/Kernel): %.2f / %.2f / %.2f\n" , Time()-t , to_seconds( tuser ) , to_seconds( tkernel ) ); + else printf( "Time: %.2f\n" , Time()-t ); + HANDLE h = GetCurrentProcess(); + PROCESS_MEMORY_COUNTERS pmc; + if( GetProcessMemoryInfo( h , &pmc , sizeof(pmc) ) ) printf( "Peak Memory (MB): %d\n" , (int)(pmc.PeakWorkingSetSize>>20) ); + } +#endif // _WIN32 || _WIN64 + return EXIT_SUCCESS; +} diff --git a/src/meshlabplugins/filter_screened_poisson/Src/SparseMatrix.h b/src/meshlabplugins/filter_screened_poisson/Src/SparseMatrix.h index 1dffc8988..f133157ff 100755 --- a/src/meshlabplugins/filter_screened_poisson/Src/SparseMatrix.h +++ b/src/meshlabplugins/filter_screened_poisson/Src/SparseMatrix.h @@ -29,10 +29,10 @@ DAMAGE. #ifndef __SPARSEMATRIX_HPP #define __SPARSEMATRIX_HPP +#define NEW_SPARSE_MATRIX 1 #define ZERO_TESTING_JACOBI 1 -#include "Vector.h" #include "Array.h" template @@ -76,43 +76,26 @@ public: SparseMatrix operator * (const T& V) const; SparseMatrix& operator *= (const T& V); - - template - Vector operator * (const Vector& V) const; - template - Vector Multiply( const Vector& V ) const; - template - void Multiply( const Vector& In , Vector& Out , int threads=1 ) const; - - static int Solve (const SparseMatrix& M,const Vector& b, int iters,Vector& solution,const T eps=1e-8); - - template - static int SolveSymmetric( const SparseMatrix& M , const Vector& b , int iters , Vector& solution , const T2 eps=1e-8 , int reset=1 , int threads=1 ); + template< class T2 > void Multiply( ConstPointer( T2 ) in , Pointer( T2 ) out , int threads=1 ) const; + template< class T2 > void MultiplyAndAddAverage( ConstPointer( T2 ) in , Pointer( T2 ) out , int threads=1 ) const; bool write( FILE* fp ) const; bool write( const char* fileName ) const; bool read( FILE* fp ); bool read( const char* fileName ); - template< class T2 > - static int SolveJacobi( const SparseMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , Vector& Mx , T2 sor , int threads=1 , int offset=0 ); - template< class T2 > - static int SolveGS( const SparseMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , bool forward , int offset=0 ); - template< class T2 > - static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , bool forward , int threads=1 , int offset=0 ); - - template< class T2 > - static int SolveJacobi( const SparseMatrix& M , const Vector& b , Vector& x , Vector& Mx , T2 sor , int threads=1 , int offset=0 ); - template< class T2 > - static int SolveGS( const SparseMatrix& M , const Vector& b , Vector& x , bool forward , int offset=0 ); - template< class T2 > - static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix& M , const Vector& b , Vector& x , bool forward , int threads=1 , int offset=0 ); - - template< class T2 > - void getDiagonal( Vector< T2 >& diagonal , int threads=1 , int offset=0 ) const; + template< class T2 > void getDiagonal( Pointer( T2 ) diagonal , int threads=1 ) const; + template< class T2 > static int SolveJacobi( const SparseMatrix& M , ConstPointer( T2 ) b , Pointer( T2 ) x , Pointer( T2 ) Mx , T2 sor , int threads=1 ); + template< class T2 > static int SolveJacobi( const SparseMatrix& M , ConstPointer( T2 ) diagonal , ConstPointer( T2 ) b , Pointer( T2 ) x , Pointer( T2 ) Mx , T2 sor , int threads=1 ); + template< class T2 > static int SolveGS( const SparseMatrix& M , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward ); + template< class T2 > static int SolveGS( const SparseMatrix& M , ConstPointer( T2 ) diagonal , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward ); + template< class T2 > static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix& M , ConstPointer( T2 ) diagonal , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward , int threads=1 ); + template< class T2 > static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix& M , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward , int threads=1 ); + template< class T2 > static int SolveCG( const SparseMatrix& M , ConstPointer( T2 ) b , int iters , Pointer( T2 ) x , T2 eps=1e-8 , int reset=1 , bool addDCTerm=false , bool solveNormal=false , int threads=1 ); }; +#if !NEW_SPARSE_MATRIX template< class T2 > struct MapReduceVector { @@ -203,6 +186,7 @@ public: template< class T2 > void getDiagonal( Vector< T2 >& diagonal , int threads=1 ) const; }; +#endif // !NEW_SPARSE_MATRIX #include "SparseMatrix.inl"