35 #ifndef AIMS_MOMENT_MOMINVARIANT_H 36 #define AIMS_MOMENT_MOMINVARIANT_H 48 std::ostream& operator << ( std::ostream&, const MomentInvariant< T >& );
78 const double *
invariant()
const {
return &_inv[ 0 ]; }
96 std::complex< double > _n0;
97 std::complex< double > _n1[ 3 ];
98 std::complex< double > _n2[ 5 ];
99 std::complex< double > _n3[ 7 ];
116 std::complex< double > _p1;
117 std::complex< double > _p2;
118 std::complex< double > _p3;
119 std::complex< double > _p4;
120 std::complex< double > _p5;
121 std::complex< double > _p6;
122 std::complex< double > _p7;
123 std::complex< double > _p8;
124 std::complex< double > _q1;
125 std::complex< double > _q2;
126 std::complex< double > _q4;
127 std::complex< double > _q5;
128 std::complex< double > _r1;
129 std::complex< double > _r2;
131 std::complex< double > _de;
135 template<
class T >
inline 143 template<
class T >
inline 150 for (
int i=0; i<12; i++ ) _inv[ i ] = other.
invariant()[ i ];
154 template<
class T >
inline 158 _c00 = 2.0 * sqrt(
M_PI ) / 3.0;
160 _c10 = 2.0 * sqrt( 3.0 *
M_PI ) / 5.0;
161 _c11 = sqrt( 6.0 *
M_PI ) / 5.0;
163 _c20 = 2.0 * sqrt(
M_PI / 5.0 ) / 3.0;
164 _c22 = sqrt( 2.0 *
M_PI / 15.0 );
166 _c30 = 2.0 * sqrt(
M_PI / 7.0 ) / 5.0;
167 _c31 = sqrt( 3.0 *
M_PI / 7.0 ) / 5.0;
168 _c32 = sqrt( 6.0 *
M_PI / 35.0 );
169 _c33 = sqrt(
M_PI / 35.0 );
172 _cc22 = 16.0 *
M_PI / 45.0 / sqrt( 5.0 );
173 _cc222 = 32.0 *
M_PI * sqrt( 2.0 *
M_PI / 7.0 ) / 675.0;
174 _cc11 = -4.0 *
M_PI * sqrt( 3.0 ) / 25.0;
175 _cc33 = -16.0 *
M_PI / 175.0 / sqrt( 7.0 );
177 _ce = 1.0 / sqrt( 5.0 );
179 _p1 = std::complex< double >( sqrt( 10.0 / 21.0 ), 0.0 );
180 _p2 = std::complex< double >( -2.0 * sqrt( 5.0 / 21.0 ), 0.0 );
181 _p3 = std::complex< double >( sqrt( 2.0 / 7.0 ), 0.0 );
182 _p4 = std::complex< double >( 5.0 / sqrt( 21.0 ), 0.0 );
183 _p5 = std::complex< double >( -sqrt( 5.0 / 7.0 ), 0.0 );
184 _p6 = std::complex< double >( sqrt( 2.0 / 21.0 ), 0.0 );
185 _p7 = std::complex< double >( -sqrt( 3.0 / 7.0 ), 0.0 );
186 _p8 = std::complex< double >( 2.0 / sqrt( 21.0 ), 0.0 );
188 _q1 = std::complex< double >( -sqrt( 5.0 / 21.0 ), 0.0 );
189 _q2 = std::complex< double >( 1.0 / sqrt( 21.0 ), 0.0 );
190 _q4 = std::complex< double >( 1.0 / sqrt( 7.0 ), 0.0 );
191 _q5 = std::complex< double >( -2.0 * sqrt( 2.0 / 21.0 ), 0.0 );
193 _r1 = std::complex< double >( sqrt( 2.0 ), 0.0 );
194 _r2 = std::complex< double >( sqrt( 2.0 / 3.0 ), 0.0 );
196 _de = std::complex< double >( 2.0, 0.0 );
200 template<
class T >
inline 205 _mom->update( x, y, z, dir );
212 template<
class T >
inline 217 _mom->doit( d, label, dir );
224 template<
class T >
inline 229 _mom->doit( surfTri );
236 template<
class T >
inline 248 template<
class T >
inline 251 _n0 = std::complex< double >( _c00 * ( _mom->m2()[ 0 ] + _mom->m2()[ 1 ] +
252 _mom->m2()[ 2 ] ) , 0.0 );
254 _n1[ 0 ] = std::complex< double >( _c11 * ( _mom->m3()[ 0 ] + _mom->m3()[ 5 ] +
256 -_c11 * ( _mom->m3()[ 1 ] + _mom->m3()[ 3 ] +
258 _n1[ 1 ] = std::complex< double >( _c10 * ( _mom->m3()[ 2 ] + _mom->m3()[ 4 ] +
259 _mom->m3()[ 6 ] ), 0.0 );
260 _n1[ 2 ] = std::complex< double >( -_n1[ 0 ].real(), _n1[ 0 ].imag() );
262 _n2[ 0 ] = std::complex< double >( _c22 * ( _mom->m2()[ 0 ] - _mom->m2()[ 1 ] ),
263 -2.0 * _c22 * _mom->m2()[ 3 ] );
264 _n2[ 1 ] = std::complex< double >( 2.0 * _c22 * _mom->m2()[ 4 ],
265 -2.0 * _c22 * _mom->m2()[ 5 ] );
266 _n2[ 2 ] = std::complex< double >( _c20 * ( 2.0 * _mom->m2()[ 2 ] -
267 _mom->m2()[ 0 ] - _mom->m2()[ 1 ] ),
269 _n2[ 3 ] = std::complex< double >( -_n2[ 1 ].real(), _n2[ 1 ].imag() );
270 _n2[ 4 ] = conj( _n2[ 0 ] );
272 _n3[ 0 ] = std::complex< double >( _c33 * ( _mom->m3()[ 0 ] -
273 3.0 * _mom->m3()[ 5 ] ),
274 _c33 * ( _mom->m3()[ 1 ] -
275 3.0 * _mom->m3()[ 3 ] ) );
276 _n3[ 1 ] = std::complex< double >( _c32 * ( _mom->m3()[ 4 ] - _mom->m3()[ 6 ] ),
277 -2.0 * _c32 * _mom->m3()[ 9 ] );
278 _n3[ 2 ] = std::complex< double >( -_c31 * ( _mom->m3()[ 0 ] + _mom->m3()[ 5 ] -
279 4.0 * _mom->m3()[ 7 ] ),
280 _c31 * ( _mom->m3()[ 1 ] + _mom->m3()[ 3 ] -
281 4.0 * _mom->m3()[ 8 ] ) );
282 _n3[ 3 ] = std::complex< double >( _c30 * ( 2.0 * _mom->m3()[ 2 ] -
283 3.0 * ( _mom->m3()[ 4 ] +
284 _mom->m3()[ 6 ] ) ), 0.0 );
285 _n3[ 4 ] = std::complex< double >( -_n3[ 2 ].real(), _n3[ 2 ].imag() );
286 _n3[ 5 ] = conj( _n3[ 1 ] );
287 _n3[ 6 ] = std::complex< double >( -_n3[ 0 ].real(), _n3[ 0 ].imag() );
291 template<
class T >
inline 294 double mo200 = _mom->m2()[ 0 ];
295 double mo020 = _mom->m2()[ 1 ];
296 double mo002 = _mom->m2()[ 2 ];
297 double mo110 = _mom->m2()[ 3 ];
298 double mo101 = _mom->m2()[ 4 ];
299 double mo011 = _mom->m2()[ 5 ];
301 double mo300 = _mom->m3()[ 0 ];
302 double mo030 = _mom->m3()[ 1 ];
303 double mo003 = _mom->m3()[ 2 ];
304 double mo210 = _mom->m3()[ 3 ];
305 double mo201 = _mom->m3()[ 4 ];
306 double mo120 = _mom->m3()[ 5 ];
307 double mo021 = _mom->m3()[ 6 ];
308 double mo102 = _mom->m3()[ 7 ];
309 double mo012 = _mom->m3()[ 8 ];
310 double mo111 = _mom->m3()[ 9 ];
312 double m200_2 = mo200 * mo200;
313 double m020_2 = mo020 * mo020;
314 double m002_2 = mo002 * mo002;
315 double m110_2 = mo110 * mo110;
316 double m101_2 = mo101 * mo101;
317 double m011_2 = mo011 * mo011;
319 double m200m020 = mo200 * mo020;
320 double m200m002 = mo200 * mo002;
321 double m020m002 = mo020 * mo002;
323 double m200_3 = m200_2 * mo200;
324 double m020_3 = m020_2 * mo020;
325 double m002_3 = m002_2 * mo002;
326 double m200_2m020 = m200_2 * mo020;
327 double m200_2m002 = m200_2 * mo002;
328 double m020_2m200 = m020_2 * mo200;
329 double m020_2m002 = m020_2 * mo002;
330 double m002_2m200 = m002_2 * mo200;
331 double m002_2m020 = m002_2 * mo020;
332 double m200m110_2 = mo200 * m110_2;
333 double m200m101_2 = mo200 * m101_2;
334 double m200m011_2 = mo200 * m011_2;
335 double m020m110_2 = mo020 * m110_2;
336 double m020m101_2 = mo020 * m101_2;
337 double m020m011_2 = mo020 * m011_2;
338 double m002m110_2 = mo002 * m110_2;
339 double m002m101_2 = mo002 * m101_2;
340 double m002m011_2 = mo002 * m011_2;
341 double m200m020m002 = m200m020 * mo002;
342 double m110m101m011 = mo110 * mo101 * mo011;
344 double m300_2 = mo300 * mo300;
345 double m030_2 = mo030 * mo030;
346 double m003_2 = mo003 * mo003;
347 double m210_2 = mo210 * mo210;
348 double m201_2 = mo201 * mo201;
349 double m120_2 = mo120 * mo120;
350 double m021_2 = mo021 * mo021;
351 double m102_2 = mo102 * mo102;
352 double m012_2 = mo012 * mo012;
353 double m111_2 = mo111 * mo111;
355 double m300m120 = mo300 * mo120;
356 double m300m102 = mo300 * mo102;
357 double m030m210 = mo030 * mo210;
358 double m030m012 = mo030 * mo012;
359 double m003m201 = mo003 * mo201;
360 double m003m021 = mo003 * mo021;
361 double m120m102 = mo120 * mo102;
362 double m210m012 = mo210 * mo012;
363 double m201m021 = mo201 * mo021;
365 double m0 = _mom->m0();
366 double d1 = pow( m0, 5.0 / 3.0 );
367 double d2 = pow( m0, 10.0 / 3.0 );
368 double d3 = pow( m0, 5.0 );
369 double d4 = pow( m0, 4.0 );
370 double d5 = pow( m0, 17.0 / 3.0 );
371 double d6 = pow( m0, 8.0 );
374 double I0 = _c00 * ( mo200 + mo020 + mo002 ) / d1;
375 double aI0 = fabs( I0 );
376 _inv[ 0 ] = I0 * sqrt( aI0 ) / aI0;
378 double I1 = _cc22 * ( m200_2 + m020_2 + m002_2 +
379 3.0 * ( m110_2 + m101_2 + m011_2 ) -
380 m200m020 - m200m002 - m020m002 ) / d2;
381 double aI1 = fabs( I1 );
382 _inv[ 1 ] = I1 * pow( aI1, 0.25 ) / aI1;
384 double I2 = _cc222 * ( -2.0 * ( m200_3 + m020_3 + m002_3 ) +
385 3 * ( m200_2m002 + m200_2m020 + m020_2m200 +
386 m020_2m002 + m002_2m200 + m002_2m020 ) -
387 9.0 * ( m200m110_2 + m200m101_2 + m020m110_2 +
388 m020m011_2 + m002m101_2 + m002m011_2 ) +
389 18.0 * ( m200m011_2 + m020m101_2 + m002m110_2 ) -
390 12.0 * m200m020m002 - 54.0 * m110m101m011 ) / d3;
391 double aI2 = fabs( I2 );
392 _inv[ 2 ] = I2 * pow( aI2, 1.0 / 6.0 ) / aI2;
395 double I3 = _cc11 * ( m300_2 + m030_2 + m003_2 + m210_2 + m201_2 + m120_2 +
396 m021_2 + m102_2 + m012_2 +
397 2.0 * ( m300m120 + m300m102 + m030m210 + m030m012 +
398 m003m201 + m003m021 + m210m012 + m201m021 +
400 double aI3 = fabs( I3 );
401 _inv[ 3 ] = I3 * pow( aI3, 1.0 / 6.0 ) / aI3;
403 double I4 = _cc33 * ( m300_2 + m030_2 + m003_2 +
404 6.0 * ( m210_2 + m201_2 + m120_2 + m021_2 + m102_2 +
405 m012_2 ) + 15.0 * m111_2 -
406 3.0 * ( m300m120 + m300m102 + m030m210 + m030m012 +
407 m003m201 + m003m021 + m210m012 + m201m021 +
409 double aI4 = fabs( I4 );
410 _inv[ 4 ] = I4 * pow( aI4, 1.0 / 6.0 ) / aI4;
412 std::complex< double > zeta[ 5 ];
413 zeta[ 0 ] = _p1 * _n3[ 0 ] * _n3[ 4 ] + _p2 * _n3[ 1 ] * _n3[ 3 ] +
414 _p3 * _n3[ 2 ] * _n3[ 2 ];
415 zeta[ 1 ] = _p4 * _n3[ 0 ] * _n3[ 5 ] + _p5 * _n3[ 1 ] * _n3[ 4 ] +
416 _p6 * _n3[ 2 ] * _n3[ 3 ];
417 zeta[ 2 ] = _p4 * _n3[ 0 ] * _n3[ 6 ] + _p7 * _n3[ 2 ] * _n3[ 4 ] +
418 _p8 * _n3[ 3 ] * _n3[ 3 ];
419 zeta[ 3 ] = _p4 * _n3[ 1 ] * _n3[ 6 ] + _p5 * _n3[ 2 ] * _n3[ 5 ] +
420 _p6 * _n3[ 3 ] * _n3[ 4 ];
421 zeta[ 4 ] = _p1 * _n3[ 2 ] * _n3[ 6 ] + _p2 * _n3[ 3 ] * _n3[ 5 ] +
422 _p3 * _n3[ 4 ] * _n3[ 4 ];
424 std::complex< double > xi[ 5 ];
425 xi[ 0 ] = _q2 * _n3[ 2 ] * _n1[ 0 ] + _q1 * _n3[ 1 ] * _n1[ 1 ] -
426 _p5 * _n3[ 0 ] * _n1[ 2 ];
427 xi[ 1 ] = _q4 * _n3[ 3 ] * _n1[ 0 ] + _q5 * _n3[ 2 ] * _n1[ 1 ] +
428 _p1 * _n3[ 1 ] * _n1[ 2 ];
429 xi[ 2 ] = _p3 * _n3[ 4 ] * _n1[ 0 ] + _p7 * _n3[ 3 ] * _n1[ 1 ] +
430 _p3 * _n3[ 2 ] * _n1[ 2 ];
431 xi[ 3 ] = _p1 * _n3[ 5 ] * _n1[ 0 ] + _q5 * _n3[ 4 ] * _n1[ 1 ] +
432 _q4 * _n3[ 3 ] * _n1[ 2 ];
433 xi[ 4 ] = -_p5 * _n3[ 6 ] * _n1[ 0 ] + _q1 * _n3[ 5 ] * _n1[ 1 ] +
434 _q2 * _n3[ 4 ] * _n1[ 2 ];
436 std::complex< double > nu[ 5 ];
437 nu[ 0 ] = _n1[ 0 ] * _n1[ 0 ];
438 nu[ 1 ] = _r1 * _n1[ 0 ] * _n1[ 1 ];
439 nu[ 2 ] = _r2 * ( _n1[ 0 ] * _n1[ 2 ] + _n1[ 1 ] * _n1[ 1 ] );
440 nu[ 3 ] = _r1 * _n1[ 1 ] * _n1[ 2 ];
441 nu[ 4 ] = _n1[ 2 ] * _n1[ 2 ];
443 std::complex< double > dI5 = xi[ 4 ] * nu[ 0 ] - xi[ 3 ] * nu[ 1 ] +
444 xi[ 2 ] * nu[ 2 ] - xi[ 1 ] * nu[ 3 ] +
446 double I5 = _ce * dI5.real() / d6;
447 double aI5 = fabs( I5 );
448 _inv[ 5 ] = I5 * pow( aI5, 1.0 / 12.0 ) / aI5;
450 std::complex< double > dI6 = _de * xi[ 0 ] * xi[ 4 ] - _de * xi[ 1 ] * xi[ 3 ] +
452 double I6 = _ce * dI6.real() / d6;
453 double aI6 = fabs( I6 );
454 _inv[ 6 ] = I6 * pow( aI6, 1.0 / 12.0 ) / aI6;
456 std::complex< double > dI7 = zeta[ 4 ] * xi[ 0 ] - zeta[ 3 ] * xi[ 1 ] +
457 zeta[ 2 ] * xi[ 2 ] - zeta[ 1 ] * xi[ 3 ] +
459 double I7 = _ce * dI7.real() / d6;
460 double aI7 = fabs( I7 );
461 _inv[ 7 ] = I7 * pow( aI7, 1.0 / 12.0 ) / aI7;
463 std::complex< double > dI8 = _de * zeta[ 0 ] * zeta[ 4 ] -
464 _de * zeta[ 1 ] * zeta[ 3 ] + zeta[ 2 ] * zeta[ 2 ];
465 double I8 = _ce * dI8.real() / d6;
466 double aI8 = fabs( I8 );
467 _inv[ 8 ] = I8 * pow( aI8, 1.0 / 12.0 ) / aI8;
469 std::complex< double > dI9 = nu[ 4 ] * _n2[ 0 ] - nu[ 3 ] * _n2[ 1 ] +
470 nu[ 2 ] * _n2[ 2 ] - nu[ 1 ] * _n2[ 3 ] +
472 double I9 = _ce * dI9.real() / d5;
473 double aI9 = fabs( I9 );
474 _inv[ 9 ] = I9 * pow( aI9, 0.125 ) / aI9;
476 std::complex< double > dI10 = xi[ 4 ] * _n2[ 0 ] - xi[ 3 ] * _n2[ 1 ] +
477 xi[ 2 ] * _n2[ 2 ] - xi[ 1 ] * _n2[ 3 ] +
479 double I10 = _ce * dI10.real() / d5;
480 double aI10 = fabs( I10 );
481 _inv[ 10 ] = I10 * pow( aI10, 0.125 ) / aI10;
483 std::complex< double > dI11 = zeta[ 4 ] * _n2[ 0 ] - zeta[ 3 ] * _n2[ 1 ] +
484 zeta[ 2 ] * _n2[ 2 ] - zeta[ 1 ] * _n2[ 3 ] +
485 zeta[ 0 ] * _n2[ 4 ];
486 double I11 = _ce * dI11.real() / d5;
487 double aI11 = fabs( I11 );
488 _inv[ 11 ] = I11 * pow( aI11, 0.125 ) / aI11;
492 template<
class T >
inline 496 for (
int i=0; i<12; i++ ) out << mom.
invariant()[ i ] <<
" ";
virtual ~MomentInvariant()
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle
void doit(AimsData< T > &, T, int dir=(int) Moment< T >::mAdd)
const double * invariant() const
Moment< T > & moment() const
MomentInvariant(Moment< T > *m=0)
std::ostream & operator<<(std::ostream &, const MomentInvariant< T > &)
void update(double, double, double, int dir=(int) Moment< T >::mAdd)