32 #include <feel/feelpoly/quadpoint.hpp>
50 template<
class Convex, u
int16_type Integration_Degree,
typename T>
53 template< u
int16_type Integration_Degree,
typename T>
60 typedef typename super::return_type return_type;
61 typedef typename super::node_type node_type;
62 typedef typename super::nodes_type nodes_type;
63 typedef typename super::weights_type weights_type;
66 static const uint32_type Npoints = 1;
77 FEELPP_DEFINE_VISITABLE();
81 template< u
int16_type Integration_Degree,
typename T>
82 class Gauss<Simplex<1,1> , Integration_Degree ,T > :
public PointSetQuadrature<Simplex<1,1> , Integration_Degree, T>
87 typedef PointSetQuadrature<Simplex<1,1> , Integration_Degree, T> super;
88 typedef typename super::return_type return_type;
89 typedef typename super::node_type node_type;
90 typedef typename super::nodes_type nodes_type;
91 typedef typename super::weights_type weights_type;
93 static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
94 static const uint32_type Npoints = Degree;
96 typedef Gauss<Simplex<0,1>,Integration_Degree, T> face_quad_type;
102 ublas::vector<T> px( Npoints );
104 details::gaussjacobi<Npoints, T, ublas::vector<T>, ublas::vector<T> >( this->M_w, px );
105 ublas::row( this->M_points, 0 ) = px;
107 boost::shared_ptr<GT_Lagrange<1,1,1,Simplex,T> > gm(
new GT_Lagrange<1, 1, 1, Simplex, T> );
108 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
110 this->constructQROnFace( Reference<Simplex<1, 1>, 1, 1>(), gm, face_qr );
117 FEELPP_DEFINE_VISITABLE();
122 template< u
int16_type Integration_Degree,
typename T>
123 class Gauss<Simplex<2,1> , Integration_Degree ,T > :
public PointSetQuadrature<Simplex<2,1> , Integration_Degree, T>
126 typedef T value_type;
128 typedef PointSetQuadrature<Simplex<2,1> , Integration_Degree, T> super;
129 typedef typename super::return_type return_type;
131 typedef typename super::node_type node_type;
132 typedef typename super::nodes_type nodes_type;
133 typedef typename super::weights_type weights_type;
135 typedef Gauss<Simplex<1,1>,Integration_Degree, T> face_quad_type;
138 static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
139 static const uint32_type Npoints = Degree*Degree;
146 weights_type wx( Degree );
147 weights_type px( Degree );
148 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
150 weights_type wy( Degree );
151 weights_type py( Degree );
152 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wy, py, 1.0, 0.0 );
157 std::cout<<
"[Debug quadpoint] Npoints = " << Npoints << std::endl ;
158 std::cout<<
"[Debug quadpoint] _pts.size2() = " << this->M_points.size2() << std::endl ;
159 std::cout<<
"[Debug quadpoint] Degree = " << Degree << std::endl ;
163 details::xi<TRIANGLE, value_type> to_xi;
165 for (
int i = 0, k = 0; i < Degree; ++i )
167 for (
int j = 0; j < Degree; ++j, ++k )
172 std::cout<<
"[Debug quadpoint] i = " << i <<
" ; j = " << j <<
" ; k = " << k << std::endl ;
177 this->M_w( k ) = 0.5 * wx( i ) * wy( j );
183 ublas::column( this->M_points, k ) = to_xi( eta );
187 boost::shared_ptr<GT_Lagrange<2,1,2,Simplex,T> > gm(
new GT_Lagrange<2, 1, 2, Simplex, T> );
188 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
190 this->constructQROnFace( Reference<Simplex<2, 1>, 2, 1>(), gm, face_qr );
196 FEELPP_DEFINE_VISITABLE();
202 template< u
int16_type Integration_Degree,
typename T>
203 class Gauss<Simplex<3,1> , Integration_Degree ,T > :
public PointSetQuadrature<Simplex<3,1> , Integration_Degree, T>
206 typedef T value_type;
208 typedef PointSetQuadrature<Simplex<3,1> , Integration_Degree, T> super;
210 typedef typename super::return_type return_type;
211 typedef typename super::node_type node_type;
212 typedef typename super::nodes_type nodes_type;
213 typedef typename super::weights_type weights_type;
215 typedef Gauss<Simplex<2,1>,Integration_Degree, T> face_quad_type;
218 static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
219 static const uint32_type Npoints = Degree*Degree*Degree;
226 weights_type wx( Degree );
227 weights_type px( Degree );
228 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
230 weights_type wy( Degree );
231 weights_type py( Degree );
232 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wy, py, 1.0, 0.0 );
234 weights_type wz( Degree );
235 weights_type pz( Degree );
236 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wz, pz, 2.0, 0.0 );
240 details::xi<TETRAHEDRON, value_type> to_xi;
242 for (
int i = 0, k = 0; i < Degree; ++i )
244 for (
int j = 0; j < Degree; ++j )
246 for (
int l = 0; l < Degree; ++l, ++k )
249 this->M_w( k ) = 0.125 * wx( i ) * wy( j ) * wz( l );
255 ublas::column( this->M_points, k ) = to_xi( eta );
260 boost::shared_ptr<GT_Lagrange<3, 1, 3, Simplex, T> > gm(
new GT_Lagrange<3, 1, 3, Simplex, T> );
261 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
263 this->constructQROnFace( Reference<Simplex<3, 1>,3,1>(), gm, face_qr );
268 FEELPP_DEFINE_VISITABLE();
273 template< u
int16_type Integration_Degree,
typename T>
274 class Gauss<Hypercube<1,1>, Integration_Degree ,T >
276 public PointSetQuadrature<Hypercube<1,1>, Integration_Degree, T>
279 typedef T value_type;
281 typedef PointSetQuadrature<Hypercube<1,1>, Integration_Degree, T> super;
282 typedef typename super::return_type return_type;
283 typedef typename super::node_type node_type;
284 typedef typename super::nodes_type nodes_type;
285 typedef typename super::weights_type weights_type;
287 static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
288 static const uint32_type Npoints = Degree;
290 typedef Gauss<Simplex<0,1>,Integration_Degree, T> face_quad_type;
297 weights_type wx( Degree );
298 weights_type px( Degree );
299 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
301 VLOG(1) <<
"[gauss<SP<2,1>] jacobi p = " << px <<
"\n";
302 VLOG(1) <<
"[gauss<SP<2,1>] jacobi w = " << wx <<
"\n";
305 for (
int i = 0; i < Degree; ++i )
308 this->M_w( i ) = wx( i );
309 this->M_points( 0, i ) = px( i );
314 VLOG(1) <<
"[gauss<SP<2,1>] p = " << this->M_points <<
"\n";
315 VLOG(1) <<
"[gauss<SP<2,1>] w = " << this->M_w <<
"\n";
319 boost::shared_ptr<GT_Lagrange<1,1,1, Hypercube,T> > gm(
new GT_Lagrange<1, 1, 1, Hypercube, T> );
320 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
322 this->constructQROnFace( Reference<Hypercube<1, 1>, 1, 1>(), gm, face_qr );
328 FEELPP_DEFINE_VISITABLE();
332 template< u
int16_type Integration_Degree,
typename T>
333 class Gauss<Hypercube<2,1>, Integration_Degree ,T >
335 public PointSetQuadrature<Hypercube<2,1>, Integration_Degree, T>
338 typedef T value_type;
340 typedef PointSetQuadrature<Hypercube<2,1>, Integration_Degree, T> super;
341 typedef typename super::return_type return_type;
342 typedef typename super::node_type node_type;
343 typedef typename super::nodes_type nodes_type;
344 typedef typename super::weights_type weights_type;
346 typedef Gauss<Hypercube<1,1>,Integration_Degree, T> face_quad_type;
348 static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
349 static const uint32_type Npoints = Degree*Degree;
356 weights_type wx( Degree );
357 weights_type px( Degree );
358 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
360 VLOG(1) <<
"[gauss<SP<2,1>] jacobi p = " << px <<
"\n";
361 VLOG(1) <<
"[gauss<SP<2,1>] jacobi w = " << wx <<
"\n";
364 for (
int i = 0, k = 0; i < Degree; ++i )
366 for (
int j = 0; j < Degree; ++j, ++k )
369 this->M_w( k ) = wx( i ) * wx( j );
370 this->M_points( 0, k ) = px( i );
371 this->M_points( 1, k ) = px( j );
376 VLOG(1) <<
"[gauss<SP<2,1>] p = " << this->M_points <<
"\n";
377 VLOG(1) <<
"[gauss<SP<2,1>] w = " << this->M_w <<
"\n";
379 boost::shared_ptr<GT_Lagrange<2, 1, 2, Hypercube, T> > gm(
new GT_Lagrange<2, 1, 2, Hypercube, T> );
380 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
382 this->constructQROnFace( Reference<Hypercube<2, 1>,2,1>(), gm, face_qr );
388 FEELPP_DEFINE_VISITABLE();
393 template< u
int16_type Integration_Degree,
typename T>
394 class Gauss<Hypercube<3,1>, Integration_Degree ,T >
396 public PointSetQuadrature<Hypercube<3,1>, Integration_Degree, T>
399 typedef T value_type;
401 typedef PointSetQuadrature<Hypercube<3,1>, Integration_Degree, T> super;
402 typedef typename super::return_type return_type;
403 typedef typename super::node_type node_type;
404 typedef typename super::nodes_type nodes_type;
405 typedef typename super::weights_type weights_type;
406 typedef Gauss<Hypercube<2,1>,Integration_Degree, T> face_quad_type;
407 static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
408 static const uint32_type Npoints = Degree*Degree*Degree;
415 weights_type wx( Degree );
416 weights_type px( Degree );
417 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
419 for (
int i = 0, k = 0; i < Degree; ++i )
421 for (
int j = 0; j < Degree; ++j )
423 for (
int l = 0; l < Degree ; ++l, ++k )
426 this->M_w( k ) = wx( i ) * wx( j ) * wx( l );
427 this->M_points( 0, k ) = px( i );
428 this->M_points( 1, k ) = px( j );
429 this->M_points( 2, k ) = px( l );
434 boost::shared_ptr<GT_Lagrange<3, 1, 3, Hypercube, T> > gm(
new GT_Lagrange<3, 1, 3, Hypercube, T> );
435 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
437 this->constructQROnFace( Reference<Hypercube<3, 1>,3,1>(), gm, face_qr );
443 FEELPP_DEFINE_VISITABLE();
447 template< u
int16_type Integration_Degree,
typename T>
448 class Gauss<Hypercube<4,1>, Integration_Degree ,T >
450 public PointSetQuadrature<Hypercube<4,1>, Integration_Degree, T>
453 typedef T value_type;
455 typedef PointSetQuadrature<Hypercube<4,1>, Integration_Degree, T> super;
456 typedef typename super::return_type return_type;
457 typedef typename super::node_type node_type;
458 typedef typename super::nodes_type nodes_type;
459 typedef typename super::weights_type weights_type;
461 static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
462 static const uint32_type Npoints = Degree*Degree*Degree*Degree;
469 weights_type wx( Degree );
470 weights_type px( Degree );
471 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
473 for (
int i = 0, k = 0; i < Degree; ++i )
475 for (
int j = 0; j < Degree; ++j )
477 for (
int l = 0; l < Degree ; ++l )
479 for (
int r = 0; r < Degree ; ++r, ++k )
482 this->M_w( k ) = wx( i ) * wx( j ) * wx( l ) * wx( r );
483 this->M_points( 0, k ) = px( i );
484 this->M_points( 1, k ) = px( j );
485 this->M_points( 2, k ) = px( l );
486 this->M_points( 3, k ) = px( r );
495 FEELPP_DEFINE_VISITABLE();
499 template< u
int16_type Integration_Degree,
typename T>
500 class Gauss<Hypercube<5,1>, Integration_Degree ,T >
502 public PointSetQuadrature<Hypercube<5,1>, Integration_Degree, T>
505 typedef T value_type;
507 typedef PointSetQuadrature<Hypercube<5,1>, Integration_Degree, T> super;
508 typedef typename super::return_type return_type;
509 typedef typename super::node_type node_type;
510 typedef typename super::nodes_type nodes_type;
511 typedef typename super::weights_type weights_type;
513 static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
514 static const uint32_type Npoints = Degree*Degree*Degree*Degree*Degree;
521 weights_type wx( Degree );
522 weights_type px( Degree );
523 details::gaussjacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
525 for (
int i = 0, k = 0; i < Degree; ++i )
527 for (
int j = 0; j < Degree; ++j )
529 for (
int l = 0; l < Degree ; ++l )
531 for (
int r = 0; r < Degree ; ++r )
533 for (
int s = 0; s < Degree ; ++s, ++k )
536 this->M_w( k ) = wx( i ) * wx( j ) * wx( l ) * wx( r ) * wx ( s );
537 this->M_points( 0, k ) = px( i );
538 this->M_points( 1, k ) = px( j );
539 this->M_points( 2, k ) = px( l );
540 this->M_points( 3, k ) = px( r );
541 this->M_points( 4, k ) = px( s );
551 FEELPP_DEFINE_VISITABLE();
simplex of dimension Dim
Definition: simplex.hpp:73
const uint16_type invalid_uint16_type_value
Definition: feelcore/feel.hpp:338
Gauss quadrature points.
Definition: gauss.hpp:51
Quadrature point set base class.
Definition: gausslobatto.hpp:58