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
gausslobatto.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-03-06
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 __GaussLobatto_H
30 #define __GaussLobatto_H 1
31 
33 
35 
36 #include <stdexcept>
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>
41 
42 
43 #include <feel/feelcore/traits.hpp>
44 
45 #include <feel/feelpoly/jacobi.hpp>
46 #include <feel/feelpoly/quadpoint.hpp>
49 #include <feel/feelalg/glas.hpp>
50 
51 
52 
53 namespace Feel
54 {
55 namespace ublas = boost::numeric::ublas;
56 
57 
58 template<class Convex, uint16_type Integration_Degree, typename T> class PointSetQuadrature;
59 template<int Dim, int Order, int RealDim, template<uint16_type,uint16_type,uint16_type> class Entity, typename T> struct GT_Lagrange;
60 
75 template<class Convex, uint16_type Integration_Degree, typename T>
76 class GaussLobatto : public PointSetQuadrature<Convex, Integration_Degree, T> {};
77 
79 template< uint16_type Integration_Degree, typename T>
80 class GaussLobatto<Simplex<1,1> , Integration_Degree ,T > : public PointSetQuadrature<Simplex<1,1> , Integration_Degree, T>
81 {
82 public :
83  typedef T value_type;
84 
85  typedef PointSetQuadrature<Simplex<1,1> , Integration_Degree, T> super;
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;
90 
91  static const uint16_type Degree = ( Integration_Degree+3 )/2+1;
92  static const uint32_type Npoints = Degree;
93 
94  GaussLobatto()
95  :
96  super( Npoints )
97  {
98  ublas::vector<T> px( Npoints );
99 
100  details::gausslobattojacobi<Npoints, T, ublas::vector<T>, ublas::vector<T> >( this->M_w, px );
101  ublas::row( this->M_points, 0 ) = px;
102  }
103 
104  ~GaussLobatto() {}
105 
106  FEELPP_DEFINE_VISITABLE();
107 };
108 
109 
112 template< uint16_type Integration_Degree, typename T>
113 class GaussLobatto<Simplex<2,1> , Integration_Degree ,T > : public PointSetQuadrature<Simplex<2,1> , Integration_Degree, T>
114 {
115 public :
116  typedef T value_type;
117 
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;
123 
124  typedef GaussLobatto<Simplex<1,1>,Integration_Degree, T> face_quad_type;
125 
126 
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;
130 
131  GaussLobatto()
132  :
133  super( Npoints )
134  {
135  // build rules in x and y direction
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 );
139 
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 );
143 
144  // coordinate in cartesian space
145  node_type eta( 2 );
146  details::xi<TRIANGLE, value_type> to_xi;
147 
148  for ( int i = 0, k = 0; i < DegreeX; ++i )
149  {
150  for ( int j = 0; j < DegreeY; ++j, ++k )
151  {
152  // computes the weight of the k-th node
153  this->M_w( k ) = 0.5 * wx( i ) * wy( j );
154  // use expansion for the collapsed triangle to compute the points
155  // coordinates (from cartesian to collapsed coordinates)
156  eta( 0 ) = px( i );
157  eta( 1 ) = py( j );
158  ublas::column( this->M_points, k ) = to_xi( eta );
159  }
160  }
161 
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 );
164  // construct face quadratures
165  this->constructQROnFace( Reference<Simplex<2, 1, 2>,2,1>(), gm, face_qr );
166 
167  }
168 
169  ~GaussLobatto() {}
170 
171  FEELPP_DEFINE_VISITABLE();
172 };
173 
176 template< uint16_type Integration_Degree, typename T>
177 class GaussLobatto<Simplex<3,1> , Integration_Degree ,T > : public PointSetQuadrature<Simplex<3,1> , Integration_Degree, T>
178 {
179 public :
180  typedef T value_type;
181 
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;
187 
188  typedef GaussLobatto<Simplex<2,1>,Integration_Degree, T> face_quad_type;
189 
190 
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;
195 
196  GaussLobatto()
197  :
198  super( Npoints )
199  {
200  // build rules in x and y direction
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 );
204 
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 );
208 
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 );
212 
213  // coordinate in cartesian space
214  node_type eta( 3 );
215  details::xi<TETRAHEDRON, value_type> to_xi;
216 
217  for ( int i = 0, k = 0; i < DegreeX; ++i )
218  {
219  for ( int j = 0; j < DegreeY; ++j )
220  {
221  for ( int l = 0; l < DegreeZ; ++l, ++k )
222  {
223  // computes the weight of the k-th node
224  this->M_w( k ) = 0.125 * wx( i ) * wy( j ) * wz( l );
225  // use expansion for the collapsed triangle to compute the points
226  // coordinates (from cartesian to collapsed coordinates)
227  eta( 0 ) = px( i );
228  eta( 1 ) = py( j );
229  eta( 2 ) = pz( l );
230  ublas::column( this->M_points, k ) = to_xi( eta );
231  }
232  }
233  }
234 
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 );
237  // construct face quadratures
238  this->constructQROnFace( Reference<Simplex<3, 1, 3>,3,1>(), gm, face_qr );
239  }
240 
241  ~GaussLobatto() {}
242 
243  FEELPP_DEFINE_VISITABLE();
244 };
245 
246 
247 
248 
253 template< uint16_type Integration_Degree, typename T>
254 class GaussLobatto<Hypercube<2,1>, Integration_Degree ,T >
255  :
256 public PointSetQuadrature<Hypercube<2,1>, Integration_Degree, T>
257 {
258 public :
259  typedef T value_type;
260 
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;
269 
270  GaussLobatto()
271  :
272  super( Npoints )
273  {
274  // build rules in x and y direction
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 );
278 
279  for ( int i = 0, k = 0; i < Degree; ++i )
280  {
281  for ( int j = 0; j < Degree; ++j, ++k )
282  {
283  // computes the weight of the k-th node
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 );
287  }
288  }
289 
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 );
292  // construct face quadratures
293  this->constructQROnFace( Reference<Hypercube<2, 1, 2>,2,1>(), gm, face_qr );
294  }
295 
296  ~GaussLobatto() {}
297 
298  FEELPP_DEFINE_VISITABLE();
299 };
300 
303 template< uint16_type Integration_Degree, typename T>
304 class GaussLobatto<Hypercube<3,1>, Integration_Degree ,T >
305  :
306 public PointSetQuadrature<Hypercube<3,1>, Integration_Degree, T>
307 {
308 public :
309  typedef T value_type;
310 
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;
319 
320  GaussLobatto()
321  :
322  super( Npoints )
323  {
324  // build rules in x and y direction
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 );
328 
329  for ( int i = 0, k = 0; i < Degree; ++i )
330  {
331  for ( int j = 0; j < Degree; ++j )
332  {
333  for ( int l = 0; l < Degree ; ++l, ++k )
334  {
335  // computes the weight of the k-th node
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 );
340  }
341  }
342  }
343 
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 );
346  // construct face quadratures
347  this->constructQROnFace( Reference<Hypercube<3, 1, 3>,3,1>(), gm, face_qr );
348  }
349 
350  ~GaussLobatto() {}
351  FEELPP_DEFINE_VISITABLE();
352 };
354 
355 
364 template< class Convex,
365  uint16_type Order,
366  typename T = double >
367 class PointSetGaussLobatto : public PointSetInterpolation<Convex::nDim, Order, T, Hypercube>
368 {
369 public:
370 
371  typedef PointSetInterpolation< Convex::nDim, Order, T, Hypercube> super;
372 
373  typedef typename super::return_type return_type;
374 
375  typedef T value_type;
376 
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;
381 
382  static const size_type Shape = Convex::Shape;
383 
384  static const bool is_simplex = Convex::is_simplex;
385  static const bool is_hypercube = Convex::is_hypercube;
386 
387  typedef Reference<Convex, Dim, Convex::nOrder, Convex::nDim/*Convex::nRealDim*/, value_type> reference_convex_type;
388 
389  typedef ublas::vector<value_type> vector_type;
390 
391  typedef typename super::node_type node_type;
392  typedef typename super::nodes_type nodes_type;
393 
394  typedef typename reference_convex_type::points_type points_type;
395 
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;
399 
400  typedef Hypercube<Dim, Order, Dim> conv_order_type;
401  static const uint32_type numPoints = conv_order_type::numPoints;
402 
403  reference_convex_type RefConv;
404 
405  typedef typename super::range_type range_type;
406  typedef typename super::index_map_type index_map_type;
407 
408  PointSetGaussLobatto( int interior = 0 )
409  {
410  FEELPP_ASSERT( is_hypercube || ( Dim == 1 ) ).error( "gauss lobatto points are just defined in simplex products" );
411 
412  nodes_type pts( Dim, numPoints );
413 
414  calculate_gl_points( mpl::bool_< ( Order > 1 )>() );
415 
416  if ( interior == 0 && Order > 0 )
417  {
418  for ( uint16_type d = 0, p = 0; d < topological_dimension+1; ++d )
419  {
420  for ( int e = RefConv.entityRange( d ).begin();
421  e < RefConv.entityRange( d ).end();
422  ++e )
423  {
424  nodes_type Gt ( makePoints( d, e ) );
425 
426  if ( Gt.size2() )
427  {
428  ublas::subrange( pts, 0, Dim, p, p+Gt.size2() ) = Gt;
429 
430  for ( size_type j = 0; j < Gt.size2(); ++j )
431  {
432  this->addToEid( d, p+j );
433  this->addToPtE( p+j, std::make_pair( d, e ) );
434  }
435 
436  p+=Gt.size2();
437  }
438  }
439  }
440 
441  this->setPoints( pts );
442  }
443 
444  else if ( interior == 1 && Order > 0 )
445  this->setPoints( makePoints( Dim, 0 ) );
446 
447  else if ( Order == 0 )
448  this->setPoints( glas::average( RefConv.vertices() ) );
449 
450  this->setName( "gausslobatto", Order );
451 
452  //std::cout << "end constructor...\n";
453  }
454 
455  ~PointSetGaussLobatto() {}
456 
457 private:
458 
459  vector_type gl_pts;
460 
461  void calculate_gl_points ( mpl::bool_<true> )
462  {
463  vector_type wx( Order+1 );
464  vector_type px( Order+1 );
465 
466  details::gausslobattojacobi<Order+1, T, vector_type, vector_type >( wx, px );
467 
468  ublas::vector_range<vector_type> inner_pts ( px, ublas::range( 1, px.size()-1 ) );
469 
470  gl_pts = inner_pts;
471  }
472 
473  void calculate_gl_points ( mpl::bool_<false> )
474  {}
475 
476  points_type makePoints( uint16_type topo_dim, uint16_type __id )
477  {
478  // vertices
479  if ( topo_dim == 0 )
480  {
481  points_type G( RefConv.vertices().size1(), 1 );
482  ublas::column( G, 0 ) = ublas::column( RefConv.vertices(), __id );
483  return G;
484  }
485 
486  // interior points of the convex
487  else if ( topo_dim == topological_dimension )
488  {
489  if ( __id == 0 )
490  return makeLattice<Shape>();
491 
492  throw std::logic_error( "cannot make those points" );
493  return points_type();
494  }
495 
496  // all the other points
497  else
498  {
499  points_type G;
500  points_type Gret;
501 
502  if ( topo_dim == 1 )
503  {
504  G = makeLattice<SHAPE_LINE>();
505  Gret.resize( nRealDim, G.size2() );
506 
507  pt_to_entity_hexahedron<Shape, 1> p_to_e( __id );
508 
509  for ( size_type i = 0; i < G.size2(); ++i )
510  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
511 
512  return Gret;
513  }
514 
515  else if ( topo_dim == 2 )
516  {
517  G = makeLattice<SHAPE_QUAD>();
518  Gret.resize( nRealDim, G.size2() );
519  pt_to_entity_hexahedron<Shape, 2> p_to_e( __id );
520 
521  for ( size_type i = 0; i < G.size2(); ++i )
522  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
523 
524  return Gret;
525  }
526  }
527 
528  return points_type();
529  }
530 
531  template<size_type shape>
532  points_type makeLattice()
533  {
534  points_type G;
535 
536  if ( Order > 0 )
537  {
538  if ( shape == SHAPE_LINE )
539  G = make_line_points();
540 
541  else if ( shape == SHAPE_QUAD )
542  return make_quad_points();
543 
544  else if ( shape == SHAPE_HEXA )
545  return make_hexa_points();
546  }
547 
548  else if ( Order == 0 )
549  G = glas::average( RefConv.vertices() );
550 
551  return G;
552  }
553 
554  int n_line_points()
555  {
556  return std::max( 0, int( Order ) - 1 );
557  }
558 
559  int n_quad_points() const
560  {
561  return std::max( 0, ( int( Order ) - 1 )*( int( Order ) - 1 ) );
562 
563  }
564 
565  int n_hexa_points() const
566  {
567  return std::max( 0, ( int( Order )-1 )*( int( Order )-1 )*( int( Order )-1 ) );
568  }
569 
570  points_type
571  make_line_points()
572  {
573  points_type G ( Dim, n_line_points() );
574 
575  if ( Order > 1 )
576  {
577  vector_type ones ( ublas::scalar_vector<value_type>( G.size2(), value_type( 1 ) ) );
578 
579  ublas::row( G, 0 ) = gl_pts;
580 
581  for ( uint16_type i=1; i<Dim; i++ )
582  ublas::row( G, i ) = - ones;
583  }
584 
585  return G;
586  }
587 
588  points_type
589  make_quad_points()
590  {
591  points_type G( Dim, n_quad_points() );
592 
593  if ( Order > 1 )
594  {
595  uint16_type numInterior = Order - 1;
596 
597  for ( uint16_type i = 0, k = 0; i < numInterior; ++i )
598  {
599  for ( uint16_type j = 0; j < numInterior; ++j, ++k )
600  {
601  G( 0, k ) = gl_pts( i );
602  G( 1, k ) = gl_pts( j );
603  }
604  }
605 
606  vector_type ones ( ublas::scalar_vector<value_type>( G.size2(), value_type( 1 ) ) );
607 
608  if ( Dim == 3 )
609  ublas::row( G, 2 ) = - ones;
610  }
611 
612  return G;
613  }
614 
615  points_type
616  make_hexa_points()
617  {
618  points_type G( Dim, n_hexa_points() );
619 
620  if ( Order > 1 )
621  {
622  uint16_type numInterior = Order - 1;
623 
624  for ( uint16_type i = 0, k = 0; i < numInterior; ++i )
625  {
626  for ( uint16_type j = 0; j < numInterior; ++j )
627  {
628  for ( uint16_type l = 0; l < numInterior; ++l, ++k )
629  {
630  G( 0, k ) = gl_pts( i );
631  G( 1, k ) = gl_pts( j );
632  G( 2, k ) = gl_pts( l );
633  }
634  }
635  }
636  }
637 
638  return G;
639  }
640 
641  template<size_type shape>
642  struct pt_to_edge
643  {
644  pt_to_edge( std::vector<uint16_type> vert_ids )
645  :
646  h( 1.0 ),
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 ] ) ),
651  diff( v-u )
652  {
653  h = 1.0/( b[0]-a[0] );
654  }
655  node_type
656  operator()( node_type const& x ) const
657  {
658  return u + h * ( x[ 0 ] - a[ 0 ] ) * diff;
659  }
660  value_type h;
661  node_type a, b;
662  node_type u, v, diff;
663  };
664 
665 
666  //
667  // pt_to_face hexa
668  //
669  template<size_type shape>
670  struct pt_to_face_hexahedron
671  {
672  pt_to_face_hexahedron( std::vector<uint16_type> vert_ids )
673  :
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 ] ) ),
677  diff( 2 )
678  {
679  diff[0] = v-u;
680  diff[1] = w-u;
681  }
682  node_type
683  operator()( node_type const& x ) const
684  {
685  return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
686  }
687  node_type u, v, w;
688  ublas::vector<node_type> diff;
689  };
690 
691  template<size_type shape>
692  struct pt_to_element
693  {
694  pt_to_element() {}
695  pt_to_element( std::vector<uint16_type> const& ) {}
696  node_type operator()( node_type const& x ) const
697  {
698  return x;
699  }
700  };
701 
702  template<size_type shape,uint16_type topo_dim>
703  struct pt_to_entity_hexahedron
704  {
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> > >
710  >::type // 2
711  >::type::type _type;
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;
714 
715  pt_to_entity_hexahedron( uint16_type entity_id )
716  :
717  mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
718  {}
719 
720  node_type operator()( node_type const& x ) const
721  {
722  return mapping( x );
723  }
724  mapping_type mapping;
725  };
726 
727 };
728 } // Feel
729 #endif /* __GaussLobatto_H */

Generated on Fri Oct 25 2013 14:24:13 for Feel++ by doxygen 1.8.4