29 #ifndef __GaussLobatto_H
30 #define __GaussLobatto_H 1
37 #include <boost/numeric/ublas/matrix.hpp>
38 #include <boost/numeric/ublas/io.hpp>
39 #include <boost/numeric/ublas/vector.hpp>
40 #include <boost/numeric/ublas/vector_proxy.hpp>
43 #include <feel/feelcore/traits.hpp>
46 #include <feel/feelpoly/quadpoint.hpp>
55 namespace ublas = boost::numeric::ublas;
59 template<
int Dim,
int Order,
int RealDim,
template<u
int16_type,u
int16_type,u
int16_type>
class Entity,
typename T>
struct GT_Lagrange;
75 template<
class Convex, u
int16_type Integration_Degree,
typename T>
79 template< u
int16_type Integration_Degree,
typename T>
86 typedef typename super::return_type return_type;
87 typedef typename super::node_type node_type;
88 typedef typename super::nodes_type nodes_type;
89 typedef typename super::weights_type weights_type;
91 static const uint16_type Degree = ( Integration_Degree+3 )/2+1;
92 static const uint32_type Npoints = Degree;
98 ublas::vector<T> px( Npoints );
100 details::gausslobattojacobi<Npoints, T, ublas::vector<T>, ublas::vector<T> >( this->M_w, px );
101 ublas::row( this->M_points, 0 ) = px;
106 FEELPP_DEFINE_VISITABLE();
112 template< u
int16_type Integration_Degree,
typename T>
113 class GaussLobatto<Simplex<2,1> , Integration_Degree ,T > :
public PointSetQuadrature<Simplex<2,1> , Integration_Degree, T>
116 typedef T value_type;
118 typedef PointSetQuadrature<Simplex<2,1> , Integration_Degree, T> super;
119 typedef typename super::return_type return_type;
120 typedef typename super::node_type node_type;
121 typedef typename super::nodes_type nodes_type;
122 typedef typename super::weights_type weights_type;
124 typedef GaussLobatto<Simplex<1,1>,Integration_Degree, T> face_quad_type;
127 static const uint16_type DegreeX = ( Integration_Degree+3 )/2+1;
128 static const uint16_type DegreeY = ( Integration_Degree+2 )/2+1;
129 static const uint32_type Npoints = DegreeX*DegreeY;
136 weights_type wx( DegreeX );
137 weights_type px( DegreeX );
138 details::gausslobattojacobi<DegreeX,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
140 weights_type wy( DegreeY );
141 weights_type py( DegreeY );
142 details::left_gaussradaujacobi<DegreeY,T, ublas::vector<T>, ublas::vector<T> >( wy, py, 1.0, 0.0 );
146 details::xi<TRIANGLE, value_type> to_xi;
148 for (
int i = 0, k = 0; i < DegreeX; ++i )
150 for (
int j = 0; j < DegreeY; ++j, ++k )
153 this->M_w( k ) = 0.5 * wx( i ) * wy( j );
158 ublas::column( this->M_points, k ) = to_xi( eta );
162 boost::shared_ptr<GT_Lagrange<2,1, 2 ,Simplex, T> > gm(
new GT_Lagrange<2, 1, 2, Simplex,T> );
163 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
165 this->constructQROnFace( Reference<Simplex<2, 1, 2>,2,1>(), gm, face_qr );
171 FEELPP_DEFINE_VISITABLE();
176 template< u
int16_type Integration_Degree,
typename T>
177 class GaussLobatto<Simplex<3,1> , Integration_Degree ,T > :
public PointSetQuadrature<Simplex<3,1> , Integration_Degree, T>
180 typedef T value_type;
182 typedef PointSetQuadrature<Simplex<3,1> , Integration_Degree, T> super;
183 typedef typename super::return_type return_type;
184 typedef typename super::node_type node_type;
185 typedef typename super::nodes_type nodes_type;
186 typedef typename super::weights_type weights_type;
188 typedef GaussLobatto<Simplex<2,1>,Integration_Degree, T> face_quad_type;
191 static const uint16_type DegreeX = ( Integration_Degree+3 )/2+1;
192 static const uint16_type DegreeY = ( Integration_Degree+2 )/2+1;
193 static const uint16_type DegreeZ = ( Integration_Degree+2 )/2+1;
194 static const uint32_type Npoints = DegreeX*DegreeY*DegreeZ;
201 weights_type wx( DegreeX );
202 weights_type px( DegreeX );
203 details::gausslobattojacobi<DegreeX,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
205 weights_type wy( DegreeY );
206 weights_type py( DegreeY );
207 details::left_gaussradaujacobi<DegreeY,T, ublas::vector<T>, ublas::vector<T> >( wy, py, 1.0, 0.0 );
209 weights_type wz( DegreeZ );
210 weights_type pz( DegreeZ );
211 details::left_gaussradaujacobi<DegreeZ,T, ublas::vector<T>, ublas::vector<T> >( wz, pz, 2.0, 0.0 );
215 details::xi<TETRAHEDRON, value_type> to_xi;
217 for (
int i = 0, k = 0; i < DegreeX; ++i )
219 for (
int j = 0; j < DegreeY; ++j )
221 for (
int l = 0; l < DegreeZ; ++l, ++k )
224 this->M_w( k ) = 0.125 * wx( i ) * wy( j ) * wz( l );
230 ublas::column( this->M_points, k ) = to_xi( eta );
235 boost::shared_ptr<GT_Lagrange<3, 1, 3, Simplex, T> > gm(
new GT_Lagrange<3, 1, 3, Simplex, T> );
236 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
238 this->constructQROnFace( Reference<Simplex<3, 1, 3>,3,1>(), gm, face_qr );
243 FEELPP_DEFINE_VISITABLE();
253 template< u
int16_type Integration_Degree,
typename T>
254 class GaussLobatto<Hypercube<2,1>, Integration_Degree ,T >
256 public PointSetQuadrature<Hypercube<2,1>, Integration_Degree, T>
259 typedef T value_type;
261 typedef PointSetQuadrature<Hypercube<2,1>, Integration_Degree, T> super;
262 typedef typename super::return_type return_type;
263 typedef typename super::node_type node_type;
264 typedef typename super::nodes_type nodes_type;
265 typedef typename super::weights_type weights_type;
266 typedef GaussLobatto<Hypercube<1,1>,Integration_Degree, T> face_quad_type;
267 static const uint16_type Degree = ( Integration_Degree+3 )/2+1;
268 static const uint32_type Npoints = Degree*Degree;
275 weights_type wx( Degree );
276 weights_type px( Degree );
277 details::gausslobattojacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
279 for (
int i = 0, k = 0; i < Degree; ++i )
281 for (
int j = 0; j < Degree; ++j, ++k )
284 this->M_w( k ) = wx( i ) * wx( j );
285 this->M_points( 0, k ) = px( i );
286 this->M_points( 1, k ) = px( j );
290 boost::shared_ptr<GT_Lagrange<2, 1, 2, Hypercube, T> > gm(
new GT_Lagrange<2, 1, 2, Hypercube, T> );
291 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
293 this->constructQROnFace( Reference<Hypercube<2, 1, 2>,2,1>(), gm, face_qr );
298 FEELPP_DEFINE_VISITABLE();
303 template< u
int16_type Integration_Degree,
typename T>
304 class GaussLobatto<Hypercube<3,1>, Integration_Degree ,T >
306 public PointSetQuadrature<Hypercube<3,1>, Integration_Degree, T>
309 typedef T value_type;
311 typedef PointSetQuadrature<Hypercube<3,1>, Integration_Degree, T> super;
312 typedef typename super::return_type return_type;
313 typedef typename super::node_type node_type;
314 typedef typename super::nodes_type nodes_type;
315 typedef typename super::weights_type weights_type;
316 typedef GaussLobatto<Hypercube<2,1>,Integration_Degree, T> face_quad_type;
317 static const uint16_type Degree = ( Integration_Degree+3 )/2+1;
318 static const uint32_type Npoints = Degree*Degree*Degree;
325 weights_type wx( Degree );
326 weights_type px( Degree );
327 details::gausslobattojacobi<Degree,T, ublas::vector<T>, ublas::vector<T> >( wx, px, 0.0, 0.0 );
329 for (
int i = 0, k = 0; i < Degree; ++i )
331 for (
int j = 0; j < Degree; ++j )
333 for (
int l = 0; l < Degree ; ++l, ++k )
336 this->M_w( k ) = wx( i ) * wx( j ) * wx( l );
337 this->M_points( 0, k ) = px( i );
338 this->M_points( 1, k ) = px( j );
339 this->M_points( 2, k ) = px( l );
344 boost::shared_ptr<GT_Lagrange<3, 1, 3, Hypercube, T> > gm(
new GT_Lagrange<3, 1, 3, Hypercube, T> );
345 boost::shared_ptr<face_quad_type> face_qr(
new face_quad_type );
347 this->constructQROnFace( Reference<Hypercube<3, 1, 3>,3,1>(), gm, face_qr );
351 FEELPP_DEFINE_VISITABLE();
364 template<
class Convex,
366 typename T =
double >
367 class PointSetGaussLobatto :
public PointSetInterpolation<Convex::nDim, Order, T, Hypercube>
371 typedef PointSetInterpolation< Convex::nDim, Order, T, Hypercube> super;
373 typedef typename super::return_type return_type;
375 typedef T value_type;
377 static const uint32_type Dim = Convex::nDim;
378 static const uint32_type nPoints = Order+1;
379 static const uint32_type topological_dimension = Convex::topological_dimension;
380 static const uint32_type nRealDim = Convex::nRealDim;
382 static const size_type Shape = Convex::Shape;
384 static const bool is_simplex = Convex::is_simplex;
385 static const bool is_hypercube = Convex::is_hypercube;
387 typedef Reference<Convex, Dim, Convex::nOrder, Convex::nDim, value_type> reference_convex_type;
389 typedef ublas::vector<value_type> vector_type;
391 typedef typename super::node_type node_type;
392 typedef typename super::nodes_type nodes_type;
394 typedef typename reference_convex_type::points_type points_type;
396 typedef typename Convex::edge_to_point_t edge_to_point_t;
397 typedef typename Convex::face_to_point_t face_to_point_t;
398 typedef typename Convex::face_to_edge_t face_to_edge_t;
400 typedef Hypercube<Dim, Order, Dim> conv_order_type;
401 static const uint32_type numPoints = conv_order_type::numPoints;
403 reference_convex_type RefConv;
405 typedef typename super::range_type range_type;
406 typedef typename super::index_map_type index_map_type;
408 PointSetGaussLobatto(
int interior = 0 )
410 FEELPP_ASSERT( is_hypercube || ( Dim == 1 ) ).error(
"gauss lobatto points are just defined in simplex products" );
412 nodes_type pts( Dim, numPoints );
414 calculate_gl_points( mpl::bool_< ( Order > 1 )>() );
416 if ( interior == 0 && Order > 0 )
418 for ( uint16_type d = 0, p = 0; d < topological_dimension+1; ++d )
420 for (
int e = RefConv.entityRange( d ).begin();
421 e < RefConv.entityRange( d ).end();
424 nodes_type Gt ( makePoints( d, e ) );
428 ublas::subrange( pts, 0, Dim, p, p+Gt.size2() ) = Gt;
430 for (
size_type j = 0; j < Gt.size2(); ++j )
432 this->addToEid( d, p+j );
433 this->addToPtE( p+j, std::make_pair( d, e ) );
444 else if ( interior == 1 && Order > 0 )
447 else if ( Order == 0 )
448 this->
setPoints( glas::average( RefConv.vertices() ) );
450 this->setName(
"gausslobatto", Order );
455 ~PointSetGaussLobatto() {}
461 void calculate_gl_points ( mpl::bool_<true> )
463 vector_type wx( Order+1 );
464 vector_type px( Order+1 );
466 details::gausslobattojacobi<Order+1, T, vector_type, vector_type >( wx, px );
468 ublas::vector_range<vector_type> inner_pts ( px, ublas::range( 1, px.size()-1 ) );
473 void calculate_gl_points ( mpl::bool_<false> )
476 points_type makePoints( uint16_type topo_dim, uint16_type __id )
481 points_type G( RefConv.vertices().size1(), 1 );
482 ublas::column( G, 0 ) = ublas::column( RefConv.vertices(), __id );
487 else if ( topo_dim == topological_dimension )
490 return makeLattice<Shape>();
492 throw std::logic_error(
"cannot make those points" );
493 return points_type();
504 G = makeLattice<SHAPE_LINE>();
505 Gret.resize( nRealDim, G.size2() );
507 pt_to_entity_hexahedron<Shape, 1> p_to_e( __id );
509 for (
size_type i = 0; i < G.size2(); ++i )
510 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
515 else if ( topo_dim == 2 )
517 G = makeLattice<SHAPE_QUAD>();
518 Gret.resize( nRealDim, G.size2() );
519 pt_to_entity_hexahedron<Shape, 2> p_to_e( __id );
521 for (
size_type i = 0; i < G.size2(); ++i )
522 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
528 return points_type();
531 template<
size_type shape>
532 points_type makeLattice()
538 if ( shape == SHAPE_LINE )
539 G = make_line_points();
541 else if ( shape == SHAPE_QUAD )
542 return make_quad_points();
544 else if ( shape == SHAPE_HEXA )
545 return make_hexa_points();
548 else if ( Order == 0 )
549 G = glas::average( RefConv.vertices() );
556 return std::max( 0,
int( Order ) - 1 );
559 int n_quad_points()
const
561 return std::max( 0, (
int( Order ) - 1 )*(
int( Order ) - 1 ) );
565 int n_hexa_points()
const
567 return std::max( 0, (
int( Order )-1 )*(
int( Order )-1 )*(
int( Order )-1 ) );
573 points_type G ( Dim, n_line_points() );
577 vector_type ones ( ublas::scalar_vector<value_type>( G.size2(), value_type( 1 ) ) );
579 ublas::row( G, 0 ) = gl_pts;
581 for ( uint16_type i=1; i<Dim; i++ )
582 ublas::row( G, i ) = - ones;
591 points_type G( Dim, n_quad_points() );
595 uint16_type numInterior = Order - 1;
597 for ( uint16_type i = 0, k = 0; i < numInterior; ++i )
599 for ( uint16_type j = 0; j < numInterior; ++j, ++k )
601 G( 0, k ) = gl_pts( i );
602 G( 1, k ) = gl_pts( j );
606 vector_type ones ( ublas::scalar_vector<value_type>( G.size2(), value_type( 1 ) ) );
609 ublas::row( G, 2 ) = - ones;
618 points_type G( Dim, n_hexa_points() );
622 uint16_type numInterior = Order - 1;
624 for ( uint16_type i = 0, k = 0; i < numInterior; ++i )
626 for ( uint16_type j = 0; j < numInterior; ++j )
628 for ( uint16_type l = 0; l < numInterior; ++l, ++k )
630 G( 0, k ) = gl_pts( i );
631 G( 1, k ) = gl_pts( j );
632 G( 2, k ) = gl_pts( l );
641 template<
size_type shape>
644 pt_to_edge( std::vector<uint16_type> vert_ids )
647 a( Entity<SHAPE_LINE, value_type>().vertex( 0 ) ),
648 b( Entity<SHAPE_LINE, value_type>().vertex( 1 ) ),
649 u( Entity<shape, value_type>().vertex( vert_ids[ 0 ] ) ),
650 v( Entity<shape, value_type>().vertex( vert_ids[ 1 ] ) ),
653 h = 1.0/( b[0]-a[0] );
656 operator()( node_type
const& x )
const
658 return u + h * ( x[ 0 ] - a[ 0 ] ) * diff;
662 node_type u, v, diff;
669 template<
size_type shape>
670 struct pt_to_face_hexahedron
672 pt_to_face_hexahedron( std::vector<uint16_type> vert_ids )
674 u( Entity<shape, value_type>().vertex( vert_ids[ 0 ] ) ),
675 v( Entity<shape, value_type>().vertex( vert_ids[ 1 ] ) ),
676 w( Entity<shape, value_type>().vertex( vert_ids[ 3 ] ) ),
683 operator()( node_type
const& x )
const
685 return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
688 ublas::vector<node_type> diff;
691 template<
size_type shape>
695 pt_to_element( std::vector<uint16_type>
const& ) {}
696 node_type operator()( node_type
const& x )
const
702 template<
size_type shape,u
int16_type topo_dim>
703 struct pt_to_entity_hexahedron
705 typedef typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_LINE> >,
706 mpl::identity<mpl::vector<boost::none_t,pt_to_edge<shape>,pt_to_edge<shape> > >,
707 typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_QUAD> >,
708 mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_element<shape> > >,
709 mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_face_hexahedron<shape>, pt_to_element<shape> > >
712 typedef typename mpl::at<_type, mpl::int_<topo_dim> >::type mapping_type;
713 typedef mpl::vector<boost::none_t, edge_to_point_t, face_to_point_t> list_v;
715 pt_to_entity_hexahedron( uint16_type entity_id )
717 mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
720 node_type operator()( node_type
const& x )
const
724 mapping_type mapping;
simplex of dimension Dim
Definition: simplex.hpp:73
Gauss quadrature points.
Definition: gausslobatto.hpp:76
Quadrature point set base class.
Definition: gausslobatto.hpp:58
void setPoints(nodes_type const &pts)
Definition: pointset.hpp:267
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319