30 #define __IMSimplex_H 1
34 #include <boost/assign/list_of.hpp>
35 #include <boost/assign/std/vector.hpp>
36 #include <boost/numeric/ublas/matrix.hpp>
37 #include <boost/numeric/ublas/matrix_proxy.hpp>
38 #include <boost/numeric/ublas/vector.hpp>
39 #include <boost/numeric/ublas/io.hpp>
41 #include <feel/feelcore/feel.hpp>
47 namespace ublas = boost::numeric::ublas;
56 using namespace boost::assign;
62 template<
int Order,
typename T>
68 struct IMTriangle<1,T>
71 static const uint16_type nDim = 2;
72 static const uint16_type nOrder = 1;
73 static const uint16_type nPoints = 1;
81 push_back( q[0] ).repeat( 3, value_type( 1 )/value_type( 3 ) );
85 std::vector<uint16_type> m;
86 std::vector<value_type> w;
87 std::vector<std::vector<value_type> > q;
90 template<
typename T>
struct IMTriangle<0,T>:
public IMTriangle<1,T> {};
93 struct IMTriangle<2,T>
96 static const uint16_type nDim = 2;
97 static const uint16_type nOrder = 2;
98 static const uint16_type nPoints = 3;
105 w += value_type( 1 )/value_type( 3 );
106 q[0] += value_type( 2 )/value_type( 3 ), repeat( 2, value_type( 1 )/value_type( 6 ) );
110 std::vector<uint16_type> m;
111 std::vector<value_type> w;
112 std::vector<std::vector<value_type> > q;
116 struct IMTriangle<3,double>
118 typedef double value_type;
119 static const uint16_type nDim = 2;
120 static const uint16_type nOrder = 3;
121 static const uint16_type nPoints = 6;
127 w += value_type( 1 )/value_type( 6 );
128 q[0] += value_type( 0.659027622374092 ), value_type( 0.231933368553031 ), value_type( 0.109039009072877 );
131 std::vector<uint16_type> m;
132 std::vector<value_type> w;
133 std::vector<std::vector<value_type> > q;
137 struct IMTriangle<4,double>
139 typedef double value_type;
140 static const uint16_type nDim = 2;
141 static const uint16_type nOrder = 4;
142 static const uint16_type nPoints = 6;
148 w += 0.109951743655322, 0.223381589678011;
149 q[0] += 0.816847572980459, 0.091576213509771, 0.091576213509771;
150 q[1] += 0.108103018168070, 0.445948490915965, 0.445948490915965;
153 std::vector<uint16_type> m;
154 std::vector<value_type> w;
155 std::vector<std::vector<value_type> > q;
159 struct IMTriangle<5,double>
161 typedef double value_type;
162 static const uint16_type nDim = 2;
163 static const uint16_type nOrder = 5;
164 static const uint16_type nPoints = 7;
170 w += 0.22500000000000, 0.125939180544827, 0.132394152788506;
171 q[0] += value_type( 1.0 )/value_type( 3.0 ), value_type( 1.0 )/value_type( 3.0 ), value_type( 1.0 )/value_type( 3.0 );
172 q[1] += 0.797426985353087, 0.101286507323456, 0.101286507323456;
173 q[2] += 0.059715871789770, 0.470142064105115, 0.470142064105115;
176 std::vector<uint16_type> m;
177 std::vector<value_type> w;
178 std::vector<std::vector<value_type> > q;
182 struct IMTriangle<6,double>
184 typedef double value_type;
185 static const uint16_type nDim = 2;
186 static const uint16_type nOrder = 6;
187 static const uint16_type nPoints = 12;
193 w += 0.050844906370207, 0.116786275726379, 0.082851075618374;
194 q[0] += 0.873821971016996, 0.063089014491502, 0.063089014491502;
195 q[1] += 0.501426509658179, 0.249286745170910, 0.249286745170910;
196 q[2] += 0.636502499121399, 0.310352451033784, 0.053145049844817;
199 std::vector<uint16_type> m;
200 std::vector<value_type> w;
201 std::vector<std::vector<value_type> > q;
210 template<
int Order,
typename T>
216 struct IMTetrahedra<1,T>
218 typedef T value_type;
219 static const uint16_type nDim = 3;
220 static const uint16_type nOrder = 1;
221 static const uint16_type nPoints = 1;
229 q[0] += value_type( 1 )/value_type( 4 );
232 std::vector<uint16_type> m;
233 std::vector<value_type> w;
234 std::vector<std::vector<value_type> > q;
236 template<
typename T>
struct IMTetrahedra<0,T>:
public IMTetrahedra<1,T> {};
238 struct IMTetrahedra<2,T>
240 typedef T value_type;
241 static const uint16_type nDim = 3;
242 static const uint16_type nOrder = 2;
243 static const uint16_type nPoints = 4;
249 w += value_type( 1 )/value_type( 4 );
250 q[0] += 0.5854101966249685, 0.1381966011250105;
253 std::vector<uint16_type> m;
254 std::vector<value_type> w;
255 std::vector<std::vector<value_type> > q;
259 struct IMTetrahedra<4,T>
261 typedef T value_type;
262 static const uint16_type nDim = 3;
263 static const uint16_type nOrder = 4;
264 static const uint16_type nPoints = 16;
270 w += 0.05037379410012282, 0.06654206863329239;
271 q[0] += 0.7716429020672371, 0.7611903264425430e-01;
272 q[1] += 0.1197005277978019, 0.7183164526766925e-01, 0.4042339134672644;
275 std::vector<uint16_type> m;
276 std::vector<value_type> w;
277 std::vector<std::vector<value_type> > q;
279 template<
typename T>
struct IMTetrahedra<3,T>:
public IMTetrahedra<4,T> {};
281 struct IMTetrahedra<6,T>
283 typedef T value_type;
284 static const uint16_type nDim = 3;
285 static const uint16_type nOrder = 6;
286 static const uint16_type nPoints = 29;
292 w += 0.9040129046014750e-01, 0.1911983427899124e-01,
293 0.4361493840666568e-01, 0.2581167596199161e-01;
296 q[1] += 0.8277192480479295, 0.5742691731735683e-01;
297 q[2] += 0.5135188412556341e-01, 0.4860510285706072, 0.2312985436519147;
298 q[3] += 0.2967538129690260, 0.6081079894015281, 0.4756909881472290e-01;
301 std::vector<uint16_type> m;
302 std::vector<value_type> w;
303 std::vector<std::vector<value_type> > q;
305 template<
typename T>
struct IMTetrahedra<5,T>:
public IMTetrahedra<6,T> {};
317 template<
int Dim,
int Order,
typename T>
320 public PointSetQuadrature<Simplex<Dim,1> , Order, T>
322 typedef PointSetQuadrature<Simplex<Dim,1> , Order, T> super;
328 typedef T value_type;
329 typedef ublas::matrix<value_type,ublas::column_major> matrix_type;
330 typedef ublas::vector<value_type> vector_type;
331 typedef typename mpl::if_<mpl::equal_to<mpl::int_<Dim>,mpl::int_<2> >,
332 mpl::identity<detail::IMTriangle<Order,T> >,
333 mpl::identity<detail::IMTetrahedra<Order,T> > >::type::type quad_type;
336 typedef typename mpl::if_<mpl::equal_to<mpl::int_<Dim>,mpl::int_<2> >,
337 mpl::identity<Gauss<Simplex<Dim-1,1>,Order,T> >,
338 mpl::identity<
IMSimplex<Dim-1,Order,T> > >::type::type face_quad_type;
340 typedef IMSimplex<Dim-1,Order,T> face_quad_type;
342 static const uint16_type nDim = Dim;
343 static const uint16_type nOrder = Order;
352 super( quad_type::nPoints ),
356 for (
size_type i=0, wi = 0; i< M_quad.q.size(); i++ )
358 permute( M_quad.m[i], M_quad.q[i], wi );
360 for (
int j=0; j<M_quad.m[i]; j++, wi++ )
362 this->M_w( wi ) = factor()*M_quad.w[i];
370 typedef GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element_type;
371 RealToReference<element_type, Simplex,value_type> real_to_ref( element() );
372 this->setPoints( real_to_ref( this->M_points ) );
373 this->M_w /= real_to_ref.J();
380 boost::shared_ptr<GT_Lagrange<Dim,1,Simplex,T> > gm(
new GT_Lagrange<Dim, 1, Simplex, T> );
381 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
383 this->constructQROnFace( Reference<Simplex<Dim, 1>, Dim, 1>(), gm, face_qr );
393 return ( Dim==2 )?T( 1 )/T( 2 ):T( 1 )/T( 6 );
400 GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element()
const
402 return element( mpl::int_<nDim>() );
405 GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element( mpl::int_<2> )
const
408 typename node<value_type>::type pt( nDim );
409 typedef GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element_type;
410 typedef typename element_type::point_type point_type;
413 pt[0] = value_type( 0 );
414 pt[1] = value_type( 0 );
416 elem.setPoint( 0, p1 );
419 pt[0] = value_type( 1 );
420 pt[1] = value_type( 0 );
422 elem.setPoint( 1, p2 );
425 pt[0] = value_type( 0 );
426 pt[1] = value_type( 1 );
428 elem.setPoint( 2, p3 );
433 GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element( mpl::int_<3> )
const
436 typename node<value_type>::type pt( nDim );
437 typedef GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element_type;
438 typedef typename element_type::point_type point_type;
441 pt[0] = value_type( 0 );
442 pt[1] = value_type( 0 );
443 pt[2] = value_type( 0 );
445 elem.setPoint( 0, p1 );
448 pt[0] = value_type( 1 );
449 pt[1] = value_type( 0 );
450 pt[2] = value_type( 0 );
452 elem.setPoint( 1, p2 );
455 pt[0] = value_type( 0 );
456 pt[1] = value_type( 1 );
457 pt[2] = value_type( 0 );
459 elem.setPoint( 2, p3 );
462 pt[0] = value_type( 0 );
463 pt[1] = value_type( 0 );
464 pt[2] = value_type( 1 );
466 elem.setPoint( 3, p4 );
471 template<
typename QVec>
473 permute(
int m, QVec
const& q,
int wi )
475 permute( m, q, wi, mpl::int_<Dim>() );
477 template<
typename QVec>
479 permute(
int m, QVec
const& q,
int wi, mpl::int_<2> )
485 this->M_points( 0, wi ) = q[ 0 ];
486 this->M_points( 1, wi ) = q[ 1 ];
493 this->M_points( 0, wi ) = q[ 0 ];
494 this->M_points( 1, wi ) = q[ 1 ];
497 this->M_points( 0, wi ) = q[ 1 ];
498 this->M_points( 1, wi ) = q[ 0 ];
501 this->M_points( 0, wi ) = q[ 2 ];
502 this->M_points( 1, wi ) = q[ 1 ];
511 this->M_points( 0, wi ) = q[ 0 ];
512 this->M_points( 1, wi ) = q[ 1 ];
516 this->M_points( 0, wi ) = q[ 0 ];
517 this->M_points( 1, wi ) = q[ 2 ];
521 this->M_points( 0, wi ) = q[ 1 ];
522 this->M_points( 1, wi ) = q[ 0 ];
526 this->M_points( 0, wi ) = q[ 1 ];
527 this->M_points( 1, wi ) = q[ 2 ];
531 this->M_points( 0, wi ) = q[ 2 ];
532 this->M_points( 1, wi ) = q[ 1 ];
536 this->M_points( 0, wi ) = q[ 2 ];
537 this->M_points( 1, wi ) = q[ 0 ];
543 template<
typename QVec>
545 permute(
int m, QVec
const& q,
int wi, mpl::int_<3> )
550 this->M_points( 0, wi ) = q[ 0 ];
551 this->M_points( 1, wi ) = q[ 0 ];
552 this->M_points( 2, wi ) = q[ 0 ];
558 this->M_points( 0, wi ) = q[ 0 ];
559 this->M_points( 1, wi ) = q[ 1 ];
560 this->M_points( 2, wi ) = q[ 1 ];
564 this->M_points( 0, wi ) = q[ 1 ];
565 this->M_points( 1, wi ) = q[ 0 ];
566 this->M_points( 2, wi ) = q[ 1 ];
570 this->M_points( 0, wi ) = q[ 1 ];
571 this->M_points( 1, wi ) = q[ 1 ];
572 this->M_points( 2, wi ) = q[ 0 ];
576 this->M_points( 0, wi ) = q[ 1 ];
577 this->M_points( 1, wi ) = q[ 1 ];
578 this->M_points( 2, wi ) = q[ 1 ];
585 this->M_points( 0, wi ) = q[ 0 ];
586 this->M_points( 1, wi ) = q[ 1 ];
587 this->M_points( 2, wi ) = q[ 2 ];
591 this->M_points( 0, wi ) = q[ 0 ];
592 this->M_points( 1, wi ) = q[ 2 ];
593 this->M_points( 2, wi ) = q[ 1 ];
597 this->M_points( 0, wi ) = q[ 0 ];
598 this->M_points( 1, wi ) = q[ 2 ];
599 this->M_points( 2, wi ) = q[ 2 ];
604 this->M_points( 0, wi ) = q[ 1 ];
605 this->M_points( 1, wi ) = q[ 0 ];
606 this->M_points( 2, wi ) = q[ 2 ];
610 this->M_points( 0, wi ) = q[ 2 ];
611 this->M_points( 1, wi ) = q[ 0 ];
612 this->M_points( 2, wi ) = q[ 1 ];
616 this->M_points( 0, wi ) = q[ 2 ];
617 this->M_points( 1, wi ) = q[ 0 ];
618 this->M_points( 2, wi ) = q[ 2 ];
622 this->M_points( 0, wi ) = q[ 1 ];
623 this->M_points( 1, wi ) = q[ 2 ];
624 this->M_points( 2, wi ) = q[ 0 ];
628 this->M_points( 0, wi ) = q[ 2 ];
629 this->M_points( 1, wi ) = q[ 1 ];
630 this->M_points( 2, wi ) = q[ 0 ];
634 this->M_points( 0, wi ) = q[ 2 ];
635 this->M_points( 1, wi ) = q[ 2 ];
636 this->M_points( 2, wi ) = q[ 0 ];
640 this->M_points( 0, wi ) = q[ 1 ];
641 this->M_points( 1, wi ) = q[ 2 ];
642 this->M_points( 2, wi ) = q[ 2 ];
646 this->M_points( 0, wi ) = q[ 2 ];
647 this->M_points( 1, wi ) = q[ 1 ];
648 this->M_points( 2, wi ) = q[ 2 ];
652 this->M_points( 0, wi ) = q[ 2 ];
653 this->M_points( 1, wi ) = q[ 2 ];
654 this->M_points( 2, wi ) = q[ 1 ];
660 FEELPP_ASSERT( 0 )( m ).error(
"invalid multiplicity in IMSimplex<3>" );
667 return test( mpl::int_<nDim>() );
669 bool test( mpl::int_<2> )
672 ublas::vector<value_type> x( ublas::row( this->M_points, 0 ) );
673 ublas::vector<value_type> y( ublas::row( this->M_points, 1 ) );
674 ublas::vector<value_type> w( this->M_w );
676 for (
int a=0; a<=nOrder; a++ )
678 for (
int b=0; b<nOrder-a; b++ )
680 int cMax = nOrder - a - b;
682 for (
int c=0; c<=cMax; c++ )
686 for (
int q=0; q<w.size(); q++ )
688 sum += 0.5*w[q] * pow( x[q], (
double ) a ) * pow( y[q], (
double ) b )
689 * pow( 1.0 - x[q] - y[q], (
double ) c );
692 double err = fabs( sum - exact( a,b,c ) );
693 bool localPass = err < 1.0e-14;
694 pass = pass && localPass;
698 fprintf( stderr,
"order=%d m (%d, %d, %d) q=%22.15g exact=%22.15g\n", nOrder, a, b, c, sum, exact( a, b, c ) );
699 std::cerr <<
"error = " << err << std::endl;
704 fprintf( stderr,
"order=%d m (%d, %d, %d) q=%22.15g exact=%22.15g\n", nOrder, a, b, c, sum, exact( a, b, c ) );
712 bool test( mpl::int_<3> )
714 ublas::vector<value_type> x( ublas::row( this->M_points, 0 ) );
715 ublas::vector<value_type> y( ublas::row( this->M_points, 1 ) );
716 ublas::vector<value_type> z( ublas::row( this->M_points, 2 ) );
717 ublas::vector<value_type> w( this->M_w );
722 for (
int a=0; a<=p; a++ )
724 for (
int b=0; b<p-a; b++ )
726 int cMax = p - a - b;
728 for (
int c=0; c<=cMax; c++ )
730 int dMax = p - a - b - c;
732 for (
int d=0; d<=dMax; d++ )
736 for (
int q=0; q<w.size(); q++ )
738 sum += ( 1.0/6.0 )*w[q] * pow( x[q], (
double ) a ) * pow( y[q], (
double ) b )
739 * pow( z[q], (
double ) c )
740 * pow( 1.0 - x[q] - y[q] - z[q], (
double ) d );
743 double err = fabs( sum - exact( a,b,c,d ) );
744 bool localPass = err < 1.0e-14;
745 pass = pass && localPass;
749 fprintf( stderr,
"order=%d m (%d, %d, %d %d) q=%22.15g exact=%22.15g\n", p, a, b, c, d, sum, exact( a, b, c, d ) );
750 std::cerr <<
"error = " << err << std::endl;
755 fprintf( stderr,
"order=%d m (%d, %d, %d %d) q=%22.15g exact=%22.15g\n", p, a, b, c, d, sum, exact( a, b, c, d ) );
766 double exact(
int a,
int b,
int c )
const
768 return fact( a )*fact( b )*fact( c )/fact( a+b+c+2 );
770 double exact(
int a,
int b,
int c,
int d )
const
772 return fact( a )*fact( b )*fact( c )*fact( d )/fact( a+b+c+d+3 );
775 double fact(
int x )
const
777 if ( x==0 )
return 1.0;
779 return x*fact( x-1 );
794 template<
typename T >
795 class TriangleQuadRule
799 typedef T value_type;
800 typedef typename node<value_type>::type node_type;
801 typedef ublas::matrix<value_type, ublas::column_major> nodes_type;
802 typedef ublas::vector<value_type> weights_type;
804 nodes_type
const&
points()
const
809 ublas::matrix_column<nodes_type const> point( uint32_type __i )
const
811 return ublas::column( _pts, __i );
814 weights_type
const& weights()
const
818 value_type
const& weight(
int q )
const
823 value_type
integrate( boost::function<value_type( node_type
const& )>
const& f )
const
825 value_type res = 0.0;
827 for ( uint16_type k = 0; k < 7 ; ++k )
829 res += this->weights()[k]*f( this->point( k ) );
836 : _pts( 2, 7 ), M_w( 7 )
838 value_type t7pt_x[3] = {value_type( 1.0 )/3.0, 0.10128650732345633, 0.47014206410511508};
839 value_type t7pt_w[3] = {0.1125, 0.062969590272413576, 0.066197076394253090};
841 for (
int i = 0; i < 3 ; ++i )
843 t7pt_x[i]= 2.0*t7pt_x[i]-1.0;
844 t7pt_w[i]= 4.0*t7pt_w[i];
847 M_w( 0 ) = t7pt_w[0];
848 M_w( 1 ) = t7pt_w[1];
849 M_w( 2 ) = t7pt_w[1];
850 M_w( 3 ) = t7pt_w[1];
851 M_w( 4 ) = t7pt_w[2];
852 M_w( 5 ) = t7pt_w[2];
853 M_w( 6 ) = t7pt_w[2];
855 _pts( 0, 0 ) = t7pt_x[0];
856 _pts( 1, 0 ) = t7pt_x[0];
858 _pts( 0, 1 ) = t7pt_x[1];
859 _pts( 1, 1 ) = t7pt_x[1];
861 _pts( 0, 2 ) = t7pt_x[1];
862 _pts( 1, 2 ) = 1.0-2.0*t7pt_x[1];
864 _pts( 0, 3 ) = 1.0-2.0*t7pt_x[1];
865 _pts( 1, 3 ) = t7pt_x[1];
867 _pts( 0, 4 ) = t7pt_x[2];
868 _pts( 1, 4 ) = t7pt_x[2];
870 _pts( 0, 5 ) = t7pt_x[2];
871 _pts( 1, 5 ) = 1.0-2.0*t7pt_x[2];
873 _pts( 0, 6 ) = 1.0-2.0*t7pt_x[2];
874 _pts( 1, 6 ) = t7pt_x[2];
877 ~TriangleQuadRule() {}