Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
quadpoint.hpp
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Gilles Steiner <gilles.steiner@epfl.ch>
6  Christophe Prud'homme <christophe.prudhomme@feelpp.org>
7  Date: 2005-07-28
8 
9  Copyright (C) 2005,2006 EPFL
10  Copyright (C) 2006-2010 Universite Joseph Fourier
11 
12  This library is free software; you can redistribute it and/or
13  modify it under the terms of the GNU Lesser General Public
14  License as published by the Free Software Foundation; either
15  version 3.0 of the License, or (at your option) any later version.
16 
17  This library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20  Lesser General Public License for more details.
21 
22  You should have received a copy of the GNU Lesser General Public
23  License along with this library; if not, write to the Free Software
24  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25 */
31 #ifndef __Quadpoint_H
32 #define __Quadpoint_H 1
33 
35 #include <feel/feelpoly/jacobi.hpp>
36 #include <feel/feelpoly/geomap.hpp>
37 
38 namespace Feel
39 {
40 namespace ublas = boost::numeric::ublas;
41 
42 enum IntegrationFaceEnum
43 {
44  ALL_FACES = -1,
45  FACE_0 = 0,
46  FACE_1 = 1,
47  FACE_2,
48  FACE_3,
49  FACE_4,
50  FACE_5,
51  FACE_6,
52  FACE_7,
53  FACE_8
54 };
55 
64 template<class Convex, uint16_type Integration_Degree, typename T>
65 class PointSetQuadrature : public PointSet<Convex,T>
66 {
67 public :
68  static const bool is_face_im = false;
69  typedef T value_type;
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;
76 
77  typedef PointSetQuadrature<Convex, Integration_Degree, T> self_type;
78 
79  typedef self_type parent_quadrature_type;
80  static const uint16_type I_deg = Integration_Degree;
81 
82  PointSetQuadrature(): super(), M_w(), M_prod(), M_exprq() {}
83 
84  PointSetQuadrature( const PointSetQuadrature& Qp )
85  : super( Qp ),
86  M_w( Qp.weights() ),
87  M_w_face( Qp.allfweights() ),
88  M_n_face( Qp.allfpoints() ),
89  M_prod( Qp.nPoints() ),
90  M_exprq( Qp.nPoints() )
91  {}
92 
93  PointSetQuadrature( uint32_type Npoints )
94  : super( Npoints ), M_w( Npoints ), M_prod( Npoints ), M_exprq( Npoints )
95  {}
96 
97  PointSetQuadrature( weights_type Wts )
98  :
99  super( Wts.size() ),
100  M_w( Wts ),
101  M_prod( Wts.size() ),
102  M_exprq( Wts.size() )
103  {}
104 
105 
106  virtual ~PointSetQuadrature()
107  {}
108 
109  virtual bool isFaceIm() const
110  {
111  return is_face_im;
112  }
113 
114  self_type& operator=( self_type const& q )
115  {
116  if ( this != &q )
117  {
118  super::operator=( q );
119  M_w = q.M_w;
120  M_w_face = q.M_w_face;
121  M_n_face = q.M_n_face;
122  M_prod = q.M_prod;
123  M_exprq = q.M_exprq;
124  }
125 
126  return *this;
127  }
128 
132  std::vector<weights_type> const& allfweights() const
133  {
134  return M_w_face;
135  }
136 
137  std::vector<nodes_type> const& allfpoints() const
138  {
139  return M_n_face;
140  }
141 
142 
146  weights_type const& weights( uint16_type __f ) const
147  {
148  if ( __f == uint16_type( -1 ) )
149  return M_w;
150 
151  return M_w_face[__f];
152  }
153 
157  nodes_type const& fpoints( uint16_type __f ) const
158  {
159  if ( __f == uint16_type( -1 ) )
160  return this->points();
161 
162  return M_n_face[__f];
163  }
164 
165 
169  value_type weight( uint16_type __f, uint32_type q ) const
170  {
171  return M_w_face[__f][q];
172  }
176  ublas::matrix_column<nodes_type const> fpoint( uint16_type __f, uint32_type __q ) const
177  {
178  return ublas::column( M_n_face[__f], __q );
179  }
180 
181 
182  // quadrature_data_type data() const { return boost::make_tuple( this->M_n, this->M_w ); }
183 
184  //size_type nFaces() const { return M_n_face.size(); }
185  //size_type nPointsOnFace( int __face = 0 ) const { return M_n_face[__face].size2(); }
186 
187  weights_type const& weights() const
188  {
189  return M_w;
190  }
191  value_type const& weight( int q ) const
192  {
193  return M_w[q];
194  }
195 
196  size_type nFaces() const
197  {
198  return M_n_face.size();
199  }
200  size_type nPointsOnFace( int __face = 0 ) const
201  {
202  return M_n_face[__face].size2();
203  }
204 
211  template<typename Expression>
212  value_type integrate( Expression f ) const
213  {
214  //BOOST_STATIC_ASSERT( ( boost::is_same<value_type,qd_real>::value ) );
215  //BOOST_MPL_ASSERT_MSG( ( boost::is_same<value_type,qd_real>::value ), INVALID_TYPE, (value_type) );
216  //node_type res( this->weights()[0]*f( this->point(0) ) );
217  value_type res( 0.0 );
218 
219  for ( uint32_type k = 0; k < this->nPoints(); ++k )
220  {
221 
222  value_type fk = f( k );
223  res += this->weights()[k]*fk;
224  }
225 
226  return res;
227  }
228 
235  template<typename Expression>
236  value_type integrateAtPoints( Expression const& f ) const
237  {
238  //node_type res( this->weights()[0]*f( this->point(0) ) );
239  value_type res( 0.0 );
240 
241  for ( uint32_type k = 0; k < this->nPoints(); ++k )
242  res += this->weights()[k]*f( this->point( k ) );
243 
244  return res;
245  }
246 
253  template<typename Expression>
254  value_type integrate( IntegrationFaceEnum __face,
255  Expression const& f ) const
256  {
257  //std::cout << "integrating using " << nNodes << "\n";
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 );
261 
262  if ( __face != ALL_FACES )
263  {
264  // integrate on face f
265  for ( uint32_type k = 0; k < this->nPointsOnFace( __face ); ++k )
266  res += this->weight( __face, k )*f( k );
267  }
268 
269  else
270  {
271  // integrate on all faces
272  for ( uint16_type i = 0; i < this->nFaces(); ++i )
273  {
274  value_type __res_face( 0.0 );
275 
276  for ( uint32_type k = 0; k < this->nPointsOnFace( i ); ++k )
277  __res_face += this->weight( i, k )*f( k );
278 
279  //std::cout << "res on face " << i << " = " << __res_face << "\n";
280  res += __res_face;
281  }
282  }
283 
284  return res;
285  }
292  template<typename Expression>
293  value_type integrateAtPoints( IntegrationFaceEnum __face,
294  Expression const& f ) const
295  {
296  DVLOG(2) << "[PointSetQuadrature] face " << int( __face )<< " integration\n";
297  //std::cout << "integrating using " << nNodes << "\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 );
303 
304  if ( __face != ALL_FACES )
305  {
306  // integrate on face f
307  for ( uint32_type k = 0; k < this->nPointsOnFace( __face ); ++k )
308  res += this->weight( __face, k )*f( this->fpoint( __face, k ) );
309  }
310 
311  else
312  {
313  // integrate on all faces
314  for ( uint16_type i = 0; i < this->nFaces(); ++i )
315  {
316  value_type __res_face( 0.0 );
317 
318  for ( uint32_type k = 0; k < this->nPointsOnFace( i ); ++k )
319  __res_face += this->weight( i, k )*f( this->fpoint( i, k ) );
320 
321  //std::cout << "res on face " << i << " = " << __res_face << "\n";
322  res += __res_face;
323  }
324  }
325 
326  return res;
327  }
328 
329  template<typename IndexTest, typename IndexTrial, typename ExprType>
330  value_type operator()( ExprType const& expr,
331  IndexTest const& indi,
332  IndexTrial const& indj,
333  uint16_type c1,
334  uint16_type c2 ) const
335  {
336  for ( uint16_type q = 0; q < this->nPoints(); ++q )
337  {
338  M_exprq[q] = expr.evalijq( indi, indj, c1, c2, q );
339  }
340 
341  return M_prod.dot( M_exprq );
342  }
343  template<typename IndexTest, typename ExprType>
344  value_type operator()( ExprType const& expr,
345  IndexTest const& indi,
346  uint16_type c1,
347  uint16_type c2 ) const
348  {
349  for ( uint16_type q = 0; q < this->nPoints(); ++q )
350  {
351  M_exprq[q] = expr.evaliq( indi, c1, c2, q );
352  }
353 
354  return M_prod.dot( M_exprq );
355  }
356 
357  template<typename ExprT>
358  value_type operator()( ExprT const& expr,
359  uint16_type c1,
360  uint16_type c2 ) const
361  {
362  for ( uint16_type q = 0; q < this->nPoints(); ++q )
363  {
364  M_exprq[q] = expr.evalq( c1, c2, q );
365  }
366 
367  return M_prod.dot( M_exprq );
368  }
369  template<typename IndexTest, typename IndexTrial, typename ExprType>
370  value_type operator()( ExprType const& expr,
371  IndexTest const& indi,
372  IndexTrial const& indj,
373  uint16_type c1,
374  uint16_type c2,
375  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad ) const
376  {
377  value_type res = value_type( 0 );
378 
379  for ( uint16_type q = 0; q < indexLocalToQuad.size(); ++q )
380  {
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;
384  }
385 
386  return res;
387  }
388  template<typename IndexTest, typename ExprType>
389  value_type operator()( ExprType const& expr,
390  IndexTest const& indi,
391  uint16_type c1,
392  uint16_type c2,
393  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad ) const
394  {
395  value_type res = value_type( 0 );
396 
397  for ( uint16_type q = 0; q < indexLocalToQuad.size(); ++q )
398  {
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;
402  }
403 
404  return res;
405  }
406 
407  template<typename GMC>
408  void update( GMC const& gmc )
409  {
410  if ( this->isFaceIm() )
411  for ( uint16_type q = 0; q < this->nPoints(); ++q )
412  {
413  M_prod[q] = M_w( q )*gmc.J( q )*gmc.normalNorm( q );
414  }
415 
416  else
417  for ( uint16_type q = 0; q < this->nPoints(); ++q )
418  {
419  M_prod[q] = M_w( q )*gmc.J( q );
420  }
421  }
422 
423  template<typename GMC>
424  void update( GMC const& gmc,
425  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad )
426  {
427 
428  if ( this->isFaceIm() )
429  for ( uint16_type q = 0; q < indexLocalToQuad.size(); ++q )
430  {
431  auto qReal = indexLocalToQuad[q].get<0>();
432  M_prod[qReal] = M_w( qReal )*gmc.J( q )*gmc.normalNorm( q );
433  }
434 
435  else
436  for ( uint16_type q = 0; q < indexLocalToQuad.size(); ++q )
437  {
438  auto qReal = indexLocalToQuad[q].get<0>();
439  M_prod[qReal] = M_w( qReal )*gmc.J( q );
440  }
441  }
442 
443 
444  class Face
445  :
446  public PointSetQuadrature<typename Convex::topological_face_type,
447  Integration_Degree,
448  T>
449  {
450  typedef PointSetQuadrature<typename Convex::topological_face_type,
451  Integration_Degree,
452  T> super;
453  public:
454  typedef super parent_quadrature_type;
455  static const bool is_face_im = true;
456  Face()
457  :
458  super(),
460  {
461  }
462  Face( self_type const& quad_elt, uint16_type f = 0 )
463  :
464  super(),
465  M_f( f )
466  {
467  this->setPoints( quad_elt.fpoints( f ) );
468  this->setWeights( quad_elt.weights( f ) );
469 
470  }
471  Face( Face const& f )
472  :
473  super( f ),
474  M_f( f.M_f )
475 
476  {}
477  Face& operator=( Face const& f )
478  {
479  if ( this != &f )
480  {
481  super::operator=( f );
482  M_f = f.M_f;
483  }
484 
485  return *this;
486  }
487  bool isFaceIm() const
488  {
489  return is_face_im;
490  }
491  uint16_type face() const
492  {
493  return M_f;
494  }
495  private:
496  uint16_type M_f;
497  };
498 
499  typedef Face face_quadrature_type;
500  face_quadrature_type face( uint16_type f ) const
501  {
502  return Face( *this, f );
503  }
504 
505  FEELPP_DEFINE_VISITABLE();
506 
507 protected:
508 
509  void setWeights( weights_type const& w )
510  {
511  M_prod.resize( w.size() );
512  M_exprq.resize( w.size() );
513  M_w = w;
514  }
515 
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 )
520  {
521  constructQROnFace( ref_convex, __gm, __qr_face, mpl::bool_<( Convex::nDim > 1 )>() );
522  }
523 
524  template<typename Elem, typename GM, typename IM>
525  void constructQROnFace( Elem const& /*ref_convex*/,
526  boost::shared_ptr<GM> const& /*__gm*/,
527  boost::shared_ptr<IM> const& /*__qr_face*/,
528  mpl::bool_<false> )
529  {
530  BOOST_STATIC_ASSERT( Elem::nDim == 1 );
531  M_n_face.resize( Elem::numTopologicalFaces );
532  M_w_face.resize( Elem::numTopologicalFaces );
533 
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;
538 
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;
543 
544  }
545 
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,
550  mpl::bool_<true> );
551 
552 protected:
553 
554  weights_type M_w;
555 
556  std::vector<weights_type> M_w_face;
557  std::vector<nodes_type> M_n_face;
558 
559  vector_type M_prod;
560  mutable vector_type M_exprq;
561 
562 };
563 
564 template<class Convex, uint16_type Integration_Degree, typename T>
565 template<typename Elem, typename GM, typename IM>
566 void
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,
570  mpl::bool_<true> )
571 {
572  M_n_face.resize( Elem::numTopologicalFaces );
573  M_w_face.resize( Elem::numTopologicalFaces );
574 
575  // FIXME: we don't handle the case where faces are not of the
576  // same type like for prism or pyramids. In that case geopc
577  // should be defined per face
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() ) );
581 
582  for ( uint16_type __f = 0; __f < Elem::numTopologicalFaces; ++__f )
583  {
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";
587  //toPython( ref_convex_face );
588 
589  typename GM::template Context<vm::JACOBIAN|vm::POINT|vm::KB,element_type> __c( __gm->boundaryMap(),
590  ref_convex_face,
591  __geopc );
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";
595 
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 "
600  << __f << " : "
601  << __qr_face->nPoints() << "\n";
602  // transform the quad nodes on the boundary _reference_
603  // element to the face of the reference element face
604 
605  for ( uint16_type __ip = 0; __ip < __qr_face->nPoints(); ++__ip )
606  {
607  ublas::column( M_n_face[__f], __ip ) = __c.xReal( __ip );
608 
609  // w = w_face * ||B*n|| * J
610  M_w_face[ __f]( __ip ) = __qr_face->weight( __ip )*__c.J( __ip );
611 
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";
616  // DVLOG(2) << "face " << __f << " ip = " << __ip << " x ref =" << __c.xRef( __ip ) << "\n";
617  // DVLOG(2) << "face " << __f << " ip = " << __ip << " x real =" << __c.xReal( __ip ) << "\n";
618  }
619 
620  DVLOG(2) << "length = " << __len << "\n";
621  }
622 }
623 
624 
625 } // Feel
626 
627 
628 namespace std
629 {
630 template<class Convex, Feel::uint16_type Integration_Degree, typename T>
631 std::ostream& operator<<( std::ostream& os,
633 {
634  os << "quadrature point set:\n"
635  << "number of points: " << qps.nPoints() << "\n"
636  << "points : " << qps.points() << "\n"
637  << "weights : " << qps.weights() << "\n";
638 
639  for ( Feel::uint16_type f = 0; f < qps.nFaces(); ++f )
640  {
641  os << " o Face " << f << "\n";
642  os << qps.fpoints( f ) << "\n";
643  os << qps.weights( f ) << "\n";
644  }
645 
646  return os;
647 }
648 
649 }
650 
651 #endif /* _QuadPoint_H */

Generated on Sun Oct 20 2013 08:25:04 for Feel++ by doxygen 1.8.4