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
gauss.hpp
Go to the documentation of this file.
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): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2006-12-30
7 
8  Copyright (C) 2006 Universite Joseph Fourier (Grenoble)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __Gauss_H
30 #define __Gauss_H 1
31 
32 #include <feel/feelpoly/quadpoint.hpp>
33 
34 namespace Feel
35 {
50 template<class Convex, uint16_type Integration_Degree, typename T>
51 class Gauss : public PointSetQuadrature<Convex, Integration_Degree, T> {};
52 
53 template< uint16_type Integration_Degree, typename T>
54 class Gauss<Simplex<0,1> , Integration_Degree ,T > : public PointSetQuadrature<Simplex<0,1> , Integration_Degree, T>
55 {
56 public :
57  typedef T value_type;
58 
59  typedef PointSetQuadrature<Simplex<0,1> , Integration_Degree, T> super;
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;
64 
65  static const uint16_type Degree = invalid_uint16_type_value;
66  static const uint32_type Npoints = 1;
67 
68  Gauss()
69  :
70  super( Npoints )
71  {
72 
73  }
74 
75  ~Gauss() {}
76 
77  FEELPP_DEFINE_VISITABLE();
78 };
79 
81 template< uint16_type Integration_Degree, typename T>
82 class Gauss<Simplex<1,1> , Integration_Degree ,T > : public PointSetQuadrature<Simplex<1,1> , Integration_Degree, T>
83 {
84 public :
85  typedef T value_type;
86 
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;
92 
93  static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
94  static const uint32_type Npoints = Degree;
95 
96  typedef Gauss<Simplex<0,1>,Integration_Degree, T> face_quad_type;
97 
98  Gauss()
99  :
100  super( Npoints )
101  {
102  ublas::vector<T> px( Npoints );
103 
104  details::gaussjacobi<Npoints, T, ublas::vector<T>, ublas::vector<T> >( this->M_w, px );
105  ublas::row( this->M_points, 0 ) = px;
106 
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 );
109  // construct face quadratures
110  this->constructQROnFace( Reference<Simplex<1, 1>, 1, 1>(), gm, face_qr );
111 
112 
113  }
114 
115  ~Gauss() {}
116 
117  FEELPP_DEFINE_VISITABLE();
118 };
119 
122 template< uint16_type Integration_Degree, typename T>
123 class Gauss<Simplex<2,1> , Integration_Degree ,T > : public PointSetQuadrature<Simplex<2,1> , Integration_Degree, T>
124 {
125 public :
126  typedef T value_type;
127 
128  typedef PointSetQuadrature<Simplex<2,1> , Integration_Degree, T> super;
129  typedef typename super::return_type return_type;
130 
131  typedef typename super::node_type node_type;
132  typedef typename super::nodes_type nodes_type;
133  typedef typename super::weights_type weights_type;
134 
135  typedef Gauss<Simplex<1,1>,Integration_Degree, T> face_quad_type;
136 
137 
138  static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
139  static const uint32_type Npoints = Degree*Degree;
140 
141  Gauss()
142  :
143  super( Npoints )
144  {
145  // build rules in x and y direction
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 );
149 
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 );
153 
154  // coordinate in cartesian space
155 
156 #if 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 ;
160 #endif
161 
162  node_type eta( 2 );
163  details::xi<TRIANGLE, value_type> to_xi;
164 
165  for ( int i = 0, k = 0; i < Degree; ++i )
166  {
167  for ( int j = 0; j < Degree; ++j, ++k )
168  {
169 #if 0
170 
171  if ( j%100==0 )
172  std::cout<<"[Debug quadpoint] i = " << i << " ; j = " << j << " ; k = " << k << std::endl ;
173 
174 #endif
175 
176  // computes the weight of the k-th node
177  this->M_w( k ) = 0.5 * wx( i ) * wy( j );
178  // use expansion for the collapsed triangle to compute the points
179  // coordinates (from cartesian to collapsed coordinates)
180  eta( 0 ) = px( i );
181  eta( 1 ) = py( j );
182 
183  ublas::column( this->M_points, k ) = to_xi( eta );
184  }
185  }
186 
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 );
189  // construct face quadratures
190  this->constructQROnFace( Reference<Simplex<2, 1>, 2, 1>(), gm, face_qr );
191 
192  }
193 
194  ~Gauss() {}
195 
196  FEELPP_DEFINE_VISITABLE();
197 };
198 
199 
202 template< uint16_type Integration_Degree, typename T>
203 class Gauss<Simplex<3,1> , Integration_Degree ,T > : public PointSetQuadrature<Simplex<3,1> , Integration_Degree, T>
204 {
205 public :
206  typedef T value_type;
207 
208  typedef PointSetQuadrature<Simplex<3,1> , Integration_Degree, T> super;
209 
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;
214 
215  typedef Gauss<Simplex<2,1>,Integration_Degree, T> face_quad_type;
216 
217 
218  static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
219  static const uint32_type Npoints = Degree*Degree*Degree;
220 
221  Gauss()
222  :
223  super( Npoints )
224  {
225  // build rules in x and y direction
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 );
229 
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 );
233 
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 );
237 
238  // coordinate in cartesian space
239  node_type eta( 3 );
240  details::xi<TETRAHEDRON, value_type> to_xi;
241 
242  for ( int i = 0, k = 0; i < Degree; ++i )
243  {
244  for ( int j = 0; j < Degree; ++j )
245  {
246  for ( int l = 0; l < Degree; ++l, ++k )
247  {
248  // computes the weight of the k-th node
249  this->M_w( k ) = 0.125 * wx( i ) * wy( j ) * wz( l );
250  // use expansion for the collapsed triangle to compute the points
251  // coordinates (from cartesian to collapsed coordinates)
252  eta( 0 ) = px( i );
253  eta( 1 ) = py( j );
254  eta( 2 ) = pz( l );
255  ublas::column( this->M_points, k ) = to_xi( eta );
256  }
257  }
258  }
259 
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 );
262  // construct face quadratures
263  this->constructQROnFace( Reference<Simplex<3, 1>,3,1>(), gm, face_qr );
264  }
265 
266  ~Gauss() {}
267 
268  FEELPP_DEFINE_VISITABLE();
269 };
270 
273 template< uint16_type Integration_Degree, typename T>
274 class Gauss<Hypercube<1,1>, Integration_Degree ,T >
275  :
276 public PointSetQuadrature<Hypercube<1,1>, Integration_Degree, T>
277 {
278 public :
279  typedef T value_type;
280 
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;
286 
287  static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
288  static const uint32_type Npoints = Degree;
289 
290  typedef Gauss<Simplex<0,1>,Integration_Degree, T> face_quad_type;
291 
292  Gauss( )
293  :
294  super( Npoints )
295  {
296  // build rules in x and y direction
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 );
300 #if 0
301  VLOG(1) << "[gauss<SP<2,1>] jacobi p = " << px << "\n";
302  VLOG(1) << "[gauss<SP<2,1>] jacobi w = " << wx << "\n";
303 #endif
304 
305  for ( int i = 0; i < Degree; ++i )
306  {
307  // computes the weight of the k-th node
308  this->M_w( i ) = wx( i );
309  this->M_points( 0, i ) = px( i );
310  }
311 
312 
313 #if 0
314  VLOG(1) << "[gauss<SP<2,1>] p = " << this->M_points << "\n";
315  VLOG(1) << "[gauss<SP<2,1>] w = " << this->M_w << "\n";
316 #endif
317 
318 
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 );
321  // construct face quadratures
322  this->constructQROnFace( Reference<Hypercube<1, 1>, 1, 1>(), gm, face_qr );
323 
324  }
325 
326  ~Gauss() {}
327 
328  FEELPP_DEFINE_VISITABLE();
329 };
332 template< uint16_type Integration_Degree, typename T>
333 class Gauss<Hypercube<2,1>, Integration_Degree ,T >
334  :
335 public PointSetQuadrature<Hypercube<2,1>, Integration_Degree, T>
336 {
337 public :
338  typedef T value_type;
339 
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;
345 
346  typedef Gauss<Hypercube<1,1>,Integration_Degree, T> face_quad_type;
347 
348  static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
349  static const uint32_type Npoints = Degree*Degree;
350 
351  Gauss()
352  :
353  super( Npoints )
354  {
355  // build rules in x and y direction
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 );
359 #if 0
360  VLOG(1) << "[gauss<SP<2,1>] jacobi p = " << px << "\n";
361  VLOG(1) << "[gauss<SP<2,1>] jacobi w = " << wx << "\n";
362 #endif
363 
364  for ( int i = 0, k = 0; i < Degree; ++i )
365  {
366  for ( int j = 0; j < Degree; ++j, ++k )
367  {
368  // computes the weight of the k-th node
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 );
372  }
373  }
374 
375 #if 0
376  VLOG(1) << "[gauss<SP<2,1>] p = " << this->M_points << "\n";
377  VLOG(1) << "[gauss<SP<2,1>] w = " << this->M_w << "\n";
378 #endif
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 );
381  // construct face quadratures
382  this->constructQROnFace( Reference<Hypercube<2, 1>,2,1>(), gm, face_qr );
383 
384  }
385 
386  ~Gauss() {}
387 
388  FEELPP_DEFINE_VISITABLE();
389 };
390 
393 template< uint16_type Integration_Degree, typename T>
394 class Gauss<Hypercube<3,1>, Integration_Degree ,T >
395  :
396 public PointSetQuadrature<Hypercube<3,1>, Integration_Degree, T>
397 {
398 public :
399  typedef T value_type;
400 
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;
409 
410  Gauss()
411  :
412  super( Npoints )
413  {
414  // build rules in x and y direction
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 );
418 
419  for ( int i = 0, k = 0; i < Degree; ++i )
420  {
421  for ( int j = 0; j < Degree; ++j )
422  {
423  for ( int l = 0; l < Degree ; ++l, ++k )
424  {
425  // computes the weight of the k-th node
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 );
430  }
431  }
432  }
433 
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 );
436  // construct face quadratures
437  this->constructQROnFace( Reference<Hypercube<3, 1>,3,1>(), gm, face_qr );
438 
439  }
440 
441  ~Gauss() {}
442 
443  FEELPP_DEFINE_VISITABLE();
444 };
445 
446 
447 template< uint16_type Integration_Degree, typename T>
448 class Gauss<Hypercube<4,1>, Integration_Degree ,T >
449  :
450 public PointSetQuadrature<Hypercube<4,1>, Integration_Degree, T>
451 {
452 public :
453  typedef T value_type;
454 
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;
460 
461  static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
462  static const uint32_type Npoints = Degree*Degree*Degree*Degree;
463 
464  Gauss()
465  :
466  super( Npoints )
467  {
468  // build rules in x and y direction
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 );
472 
473  for ( int i = 0, k = 0; i < Degree; ++i )
474  {
475  for ( int j = 0; j < Degree; ++j )
476  {
477  for ( int l = 0; l < Degree ; ++l )
478  {
479  for ( int r = 0; r < Degree ; ++r, ++k )
480  {
481  // computes the weight of the k-th node
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 );
487  }
488  }
489  }
490  }
491  }
492 
493  ~Gauss() {}
494 
495  FEELPP_DEFINE_VISITABLE();
496 };
497 
498 
499 template< uint16_type Integration_Degree, typename T>
500 class Gauss<Hypercube<5,1>, Integration_Degree ,T >
501  :
502 public PointSetQuadrature<Hypercube<5,1>, Integration_Degree, T>
503 {
504 public :
505  typedef T value_type;
506 
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;
512 
513  static const uint16_type Degree = ( Integration_Degree+1 )/2+1;
514  static const uint32_type Npoints = Degree*Degree*Degree*Degree*Degree;
515 
516  Gauss()
517  :
518  super( Npoints )
519  {
520  // build rules in x and y direction
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 );
524 
525  for ( int i = 0, k = 0; i < Degree; ++i )
526  {
527  for ( int j = 0; j < Degree; ++j )
528  {
529  for ( int l = 0; l < Degree ; ++l )
530  {
531  for ( int r = 0; r < Degree ; ++r )
532  {
533  for ( int s = 0; s < Degree ; ++s, ++k )
534  {
535  // computes the weight of the k-th node
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 );
542  }
543  }
544  }
545  }
546  }
547  }
548 
549  ~Gauss() {}
550 
551  FEELPP_DEFINE_VISITABLE();
552 };
554 } // Feel
555 
556 #endif /* __Gauss_H */
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

Generated on Sun Dec 22 2013 13:11:03 for Feel++ by doxygen 1.8.5