29 #include "poisson_exceptions.h"
30 #include "binary_node.h"
59 template<
int Degree >
inline bool LeftOverlap(
unsigned int depth ,
int offset )
62 if( Degree & 1 )
return (offset < 1+Degree) && (offset > -1-Degree );
63 else return (offset < Degree) && (offset > -2-Degree );
65 template<
int Degree >
inline bool RightOverlap(
unsigned int depth ,
int offset )
69 if( Degree & 1 )
return (offset > 2-1-Degree) && (offset < 2+1+Degree );
70 else return (offset > 2-2-Degree) && (offset < 2+ Degree );
72 template<
int Degree >
inline int ReflectLeft(
unsigned int depth ,
int offset )
74 if( Degree&1 )
return -offset;
75 else return -1-offset;
77 template<
int Degree >
inline int ReflectRight(
unsigned int depth ,
int offset )
80 if( Degree&1 )
return r -offset;
81 else return r-1-offset;
84 template<
int Degree ,
class Real >
87 vvDotTable = dvDotTable = ddDotTable = NULL;
88 valueTables = dValueTables = NULL;
91 functionCount = sampleCount = 0;
94 template<
int Degree ,
class Real >
103 delete[] valueTables;
104 delete[] dValueTables;
106 delete[] baseFunctions;
107 delete[] baseBSplines;
109 vvDotTable = dvDotTable = ddDotTable = NULL;
110 valueTables = dValueTables=NULL;
111 baseFunctions = NULL;
116 template<
int Degree,
class Real>
119 this->useDotRatios = useDotRatios;
120 this->reflectBoundary = reflectBoundary;
130 for(
int i=0 ; i<=Degree ; i++ ) baseBSpline[i] = Polynomial< Degree >::BSplineComponent( i ).shift(
double(-(Degree+1)/2) + i - 0.5 );
131 dBaseFunction = baseFunction.derivative();
134 for(
int i=0 ; i<Degree+3 ; i++ )
136 sPolys[i].
start = double(-(Degree+1)/2) + i - 1.5;
138 if( i<=Degree ) sPolys[i].
p += baseBSpline[i ].shift( -1 );
139 if( i>=1 && i<=Degree+1 ) sPolys[i].
p += baseBSpline[i-1];
140 for(
int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
142 leftBaseFunction.set( sPolys , Degree+3 );
143 for(
int i=0 ; i<Degree+3 ; i++ )
145 sPolys[i].
start = double(-(Degree+1)/2) + i - 0.5;
147 if( i<=Degree ) sPolys[i].
p += baseBSpline[i ];
148 if( i>=1 && i<=Degree+1 ) sPolys[i].
p += baseBSpline[i-1].shift( 1 );
149 for(
int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
151 rightBaseFunction.set( sPolys , Degree+3 );
152 dLeftBaseFunction = leftBaseFunction.derivative();
153 dRightBaseFunction = rightBaseFunction.derivative();
154 leftBSpline = rightBSpline = baseBSpline;
155 leftBSpline [1] += leftBSpline[2].
shift( -1 ) , leftBSpline[0] *= 0;
156 rightBSpline[1] += rightBSpline[0].
shift( 1 ) , rightBSpline[2] *= 0;
158 for(
int i=0 ; i<functionCount ; i++ )
161 baseFunctions[i] = baseFunction.scale(w).shift(c);
162 baseBSplines[i] = baseBSpline.scale(w).shift(c);
163 if( reflectBoundary )
168 if ( off==0 ) baseFunctions[i] = leftBaseFunction.scale(w).shift(c);
169 else if( off==r-1 ) baseFunctions[i] = rightBaseFunction.scale(w).shift(c);
170 if ( off==0 ) baseBSplines[i] = leftBSpline.scale(w).shift(c);
171 else if( off==r-1 ) baseBSplines[i] = rightBSpline.scale(w).shift(c);
175 template<
int Degree,
class Real>
178 clearDotTables( flags );
179 int size = ( functionCount*functionCount + functionCount )>>1;
180 int fullSize = functionCount*functionCount;
181 if( flags & VV_DOT_FLAG )
183 vvDotTable =
new Real[size];
184 memset( vvDotTable , 0 ,
sizeof(
Real)*size );
186 if( flags & DV_DOT_FLAG )
188 dvDotTable =
new Real[fullSize];
189 memset( dvDotTable , 0 ,
sizeof(
Real)*fullSize );
191 if( flags & DD_DOT_FLAG )
193 ddDotTable =
new Real[size];
194 memset( ddDotTable , 0 ,
sizeof(
Real)*size );
196 double vvIntegrals[Degree+1][Degree+1];
197 double vdIntegrals[Degree+1][Degree ];
198 double dvIntegrals[Degree ][Degree+1];
199 double ddIntegrals[Degree ][Degree ];
200 int vvSums[Degree+1][Degree+1];
201 int vdSums[Degree+1][Degree ];
202 int dvSums[Degree ][Degree+1];
203 int ddSums[Degree ][Degree ];
204 SetBSplineElementIntegrals< Degree , Degree >( vvIntegrals );
209 for(
int d1=0 ; d1<=depth ; d1++ )
210 for(
int off1=0 ; off1<(1<<d1) ; off1++ )
220 for(
int i=0 ; i<int(b1.size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
222 if( b1[i][j] && start1==-1 ) start1 = i;
223 if( b1[i][j] ) end1 = i+1;
225 for(
int d2=d1 ; d2<=depth ; d2++ )
227 for(
int off2=0 ; off2<(1<<d2) ; off2++ )
229 int start2 = off2-Degree;
230 int end2 = off2+Degree+1;
231 if( start2>=end1 || start1>=end2 )
continue;
232 start2 = std::max< int >( start1 , start2 );
233 end2 = std::min< int >( end1 , end2 );
234 if( d1==d2 && off2<off1 )
continue;
240 int idx = SymmetricIndex( ii , jj );
241 int idx1 = Index( ii , jj ) , idx2 = Index( jj , ii );
243 memset( vvSums , 0 ,
sizeof(
int ) * ( Degree+1 ) * ( Degree+1 ) );
244 memset( vdSums , 0 ,
sizeof(
int ) * ( Degree+1 ) * ( Degree ) );
245 memset( dvSums , 0 ,
sizeof(
int ) * ( Degree ) * ( Degree+1 ) );
246 memset( ddSums , 0 ,
sizeof(
int ) * ( Degree ) * ( Degree ) );
247 for(
int i=start2 ; i<end2 ; i++ )
249 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) vvSums[j][k] += b1[i][j] * b2[i][k];
250 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) vdSums[j][k] += b1[i][j] * db2[i][k];
251 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) dvSums[j][k] += db1[i][j] * b2[i][k];
252 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) ddSums[j][k] += db1[i][j] * db2[i][k];
254 double vvDot = 0 , dvDot = 0 , vdDot = 0 , ddDot = 0;
255 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) vvDot += vvIntegrals[j][k] * vvSums[j][k];
256 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) vdDot += vdIntegrals[j][k] * vdSums[j][k];
257 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) dvDot += dvIntegrals[j][k] * dvSums[j][k];
258 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) ddDot += ddIntegrals[j][k] * ddSums[j][k];
265 if( fabs(vvDot)<1e-15 )
continue;
266 if( flags & VV_DOT_FLAG ) vvDotTable [idx] =
Real( vvDot );
269 if( flags & DV_DOT_FLAG ) dvDotTable[idx1] =
Real( dvDot / vvDot );
270 if( flags & DV_DOT_FLAG ) dvDotTable[idx2] =
Real( vdDot / vvDot );
271 if( flags & DD_DOT_FLAG ) ddDotTable[idx ] =
Real( ddDot / vvDot );
275 if( flags & DV_DOT_FLAG ) dvDotTable[idx1] =
Real( dvDot );
276 if( flags & DV_DOT_FLAG ) dvDotTable[idx2] =
Real( dvDot );
277 if( flags & DD_DOT_FLAG ) ddDotTable[idx ] =
Real( ddDot );
285 for(
int i=0 ; i<int(b1.size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
287 if( b1[i][j] && start1==-1 ) start1 = i;
288 if( b1[i][j] ) end1 = i+1;
293 template<
int Degree,
class Real>
296 if (flags & VV_DOT_FLAG) {
297 delete[] vvDotTable ; vvDotTable = NULL;
298 delete[] dvDotTable ; dvDotTable = NULL;
299 delete[] ddDotTable ; ddDotTable = NULL;
302 template<
int Degree ,
class Real >
308 double _start = ( off + 0.5 - 0.5*(Degree+1) ) / res - smooth;
309 double _end = ( off + 0.5 + 0.5*(Degree+1) ) / res + smooth;
313 start = int( floor( _start * (sampleCount-1) + 1 ) );
314 if( start<0 ) start = 0;
318 end = int( ceil( _end * (sampleCount-1) - 1 ) );
319 if( end>=sampleCount ) end = sampleCount-1;
321 template<
int Degree,
class Real>
325 if( flags & VALUE_FLAG ) valueTables =
new Real[functionCount*sampleCount];
326 if( flags & D_VALUE_FLAG ) dValueTables =
new Real[functionCount*sampleCount];
329 for(
int i=0 ; i<functionCount ; i++ )
338 function = baseFunctions[i];
341 for(
int j=0 ; j<sampleCount ; j++ )
343 double x=double(j)/(sampleCount-1);
344 if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=
Real(
function(x));}
345 if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=
Real(dFunction(x));}
349 template<
int Degree,
class Real>
352 if(flags & VALUE_FLAG){ valueTables=
new Real[functionCount*sampleCount];}
353 if(flags & D_VALUE_FLAG){dValueTables=
new Real[functionCount*sampleCount];}
356 for(
int i=0;i<functionCount;i++){
357 if(valueSmooth>0) {
function=baseFunctions[i].
MovingAverage(valueSmooth);}
358 else {
function=baseFunctions[i];}
360 else {dFunction=baseFunctions[i].
derivative();}
362 for(
int j=0;j<sampleCount;j++){
363 double x=double(j)/(sampleCount-1);
364 if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=
Real(
function(x));}
365 if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=
Real(dFunction(x));}
371 template<
int Degree,
class Real>
373 delete[] valueTables;
374 delete[] dValueTables;
375 valueTables=dValueTables=NULL;
378 template<
int Degree,
class Real>
380 template<
int Degree,
class Real>
383 if( i1>i2 )
return ((i1*i1+i1)>>1)+i2;
384 else return ((i2*i2+i2)>>1)+i1;
386 template<
int Degree,
class Real>
391 index = ((i2*i2+i2)>>1)+i1;
396 index = ((i1*i1+i1)>>1)+i2;
405 template<
int Degree >
411 for(
int i=0 ; i<=Degree ; i++ )
413 int idx = -_off + offset + i;
414 if( idx>=0 && idx<res ) (*this)[idx][i] = 1;
418 _addLeft( offset-2*res , boundary ) , _addRight( offset+2*res , boundary );
419 if( Degree&1 ) _addLeft( offset-res , boundary ) , _addRight( offset+res , boundary );
420 else _addLeft( -offset-1 , boundary ) , _addRight( -offset-1+2*res , boundary );
423 template<
int Degree >
426 int res = int( this->size() );
428 for(
int i=0 ; i<=Degree ; i++ )
430 int idx = -_off + offset + i;
431 if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set =
true;
433 if( set ) _addLeft( offset-2*res , boundary );
435 template<
int Degree >
438 int res = int( this->size() );
440 for(
int i=0 ; i<=Degree ; i++ )
442 int idx = -_off + offset + i;
443 if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set =
true;
445 if( set ) _addRight( offset+2*res , boundary );
447 template<
int Degree >
458 template<
int Degree >
461 d.resize( this->size() );
463 for(
int i=0 ; i<int(this->size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
465 if( j-1>=0 ) d[i][j-1] -= (*this)[i][j];
466 if( j<Degree ) d[i][j ] += (*this)[i][j];
472 template<
int Degree1 ,
int Degree2 >
475 for(
int i=0 ; i<=Degree1 ; i++ )
478 for(
int j=0 ; j<=Degree2 ; j++ )
481 integrals[i][j] = ( p1 * p2 ).integral( 0 , 1 );