32 #define __Quadpoint_H 1
36 #include <feel/feelpoly/geomap.hpp>
40 namespace ublas = boost::numeric::ublas;
42 enum IntegrationFaceEnum
64 template<
class Convex, u
int16_type Integration_Degree,
typename T>
65 class PointSetQuadrature :
public PointSet<Convex,T>
68 static const bool is_face_im =
false;
70 typedef PointSet<Convex,value_type> super;
71 typedef typename super::return_type return_type;
72 typedef typename super::node_type node_type;
73 typedef typename super::nodes_type nodes_type;
74 typedef Eigen::Matrix<value_type,Eigen::Dynamic,1> vector_type;
75 typedef ublas::vector<value_type> weights_type;
77 typedef PointSetQuadrature<Convex, Integration_Degree, T> self_type;
79 typedef self_type parent_quadrature_type;
80 static const uint16_type I_deg = Integration_Degree;
82 PointSetQuadrature(): super(), M_w(), M_prod(), M_exprq() {}
84 PointSetQuadrature(
const PointSetQuadrature& Qp )
88 M_n_face( Qp.allfpoints() ),
89 M_prod( Qp.nPoints() ),
90 M_exprq( Qp.nPoints() )
93 PointSetQuadrature( uint32_type Npoints )
94 : super( Npoints ), M_w( Npoints ), M_prod( Npoints ), M_exprq( Npoints )
97 PointSetQuadrature( weights_type Wts )
101 M_prod( Wts.size() ),
102 M_exprq( Wts.size() )
106 virtual ~PointSetQuadrature()
109 virtual bool isFaceIm()
const
114 self_type& operator=( self_type
const& q )
118 super::operator=( q );
120 M_w_face = q.M_w_face;
121 M_n_face = q.M_n_face;
137 std::vector<nodes_type>
const& allfpoints()
const
146 weights_type
const&
weights( uint16_type __f )
const
148 if ( __f == uint16_type( -1 ) )
151 return M_w_face[__f];
157 nodes_type
const&
fpoints( uint16_type __f )
const
159 if ( __f == uint16_type( -1 ) )
160 return this->points();
162 return M_n_face[__f];
169 value_type
weight( uint16_type __f, uint32_type q )
const
171 return M_w_face[__f][q];
176 ublas::matrix_column<nodes_type const>
fpoint( uint16_type __f, uint32_type __q )
const
178 return ublas::column( M_n_face[__f], __q );
187 weights_type
const& weights()
const
191 value_type
const&
weight(
int q )
const
198 return M_n_face.size();
200 size_type nPointsOnFace(
int __face = 0 )
const
202 return M_n_face[__face].size2();
211 template<
typename Expression>
217 value_type res( 0.0 );
219 for ( uint32_type k = 0; k < this->nPoints(); ++k )
222 value_type fk = f( k );
223 res += this->weights()[k]*fk;
235 template<
typename Expression>
239 value_type res( 0.0 );
241 for ( uint32_type k = 0; k < this->nPoints(); ++k )
242 res += this->weights()[k]*f( this->point( k ) );
253 template<
typename Expression>
255 Expression
const& f )
const
258 FEELPP_ASSERT( __face == ALL_FACES || ( __face >= 0 && __face < nFaces() ) )
259 ( __face )( nFaces() ).error(
"invalid face index (can be ALL_FACE or FACE_0 <= f < nFaces()" );
260 value_type res( 0.0 );
262 if ( __face != ALL_FACES )
265 for ( uint32_type k = 0; k < this->nPointsOnFace( __face ); ++k )
266 res += this->
weight( __face, k )*f( k );
272 for ( uint16_type i = 0; i < this->nFaces(); ++i )
274 value_type __res_face( 0.0 );
276 for ( uint32_type k = 0; k < this->nPointsOnFace( i ); ++k )
277 __res_face += this->
weight( i, k )*f( k );
292 template<
typename Expression>
294 Expression
const& f )
const
296 DVLOG(2) <<
"[PointSetQuadrature] face " << int( __face )<<
" integration\n";
298 FEELPP_ASSERT(
int( __face ) ==
int( ALL_FACES ) ||
299 (
int( __face ) >= 0 &&
300 int( __face ) <
int( nFaces() ) ) )
301 ( __face )( nFaces() ).error(
"invalid face index (can be ALL_FACE or FACE_0 <= f < nFaces()" );
302 value_type res( 0.0 );
304 if ( __face != ALL_FACES )
307 for ( uint32_type k = 0; k < this->nPointsOnFace( __face ); ++k )
308 res += this->
weight( __face, k )*f( this->
fpoint( __face, k ) );
314 for ( uint16_type i = 0; i < this->nFaces(); ++i )
316 value_type __res_face( 0.0 );
318 for ( uint32_type k = 0; k < this->nPointsOnFace( i ); ++k )
319 __res_face += this->
weight( i, k )*f( this->
fpoint( i, k ) );
329 template<
typename IndexTest,
typename IndexTrial,
typename ExprType>
330 value_type operator()( ExprType
const& expr,
331 IndexTest
const& indi,
332 IndexTrial
const& indj,
334 uint16_type c2 )
const
336 for ( uint16_type q = 0; q < this->nPoints(); ++q )
338 M_exprq[q] = expr.evalijq( indi, indj, c1, c2, q );
341 return M_prod.dot( M_exprq );
343 template<
typename IndexTest,
typename ExprType>
344 value_type operator()( ExprType
const& expr,
345 IndexTest
const& indi,
347 uint16_type c2 )
const
349 for ( uint16_type q = 0; q < this->nPoints(); ++q )
351 M_exprq[q] = expr.evaliq( indi, c1, c2, q );
354 return M_prod.dot( M_exprq );
357 template<
typename ExprT>
358 value_type operator()( ExprT
const& expr,
360 uint16_type c2 )
const
362 for ( uint16_type q = 0; q < this->nPoints(); ++q )
364 M_exprq[q] = expr.evalq( c1, c2, q );
367 return M_prod.dot( M_exprq );
369 template<
typename IndexTest,
typename IndexTrial,
typename ExprType>
370 value_type operator()( ExprType
const& expr,
371 IndexTest
const& indi,
372 IndexTrial
const& indj,
375 std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad )
const
377 value_type res = value_type( 0 );
379 for ( uint16_type q = 0; q < indexLocalToQuad.size(); ++q )
381 auto qReal = indexLocalToQuad[q].get<0>();
382 const value_type val_expr = expr.evalijq( indi, indj, c1, c2, q );
383 res += M_prod[qReal]*val_expr;
388 template<
typename IndexTest,
typename ExprType>
389 value_type operator()( ExprType
const& expr,
390 IndexTest
const& indi,
393 std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad )
const
395 value_type res = value_type( 0 );
397 for ( uint16_type q = 0; q < indexLocalToQuad.size(); ++q )
399 auto qReal = indexLocalToQuad[q].get<0>();
400 const value_type val_expr = expr.evaliq( indi, c1, c2, q );
401 res += M_prod[qReal]*val_expr;
407 template<
typename GMC>
408 void update( GMC
const& gmc )
410 if ( this->isFaceIm() )
411 for ( uint16_type q = 0; q < this->nPoints(); ++q )
413 M_prod[q] = M_w( q )*gmc.J( q )*gmc.normalNorm( q );
417 for ( uint16_type q = 0; q < this->nPoints(); ++q )
419 M_prod[q] = M_w( q )*gmc.J( q );
423 template<
typename GMC>
424 void update( GMC
const& gmc,
425 std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad )
428 if ( this->isFaceIm() )
429 for ( uint16_type q = 0; q < indexLocalToQuad.size(); ++q )
431 auto qReal = indexLocalToQuad[q].get<0>();
432 M_prod[qReal] = M_w( qReal )*gmc.J( q )*gmc.normalNorm( q );
436 for ( uint16_type q = 0; q < indexLocalToQuad.size(); ++q )
438 auto qReal = indexLocalToQuad[q].get<0>();
439 M_prod[qReal] = M_w( qReal )*gmc.J( q );
446 public PointSetQuadrature<typename Convex::topological_face_type,
450 typedef PointSetQuadrature<
typename Convex::topological_face_type,
454 typedef super parent_quadrature_type;
455 static const bool is_face_im =
true;
462 Face( self_type
const& quad_elt, uint16_type f = 0 )
467 this->
setPoints( quad_elt.fpoints( f ) );
468 this->setWeights( quad_elt.weights( f ) );
471 Face( Face
const& f )
477 Face& operator=( Face
const& f )
481 super::operator=( f );
487 bool isFaceIm()
const
491 uint16_type face()
const
499 typedef Face face_quadrature_type;
500 face_quadrature_type face( uint16_type f )
const
502 return Face( *
this, f );
505 FEELPP_DEFINE_VISITABLE();
509 void setWeights( weights_type
const& w )
511 M_prod.resize( w.size() );
512 M_exprq.resize( w.size() );
516 template<
typename Elem,
typename GM,
typename IM>
517 void constructQROnFace( Elem
const& ref_convex,
518 boost::shared_ptr<GM>
const& __gm,
519 boost::shared_ptr<IM>
const& __qr_face )
521 constructQROnFace( ref_convex, __gm, __qr_face, mpl::bool_<( Convex::nDim > 1 )>() );
524 template<
typename Elem,
typename GM,
typename IM>
525 void constructQROnFace( Elem
const& ,
526 boost::shared_ptr<GM>
const& ,
527 boost::shared_ptr<IM>
const& ,
530 BOOST_STATIC_ASSERT( Elem::nDim == 1 );
531 M_n_face.resize( Elem::numTopologicalFaces );
532 M_w_face.resize( Elem::numTopologicalFaces );
534 M_n_face[0].resize( Elem::nDim, 1 );
535 M_w_face[0].resize( 1 );
536 M_n_face[0]( 0, 0 ) = -1;
537 M_w_face[0]( 0 ) = 1;
539 M_n_face[1].resize( Elem::nDim, 1 );
540 M_w_face[1].resize( 1 );
541 M_n_face[1]( 0, 0 ) = 1;
542 M_w_face[1]( 0 ) = 1;
546 template<
typename Elem,
typename GM,
typename IM>
547 void constructQROnFace( Elem
const& ref_convex,
548 boost::shared_ptr<GM>
const& __gm,
549 boost::shared_ptr<IM>
const& __qr_face,
556 std::vector<weights_type> M_w_face;
557 std::vector<nodes_type> M_n_face;
560 mutable vector_type M_exprq;
564 template<
class Convex, u
int16_type Integration_Degree,
typename T>
565 template<
typename Elem,
typename GM,
typename IM>
567 PointSetQuadrature<Convex,Integration_Degree,T>::constructQROnFace( Elem
const& ref_convex,
568 boost::shared_ptr<GM>
const& __gm,
569 boost::shared_ptr<IM>
const& __qr_face,
572 M_n_face.resize( Elem::numTopologicalFaces );
573 M_w_face.resize( Elem::numTopologicalFaces );
578 typedef typename GM::face_gm_type::precompute_type face_pc_type;
579 typedef typename GM::face_gm_type::precompute_ptrtype face_pc_ptrtype;
580 face_pc_ptrtype __geopc(
new face_pc_type( __gm->boundaryMap(),__qr_face->points() ) );
582 for ( uint16_type __f = 0; __f < Elem::numTopologicalFaces; ++__f )
584 typedef typename Elem::topological_face_type element_type;
585 element_type ref_convex_face = ref_convex.topologicalFace( __f );
586 DVLOG(2) <<
"[quadpt] ref_convex_face " << __f <<
"=" << ref_convex_face.points() <<
"\n";
589 typename GM::template Context<vm::JACOBIAN|vm::POINT|vm::KB,element_type> __c( __gm->boundaryMap(),
592 __c.update( ref_convex_face );
593 DVLOG(2) <<
"[quadpt] ref_convex_face " << __f <<
" xref" << __c.xRefs() <<
"\n";
594 DVLOG(2) <<
"[quadpt] ref_convex_face " << __f <<
" xreal" << __c.xReal() <<
"\n";
596 value_type __len = 0.0;
597 M_n_face[__f].resize( Elem::nDim, __qr_face->nPoints() );
598 M_w_face[__f].resize( __qr_face->nPoints() );
599 DVLOG(2) <<
"[PointSetQuadrature::constructQROnFace] npoints on face "
601 << __qr_face->nPoints() <<
"\n";
605 for ( uint16_type __ip = 0; __ip < __qr_face->nPoints(); ++__ip )
607 ublas::column( M_n_face[__f], __ip ) = __c.xReal( __ip );
610 M_w_face[ __f]( __ip ) = __qr_face->weight( __ip )*__c.J( __ip );
612 __len += M_w_face[ __f]( __ip );
613 DVLOG(2) <<
"face " << __f <<
" ip = " << __ip <<
" J =" << __c.J( __ip ) <<
"\n";
614 DVLOG(2) <<
"face " << __f <<
" ip = " << __ip <<
" weight =" << __qr_face->weight( __ip ) <<
"\n";
615 DVLOG(2) <<
"face " << __f <<
" ip = " << __ip <<
" weight =" << M_w_face[ __f]( __ip ) <<
"\n";
620 DVLOG(2) <<
"length = " << __len <<
"\n";
630 template<
class Convex, Feel::u
int16_type Integration_Degree,
typename T>
631 std::ostream& operator<<( std::ostream& os,
634 os <<
"quadrature point set:\n"
635 <<
"number of points: " << qps.nPoints() <<
"\n"
636 <<
"points : " << qps.points() <<
"\n"
637 <<
"weights : " << qps.
weights() <<
"\n";
639 for ( Feel::uint16_type f = 0; f < qps.nFaces(); ++f )
641 os <<
" o Face " << f <<
"\n";
642 os << qps.
fpoints( f ) <<
"\n";
643 os << qps.
weights( f ) <<
"\n";
value_type integrateAtPoints(Expression const &f) const
Definition: quadpoint.hpp:236
const uint16_type invalid_uint16_type_value
Definition: feelcore/feel.hpp:338
value_type integrate(IntegrationFaceEnum __face, Expression const &f) const
Definition: quadpoint.hpp:254
value_type weight(uint16_type __f, uint32_type q) const
Definition: quadpoint.hpp:169
Quadrature point set base class.
Definition: gausslobatto.hpp:58
std::vector< weights_type > const & allfweights() const
Definition: quadpoint.hpp:132
void setPoints(nodes_type const &pts)
Definition: pointset.hpp:267
weights_type const & weights(uint16_type __f) const
Definition: quadpoint.hpp:146
value_type integrate(Expression f) const
Definition: quadpoint.hpp:212
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
value_type integrateAtPoints(IntegrationFaceEnum __face, Expression const &f) const
Definition: quadpoint.hpp:293
nodes_type const & fpoints(uint16_type __f) const
Definition: quadpoint.hpp:157
ublas::matrix_column< nodes_type const > fpoint(uint16_type __f, uint32_type __q) const
Definition: quadpoint.hpp:176