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
imsimplex2.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): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2006-06-23
7 
8  Copyright (C) 2006 EPFL
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 __IMSimplex_H
30 #define __IMSimplex_H 1
31 
32 #include <vector>
33 
34 #include <boost/assign/list_of.hpp>
35 #include <boost/assign/std/vector.hpp>
36 #include <boost/numeric/ublas/matrix.hpp>
37 #include <boost/numeric/ublas/matrix_proxy.hpp>
38 #include <boost/numeric/ublas/vector.hpp>
39 #include <boost/numeric/ublas/io.hpp>
40 
41 #include <feel/feelcore/feel.hpp>
42 #include <feel/feelpoly/gauss.hpp>
43 
44 
45 namespace Feel
46 {
47 namespace ublas = boost::numeric::ublas;
48 
50 
54 namespace detail
55 {
56 using namespace boost::assign;
57 
62 template<int Order, typename T>
63 struct IMTriangle
64 {
65 };
66 
67 template<typename T>
68 struct IMTriangle<1,T>
69 {
70  typedef T value_type;
71  static const uint16_type nDim = 2;
72  static const uint16_type nOrder = 1;
73  static const uint16_type nPoints = 1;
74 
75  IMTriangle()
76  :
77  q ( 1 )
78  {
79  m += 1;
80  w += 1;
81  push_back( q[0] ).repeat( 3, value_type( 1 )/value_type( 3 ) );
82  }
83 
84 
85  std::vector<uint16_type> m;
86  std::vector<value_type> w;
87  std::vector<std::vector<value_type> > q;
88 };
89 
90 template<typename T> struct IMTriangle<0,T>: public IMTriangle<1,T> {};
91 
92 template<typename T>
93 struct IMTriangle<2,T>
94 {
95  typedef T value_type;
96  static const uint16_type nDim = 2;
97  static const uint16_type nOrder = 2;
98  static const uint16_type nPoints = 3;
99 
100  IMTriangle()
101  :
102  q( 1 )
103  {
104  m += 3;
105  w += value_type( 1 )/value_type( 3 );
106  q[0] += value_type( 2 )/value_type( 3 ), repeat( 2, value_type( 1 )/value_type( 6 ) );
107  }
108 
109 
110  std::vector<uint16_type> m;
111  std::vector<value_type> w;
112  std::vector<std::vector<value_type> > q;
113 };
114 
115 template<>
116 struct IMTriangle<3,double>
117 {
118  typedef double value_type;
119  static const uint16_type nDim = 2;
120  static const uint16_type nOrder = 3;
121  static const uint16_type nPoints = 6;
122  IMTriangle()
123  :
124  q( 1 )
125  {
126  m += 6;
127  w += value_type( 1 )/value_type( 6 );
128  q[0] += value_type( 0.659027622374092 ), value_type( 0.231933368553031 ), value_type( 0.109039009072877 );
129  }
130 
131  std::vector<uint16_type> m;
132  std::vector<value_type> w;
133  std::vector<std::vector<value_type> > q;
134 };
135 
136 template<>
137 struct IMTriangle<4,double>
138 {
139  typedef double value_type;
140  static const uint16_type nDim = 2;
141  static const uint16_type nOrder = 4;
142  static const uint16_type nPoints = 6;
143  IMTriangle()
144  :
145  q( 2 )
146  {
147  m += 3,3;
148  w += 0.109951743655322, 0.223381589678011;
149  q[0] += 0.816847572980459, 0.091576213509771, 0.091576213509771;
150  q[1] += 0.108103018168070, 0.445948490915965, 0.445948490915965;
151  }
152 
153  std::vector<uint16_type> m;
154  std::vector<value_type> w;
155  std::vector<std::vector<value_type> > q;
156 };
157 
158 template<>
159 struct IMTriangle<5,double>
160 {
161  typedef double value_type;
162  static const uint16_type nDim = 2;
163  static const uint16_type nOrder = 5;
164  static const uint16_type nPoints = 7;
165  IMTriangle()
166  :
167  q( 3 )
168  {
169  m += 1, 3, 3;
170  w += 0.22500000000000, 0.125939180544827, 0.132394152788506;
171  q[0] += value_type( 1.0 )/value_type( 3.0 ), value_type( 1.0 )/value_type( 3.0 ), value_type( 1.0 )/value_type( 3.0 );
172  q[1] += 0.797426985353087, 0.101286507323456, 0.101286507323456;
173  q[2] += 0.059715871789770, 0.470142064105115, 0.470142064105115;
174  }
175 
176  std::vector<uint16_type> m;
177  std::vector<value_type> w;
178  std::vector<std::vector<value_type> > q;
179 };
180 
181 template<>
182 struct IMTriangle<6,double>
183 {
184  typedef double value_type;
185  static const uint16_type nDim = 2;
186  static const uint16_type nOrder = 6;
187  static const uint16_type nPoints = 12;
188  IMTriangle()
189  :
190  q( 3 )
191  {
192  m += 3, 3, 6;
193  w += 0.050844906370207, 0.116786275726379, 0.082851075618374;
194  q[0] += 0.873821971016996, 0.063089014491502, 0.063089014491502;
195  q[1] += 0.501426509658179, 0.249286745170910, 0.249286745170910;
196  q[2] += 0.636502499121399, 0.310352451033784, 0.053145049844817;
197  }
198 
199  std::vector<uint16_type> m;
200  std::vector<value_type> w;
201  std::vector<std::vector<value_type> > q;
202 };
203 
205 
210 template<int Order, typename T>
211 struct IMTetrahedra
212 {
213 };
214 
215 template<typename T>
216 struct IMTetrahedra<1,T>
217 {
218  typedef T value_type;
219  static const uint16_type nDim = 3;
220  static const uint16_type nOrder = 1;
221  static const uint16_type nPoints = 1;
222 
223  IMTetrahedra()
224  :
225  q ( 1 )
226  {
227  m += 1;
228  w += 1;
229  q[0] += value_type( 1 )/value_type( 4 );
230  }
231 
232  std::vector<uint16_type> m;
233  std::vector<value_type> w;
234  std::vector<std::vector<value_type> > q;
235 };
236 template<typename T> struct IMTetrahedra<0,T>: public IMTetrahedra<1,T> {};
237 template<typename T>
238 struct IMTetrahedra<2,T>
239 {
240  typedef T value_type;
241  static const uint16_type nDim = 3;
242  static const uint16_type nOrder = 2;
243  static const uint16_type nPoints = 4;
244  IMTetrahedra()
245  :
246  q( 1 )
247  {
248  m += 4;
249  w += value_type( 1 )/value_type( 4 );
250  q[0] += 0.5854101966249685, 0.1381966011250105;
251  }
252 
253  std::vector<uint16_type> m;
254  std::vector<value_type> w;
255  std::vector<std::vector<value_type> > q;
256 };
257 
258 template<typename T>
259 struct IMTetrahedra<4,T>
260 {
261  typedef T value_type;
262  static const uint16_type nDim = 3;
263  static const uint16_type nOrder = 4;
264  static const uint16_type nPoints = 16;
265  IMTetrahedra()
266  :
267  q( 2 )
268  {
269  m += 4, 12;
270  w += 0.05037379410012282, 0.06654206863329239;
271  q[0] += 0.7716429020672371, 0.7611903264425430e-01;
272  q[1] += 0.1197005277978019, 0.7183164526766925e-01, 0.4042339134672644;
273  }
274 
275  std::vector<uint16_type> m;
276  std::vector<value_type> w;
277  std::vector<std::vector<value_type> > q;
278 };
279 template<typename T> struct IMTetrahedra<3,T>: public IMTetrahedra<4,T> {};
280 template<typename T>
281 struct IMTetrahedra<6,T>
282 {
283  typedef T value_type;
284  static const uint16_type nDim = 3;
285  static const uint16_type nOrder = 6;
286  static const uint16_type nPoints = 29;
287  IMTetrahedra()
288  :
289  q( 4 )
290  {
291  m += 1, 4, 12, 12;
292  w += 0.9040129046014750e-01, 0.1911983427899124e-01,
293  0.4361493840666568e-01, 0.2581167596199161e-01;
294 
295  q[0] += 0.25;
296  q[1] += 0.8277192480479295, 0.5742691731735683e-01;
297  q[2] += 0.5135188412556341e-01, 0.4860510285706072, 0.2312985436519147;
298  q[3] += 0.2967538129690260, 0.6081079894015281, 0.4756909881472290e-01;
299  }
300 
301  std::vector<uint16_type> m;
302  std::vector<value_type> w;
303  std::vector<std::vector<value_type> > q;
304 };
305 template<typename T> struct IMTetrahedra<5,T>: public IMTetrahedra<6,T> {};
306 } // detail
307 
308 
317 template<int Dim,int Order, typename T>
318 class IMSimplex
319  :
320 public PointSetQuadrature<Simplex<Dim,1> , Order, T>
321 {
322  typedef PointSetQuadrature<Simplex<Dim,1> , Order, T> super;
323 public:
324 
328  typedef T value_type;
329  typedef ublas::matrix<value_type,ublas::column_major> matrix_type;
330  typedef ublas::vector<value_type> vector_type;
331  typedef typename mpl::if_<mpl::equal_to<mpl::int_<Dim>,mpl::int_<2> >,
332  mpl::identity<detail::IMTriangle<Order,T> >,
333  mpl::identity<detail::IMTetrahedra<Order,T> > >::type::type quad_type;
334 
335 #if 1
336  typedef typename mpl::if_<mpl::equal_to<mpl::int_<Dim>,mpl::int_<2> >,
337  mpl::identity<Gauss<Simplex<Dim-1,1>,Order,T> >,
338  mpl::identity<IMSimplex<Dim-1,Order,T> > >::type::type face_quad_type;
339 #else
340  typedef IMSimplex<Dim-1,Order,T> face_quad_type;
341 #endif
342  static const uint16_type nDim = Dim;
343  static const uint16_type nOrder = Order;
344 
346 
350  IMSimplex( bool transform = true )
351  :
352  super( quad_type::nPoints ),
353  M_quad()
354  {
355  //std::cout << "npoints=" << quad_type::nPoints << "\n";
356  for ( size_type i=0, wi = 0; i< M_quad.q.size(); i++ )
357  {
358  permute( M_quad.m[i], M_quad.q[i], wi );
359 
360  for ( int j=0; j<M_quad.m[i]; j++, wi++ )
361  {
362  this->M_w( wi ) = factor()*M_quad.w[i];
363  }
364  }
365 
366 
367  if ( transform )
368  {
369  //std::cout << "order = " << Order << " | " << this->M_points << "\n";
370  typedef GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element_type;
371  RealToReference<element_type, Simplex,value_type> real_to_ref( element() );
372  this->setPoints( real_to_ref( this->M_points ) );
373  this->M_w /= real_to_ref.J();
374  //std::cout << "jac = " << real_to_ref.J() << "\n";
375  }
376 
377  //std::cout << "this->M_points=" << this->M_points << "\n"
378  //<< " w=" << this->M_w << "\n";
379 
380  boost::shared_ptr<GT_Lagrange<Dim,1,Simplex,T> > gm( new GT_Lagrange<Dim, 1, Simplex, T> );
381  boost::shared_ptr<face_quad_type> face_qr( new face_quad_type );
382  // construct face quadratures
383  this->constructQROnFace( Reference<Simplex<Dim, 1>, Dim, 1>(), gm, face_qr );
384  }
385 
386  ~IMSimplex()
387  {}
388 
390 
391  T factor() const
392  {
393  return ( Dim==2 )?T( 1 )/T( 2 ):T( 1 )/T( 6 );
394  }
395 
396 
400  GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element() const
401  {
402  return element( mpl::int_<nDim>() );
403  }
404 
405  GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element( mpl::int_<2> ) const
406  {
407  // setup the convex
408  typename node<value_type>::type pt( nDim );
409  typedef GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element_type;
410  typedef typename element_type::point_type point_type;
411  element_type elem;
412  // lower left
413  pt[0] = value_type( 0 );
414  pt[1] = value_type( 0 );
415  point_type p1( pt );
416  elem.setPoint( 0, p1 );
417 
418  // lower right
419  pt[0] = value_type( 1 );
420  pt[1] = value_type( 0 );
421  point_type p2( pt );
422  elem.setPoint( 1, p2 );
423 
424  // upper
425  pt[0] = value_type( 0 );
426  pt[1] = value_type( 1 );
427  point_type p3( pt );
428  elem.setPoint( 2, p3 );
429  return elem;
430 
431  }
432 
433  GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element( mpl::int_<3> ) const
434  {
435  // setup the convex
436  typename node<value_type>::type pt( nDim );
437  typedef GeoND<nDim,GeoEntity<Simplex<nDim,1> >, value_type > element_type;
438  typedef typename element_type::point_type point_type;
439  element_type elem;
440  // lower left
441  pt[0] = value_type( 0 );
442  pt[1] = value_type( 0 );
443  pt[2] = value_type( 0 );
444  point_type p1( pt );
445  elem.setPoint( 0, p1 );
446 
447  // lower right
448  pt[0] = value_type( 1 );
449  pt[1] = value_type( 0 );
450  pt[2] = value_type( 0 );
451  point_type p2( pt );
452  elem.setPoint( 1, p2 );
453 
454  //
455  pt[0] = value_type( 0 );
456  pt[1] = value_type( 1 );
457  pt[2] = value_type( 0 );
458  point_type p3( pt );
459  elem.setPoint( 2, p3 );
460 
461  //
462  pt[0] = value_type( 0 );
463  pt[1] = value_type( 0 );
464  pt[2] = value_type( 1 );
465  point_type p4( pt );
466  elem.setPoint( 3, p4 );
467  return elem;
468 
469  }
470 
471  template<typename QVec>
472  void
473  permute( int m, QVec const& q, int wi )
474  {
475  permute( m, q, wi, mpl::int_<Dim>() );
476  }
477  template<typename QVec>
478  void
479  permute( int m, QVec const& q, int wi, mpl::int_<2> )
480  {
481  switch ( m )
482  {
483  case 1:
484  {
485  this->M_points( 0, wi ) = q[ 0 ];
486  this->M_points( 1, wi ) = q[ 1 ];
487  wi++;
488  }
489  break;
490 
491  case 3:
492  {
493  this->M_points( 0, wi ) = q[ 0 ];
494  this->M_points( 1, wi ) = q[ 1 ];
495  wi++;
496 
497  this->M_points( 0, wi ) = q[ 1 ];
498  this->M_points( 1, wi ) = q[ 0 ];
499  wi++;
500 
501  this->M_points( 0, wi ) = q[ 2 ];
502  this->M_points( 1, wi ) = q[ 1 ];
503  wi++;
504 
505  }
506  break;
507 
508  case 6:
509  {
510  //qPerm[0] = tuple(q[0], q[1], q[2]);
511  this->M_points( 0, wi ) = q[ 0 ];
512  this->M_points( 1, wi ) = q[ 1 ];
513  wi++;
514 
515  //qPerm[1] = tuple(q[0], q[2], q[1]);
516  this->M_points( 0, wi ) = q[ 0 ];
517  this->M_points( 1, wi ) = q[ 2 ];
518  wi++;
519 
520  //qPerm[2] = tuple(q[1], q[0], q[2]);
521  this->M_points( 0, wi ) = q[ 1 ];
522  this->M_points( 1, wi ) = q[ 0 ];
523  wi++;
524 
525  //qPerm[3] = tuple(q[1], q[2], q[0]);
526  this->M_points( 0, wi ) = q[ 1 ];
527  this->M_points( 1, wi ) = q[ 2 ];
528  wi++;
529 
530  //qPerm[4] = tuple(q[2], q[1], q[0]);
531  this->M_points( 0, wi ) = q[ 2 ];
532  this->M_points( 1, wi ) = q[ 1 ];
533  wi++;
534 
535  //qPerm[5] = tuple(q[2], q[0], q[1]);
536  this->M_points( 0, wi ) = q[ 2 ];
537  this->M_points( 1, wi ) = q[ 0 ];
538  wi++;
539  }
540  break;
541  }
542  }
543  template<typename QVec>
544  void
545  permute( int m, QVec const& q, int wi, mpl::int_<3> )
546  {
547  if ( m==1 )
548  {
549  //qPerm[0] = boost::make_tuple(q[0], q[0], q[0], q[0]);
550  this->M_points( 0, wi ) = q[ 0 ];
551  this->M_points( 1, wi ) = q[ 0 ];
552  this->M_points( 2, wi ) = q[ 0 ];
553  }
554 
555  else if ( m==4 )
556  {
557  //qPerm[0] = boost::make_tuple(q[0], q[1], q[1], q[1]);
558  this->M_points( 0, wi ) = q[ 0 ];
559  this->M_points( 1, wi ) = q[ 1 ];
560  this->M_points( 2, wi ) = q[ 1 ];
561  wi++;
562 
563  //qPerm[1] = boost::make_tuple(q[1], q[0], q[1], q[1]);
564  this->M_points( 0, wi ) = q[ 1 ];
565  this->M_points( 1, wi ) = q[ 0 ];
566  this->M_points( 2, wi ) = q[ 1 ];
567  wi++;
568 
569  //qPerm[2] = boost::make_tuple(q[1], q[1], q[0], q[1]);
570  this->M_points( 0, wi ) = q[ 1 ];
571  this->M_points( 1, wi ) = q[ 1 ];
572  this->M_points( 2, wi ) = q[ 0 ];
573  wi++;
574 
575  //qPerm[3] = boost::make_tuple(q[1], q[1], q[1], q[0]);
576  this->M_points( 0, wi ) = q[ 1 ];
577  this->M_points( 1, wi ) = q[ 1 ];
578  this->M_points( 2, wi ) = q[ 1 ];
579  wi++;
580  }
581 
582  else if ( m==12 )
583  {
584  //qPerm[0] = boost::make_tuple(q[0], q[1], q[2], q[2]);
585  this->M_points( 0, wi ) = q[ 0 ];
586  this->M_points( 1, wi ) = q[ 1 ];
587  this->M_points( 2, wi ) = q[ 2 ];
588  wi++;
589 
590  //qPerm[1] = boost::make_tuple(q[0], q[2], q[1], q[2]);
591  this->M_points( 0, wi ) = q[ 0 ];
592  this->M_points( 1, wi ) = q[ 2 ];
593  this->M_points( 2, wi ) = q[ 1 ];
594  wi++;
595 
596  //qPerm[2] = boost::make_tuple(q[0], q[2], q[2], q[1]);
597  this->M_points( 0, wi ) = q[ 0 ];
598  this->M_points( 1, wi ) = q[ 2 ];
599  this->M_points( 2, wi ) = q[ 2 ];
600  wi++;
601 
602 
603  //qPerm[3] = boost::make_tuple(q[1], q[0], q[2], q[2]);
604  this->M_points( 0, wi ) = q[ 1 ];
605  this->M_points( 1, wi ) = q[ 0 ];
606  this->M_points( 2, wi ) = q[ 2 ];
607  wi++;
608 
609  //qPerm[4] = boost::make_tuple(q[2], q[0], q[1], q[2]);
610  this->M_points( 0, wi ) = q[ 2 ];
611  this->M_points( 1, wi ) = q[ 0 ];
612  this->M_points( 2, wi ) = q[ 1 ];
613  wi++;
614 
615  //qPerm[5] = boost::make_tuple(q[2], q[0], q[2], q[1]);
616  this->M_points( 0, wi ) = q[ 2 ];
617  this->M_points( 1, wi ) = q[ 0 ];
618  this->M_points( 2, wi ) = q[ 2 ];
619  wi++;
620 
621  //qPerm[6] = boost::make_tuple(q[1], q[2], q[0], q[2]);
622  this->M_points( 0, wi ) = q[ 1 ];
623  this->M_points( 1, wi ) = q[ 2 ];
624  this->M_points( 2, wi ) = q[ 0 ];
625  wi++;
626 
627  //qPerm[7] = boost::make_tuple(q[2], q[1], q[0], q[2]);
628  this->M_points( 0, wi ) = q[ 2 ];
629  this->M_points( 1, wi ) = q[ 1 ];
630  this->M_points( 2, wi ) = q[ 0 ];
631  wi++;
632 
633  //qPerm[8] = boost::make_tuple(q[2], q[2], q[0], q[1]);
634  this->M_points( 0, wi ) = q[ 2 ];
635  this->M_points( 1, wi ) = q[ 2 ];
636  this->M_points( 2, wi ) = q[ 0 ];
637  wi++;
638 
639  //qPerm[9] = boost::make_tuple(q[1], q[2], q[2], q[0]);
640  this->M_points( 0, wi ) = q[ 1 ];
641  this->M_points( 1, wi ) = q[ 2 ];
642  this->M_points( 2, wi ) = q[ 2 ];
643  wi++;
644 
645  //qPerm[10] = boost::make_tuple(q[2], q[1], q[2], q[0]);
646  this->M_points( 0, wi ) = q[ 2 ];
647  this->M_points( 1, wi ) = q[ 1 ];
648  this->M_points( 2, wi ) = q[ 2 ];
649  wi++;
650 
651  //qPerm[11] = boost::make_tuple(q[2], q[2], q[1], q[0]);
652  this->M_points( 0, wi ) = q[ 2 ];
653  this->M_points( 1, wi ) = q[ 2 ];
654  this->M_points( 2, wi ) = q[ 1 ];
655  wi++;
656  }
657 
658  else
659  {
660  FEELPP_ASSERT( 0 )( m ).error( "invalid multiplicity in IMSimplex<3>" );
661  }
662 
663  }
664 
665  bool test()
666  {
667  return test( mpl::int_<nDim>() );
668  }
669  bool test( mpl::int_<2> )
670  {
671  bool pass = true;
672  ublas::vector<value_type> x( ublas::row( this->M_points, 0 ) );
673  ublas::vector<value_type> y( ublas::row( this->M_points, 1 ) );
674  ublas::vector<value_type> w( this->M_w );
675 
676  for ( int a=0; a<=nOrder; a++ )
677  {
678  for ( int b=0; b<nOrder-a; b++ )
679  {
680  int cMax = nOrder - a - b;
681 
682  for ( int c=0; c<=cMax; c++ )
683  {
684  double sum = 0.0;
685 
686  for ( int q=0; q<w.size(); q++ )
687  {
688  sum += 0.5*w[q] * pow( x[q], ( double ) a ) * pow( y[q], ( double ) b )
689  * pow( 1.0 - x[q] - y[q], ( double ) c );
690  }
691 
692  double err = fabs( sum - exact( a,b,c ) );
693  bool localPass = err < 1.0e-14;
694  pass = pass && localPass;
695 
696  if ( !localPass )
697  {
698  fprintf( stderr, "order=%d m (%d, %d, %d) q=%22.15g exact=%22.15g\n", nOrder, a, b, c, sum, exact( a, b, c ) );
699  std::cerr << "error = " << err << std::endl;
700  }
701 
702  else
703  {
704  fprintf( stderr, "order=%d m (%d, %d, %d) q=%22.15g exact=%22.15g\n", nOrder, a, b, c, sum, exact( a, b, c ) );
705  }
706  }
707  }
708  }
709 
710  return pass;
711  }
712  bool test( mpl::int_<3> )
713  {
714  ublas::vector<value_type> x( ublas::row( this->M_points, 0 ) );
715  ublas::vector<value_type> y( ublas::row( this->M_points, 1 ) );
716  ublas::vector<value_type> z( ublas::row( this->M_points, 2 ) );
717  ublas::vector<value_type> w( this->M_w );
718 
719  bool pass = true;
720  int p = nOrder;
721 
722  for ( int a=0; a<=p; a++ )
723  {
724  for ( int b=0; b<p-a; b++ )
725  {
726  int cMax = p - a - b;
727 
728  for ( int c=0; c<=cMax; c++ )
729  {
730  int dMax = p - a - b - c;
731 
732  for ( int d=0; d<=dMax; d++ )
733  {
734  double sum = 0.0;
735 
736  for ( int q=0; q<w.size(); q++ )
737  {
738  sum += ( 1.0/6.0 )*w[q] * pow( x[q], ( double ) a ) * pow( y[q], ( double ) b )
739  * pow( z[q], ( double ) c )
740  * pow( 1.0 - x[q] - y[q] - z[q], ( double ) d );
741  }
742 
743  double err = fabs( sum - exact( a,b,c,d ) );
744  bool localPass = err < 1.0e-14;
745  pass = pass && localPass;
746 
747  if ( !localPass )
748  {
749  fprintf( stderr, "order=%d m (%d, %d, %d %d) q=%22.15g exact=%22.15g\n", p, a, b, c, d, sum, exact( a, b, c, d ) );
750  std::cerr << "error = " << err << std::endl;
751  }
752 
753  else
754  {
755  fprintf( stderr, "order=%d m (%d, %d, %d %d) q=%22.15g exact=%22.15g\n", p, a, b, c, d, sum, exact( a, b, c, d ) );
756  }
757  }
758  }
759  }
760  }
761 
762  return pass;
763 
764  }
765 
766  double exact( int a, int b, int c ) const
767  {
768  return fact( a )*fact( b )*fact( c )/fact( a+b+c+2 );
769  }
770  double exact( int a, int b, int c, int d ) const
771  {
772  return fact( a )*fact( b )*fact( c )*fact( d )/fact( a+b+c+d+3 );
773  }
774 
775  double fact( int x ) const
776  {
777  if ( x==0 ) return 1.0;
778 
779  return x*fact( x-1 );
780  }
782 
783 
784 private:
785  quad_type M_quad;
786 };
787 
788 #if 0
789 
790 /******* Particular Quadrature Rule from old Feel *******/
791 
794 template< typename T >
795 class TriangleQuadRule
796 {
797 public:
798 
799  typedef T value_type;
800  typedef typename node<value_type>::type node_type;
801  typedef ublas::matrix<value_type, ublas::column_major> nodes_type;
802  typedef ublas::vector<value_type> weights_type;
803 
804  nodes_type const& points() const
805  {
806  return _pts;
807  }
808 
809  ublas::matrix_column<nodes_type const> point( uint32_type __i ) const
810  {
811  return ublas::column( _pts, __i );
812  }
813 
814  weights_type const& weights() const
815  {
816  return M_w;
817  }
818  value_type const& weight( int q ) const
819  {
820  return M_w[q];
821  }
822 
823  value_type integrate( boost::function<value_type( node_type const& )> const& f ) const
824  {
825  value_type res = 0.0;
826 
827  for ( uint16_type k = 0; k < 7 ; ++k )
828  {
829  res += this->weights()[k]*f( this->point( k ) );
830  }
831 
832  return res;
833  }
834 
835  TriangleQuadRule()
836  : _pts( 2, 7 ), M_w( 7 )
837  {
838  value_type t7pt_x[3] = {value_type( 1.0 )/3.0, 0.10128650732345633, 0.47014206410511508};
839  value_type t7pt_w[3] = {0.1125, 0.062969590272413576, 0.066197076394253090};
840 
841  for ( int i = 0; i < 3 ; ++i )
842  {
843  t7pt_x[i]= 2.0*t7pt_x[i]-1.0;
844  t7pt_w[i]= 4.0*t7pt_w[i];
845  }
846 
847  M_w( 0 ) = t7pt_w[0];
848  M_w( 1 ) = t7pt_w[1];
849  M_w( 2 ) = t7pt_w[1];
850  M_w( 3 ) = t7pt_w[1];
851  M_w( 4 ) = t7pt_w[2];
852  M_w( 5 ) = t7pt_w[2];
853  M_w( 6 ) = t7pt_w[2];
854 
855  _pts( 0, 0 ) = t7pt_x[0];
856  _pts( 1, 0 ) = t7pt_x[0];
857 
858  _pts( 0, 1 ) = t7pt_x[1];
859  _pts( 1, 1 ) = t7pt_x[1];
860 
861  _pts( 0, 2 ) = t7pt_x[1];
862  _pts( 1, 2 ) = 1.0-2.0*t7pt_x[1];
863 
864  _pts( 0, 3 ) = 1.0-2.0*t7pt_x[1];
865  _pts( 1, 3 ) = t7pt_x[1];
866 
867  _pts( 0, 4 ) = t7pt_x[2];
868  _pts( 1, 4 ) = t7pt_x[2];
869 
870  _pts( 0, 5 ) = t7pt_x[2];
871  _pts( 1, 5 ) = 1.0-2.0*t7pt_x[2];
872 
873  _pts( 0, 6 ) = 1.0-2.0*t7pt_x[2];
874  _pts( 1, 6 ) = t7pt_x[2];
875  }
876 
877  ~TriangleQuadRule() {}
878 
879 protected:
880  nodes_type _pts;
881  weights_type M_w;
882 };
883 
884 
885 #endif // 0
886 }
887 
888 #endif /* __IMSimplex_H */

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