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
equispaced.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 */
30 #ifndef __PointSetEquiSpaced_H
31 #define __PointSetEquiSpaced_H 1
32 
37 
39 #include <feel/feelcore/traits.hpp>
40 
41 #include <stdexcept>
42 #include <boost/numeric/ublas/matrix.hpp>
43 #include <boost/numeric/ublas/io.hpp>
44 
45 #include <feel/feelalg/glas.hpp>
46 
47 
48 namespace Feel
49 {
50 namespace ublas = boost::numeric::ublas;
51 
52 template<class Convex, uint16_type Order, typename T>
53 class PointSetEquiSpaced : public PointSet< Convex, T>
54 {
55 
56 public :
57 
58  typedef PointSet<Convex, T> super;
59 
60  typedef T value_type;
61 
62  typedef typename super::return_type return_type;
63  typedef typename super::self_type self_type;
64 
65  typedef typename super::node_type node_type;
66 
67  typedef typename super::nodes_type nodes_type;
68 
69  typedef typename matrix_node<value_type>::type points_type;
70 
71  static const uint32_type Dim = Convex::nDim;
72  static const uint32_type convexOrder = Convex::nOrder;
73  static const uint32_type topological_dimension = Convex::topological_dimension;
74  static const uint32_type nRealDim = Convex::nRealDim;
75 
76  static const size_type Shape = Convex::Shape;
77 
78  static const bool is_simplex = Convex::is_simplex;
79  static const bool is_hypercube = Convex::is_hypercube;
80 
81  typedef mpl::if_< mpl::bool_< is_simplex >,
82  Simplex<Dim, Order, /*nRealDim*/Dim> ,
83  Hypercube<Dim, Order, /*nRealDim*/Dim> > conv_order_type;
84 
85  typedef Reference<Convex, Dim, convexOrder, Dim/*nRealDim*/, value_type> RefElem;
86 
87  static const uint32_type numPoints = conv_order_type::type::numPoints;
88  static const uint32_type nbPtsPerVertex = conv_order_type::type::nbPtsPerVertex;
89  static const uint32_type nbPtsPerEdge = conv_order_type::type::nbPtsPerEdge;
90  static const uint32_type nbPtsPerFace = conv_order_type::type::nbPtsPerFace;
91  static const uint32_type nbPtsPerVolume = conv_order_type::type::nbPtsPerVolume;
92 
93  typedef typename Convex::edge_to_point_t edge_to_point_t;
94  typedef typename Convex::face_to_point_t face_to_point_t;
95  typedef typename Convex::face_to_edge_t face_to_edge_t;
96 
97  typedef std::pair<uint16_type, uint16_type> range_type;
98  typedef std::vector< std::vector<size_type> > index_map_type;
99 
100  RefElem RefConv;
101 
102  PointSetEquiSpaced( int interior = 0 )
103  :
104  super( numPoints, Dim ),
105  M_eid()
106  {
107  M_eid.resize( topological_dimension + 1 );
108  M_pt_to_entity.resize( numPoints );
109 
110  nodes_type pts( Dim, numPoints );
111 
112  if ( interior == 0 && Order > 0 )
113  {
114  // loop on each convex of topological dimension <= to the current convex
115  // where we build the polynomial set
116  for ( uint16_type d = 0, p = 0; d < topological_dimension+1; ++d )
117  {
118  // loop on each entity forming the convex of topological
119  // dimension d
120  for ( int e = RefConv.entityRange( d ).begin();
121  e < RefConv.entityRange( d ).end();
122  ++e )
123  {
124  nodes_type Gt ( makePoints( d, e ) );
125 
126  if ( Gt.size2() )
127  {
128  ublas::subrange( pts, 0, Dim, p, p+Gt.size2() ) = Gt;
129 
130  for ( size_type j = 0; j < Gt.size2(); ++j )
131  {
132  addToEid( d, p+j );
133  addToPtE( p+j, std::make_pair( d, e ) );
134  }
135 
136  p+=Gt.size2();
137  }
138  }
139  }
140 
141  this->setPoints( pts );
142  }
143 
144  else if ( interior == 1 && Order > 0 )
145  this->setPoints( makePoints( Dim, 0 ) );
146 
147  else if ( Order == 0 )
148  this->setPoints( glas::average( RefConv.vertices() ) );
149 
150  this->setName( "equispaced", Order );
151  }
152 
153  ~PointSetEquiSpaced() {}
154 
155  ublas::matrix_range<nodes_type const> pointsByEntity( uint16_type e ) const
156  {
157  FEELPP_ASSERT( M_eid[e].size() )( e ).error( "no points defined on this entity" );
158 
159  return ublas::project( this->points(),
160  ublas::range( 0,Dim ),
161  ublas::range( *M_eid[e].begin(), *M_eid[e].rbegin() + 1 ) );
162  }
163 
164  std::pair<uint16_type, uint16_type> interiorRangeById( uint16_type e, uint16_type id ) const
165  {
166  uint16_type numEntities = 1;
167 
168  if ( e == 0 )
169  numEntities = Convex::numVertices;
170 
171  else if ( e == 1 )
172  numEntities = Convex::numEdges;
173 
174  else if ( e == 2 )
175  numEntities = Convex::numFaces;
176 
177  uint16_type N = M_eid[e].size()/numEntities;
178 
179  return std::make_pair( *M_eid[e].begin() + id*N, *M_eid[e].begin() + ( id+1 )*N );
180  }
181 
182  ublas::matrix_range<nodes_type const> interiorPointsById( uint16_type e, uint16_type id ) const
183  {
184  std::pair<uint16_type, uint16_type> position = interiorRangeById( e, id );
185 
186  ublas::matrix_range<nodes_type const> G = ublas::project( this->points(),
187  ublas::range( 0, Dim ),
188  ublas::range( position.first, position.second ) );
189 
190  return G;
191  }
192 
193  uint32_type entityIds( int i, int j ) const
194  {
195  return M_eid[i][j];
196  }
197 
198  uint32_type numEntities( int i ) const
199  {
200  return M_eid[i].size();
201  }
202 
203  std::pair<uint16_type, uint16_type> const& pointToEntity( int p ) const
204  {
205  return M_pt_to_entity[p];
206  }
207 
208  //Returns the local indices of all the subentities that compose the entity
209  index_map_type entityToLocal ( uint16_type top_dim, uint16_type local_id, bool boundary = 0 ) const
210  {
211  index_map_type indices( top_dim+1 );
212 
213  if ( top_dim == 0 && boundary )
214  {
215  range_type pair = interiorRangeById( top_dim, local_id );
216 
217  indices[0].push_back( pair.first );
218  }
219 
220  else
221  {
222  if ( !boundary && top_dim > 0 )
223  {
224  if ( numEntities( top_dim ) != 0 )
225  indices[top_dim].push_back( local_id );
226  }
227 
228  else
229  {
230  if ( top_dim > 0 )
231  {
232  //For the time being, this definition works, but in a more general framework,
233  //a entity shoulkd be created here with the respective top_dim
234  //in order to retrieve the number of vertices, edges, faces.
235 
236  //number of vertices, number of edges and number of faces in the volume
237  std::map<uint16_type, uint16_type> numPointsInEntity;
238 
239  numPointsInEntity[0] = ( uint16_type ) ( is_simplex )?( top_dim+1 ):( 2 + ( top_dim-1 )*( 2+3*( top_dim-2 ) ) );
240  numPointsInEntity[1] = ( uint16_type ) ( 2 + ( top_dim-1 )*( 1+2*is_simplex + ( 5-4*is_simplex )*( top_dim-1 ) ) )/2;
241  numPointsInEntity[2] = ( uint16_type ) 6-2*is_simplex;
242 
243  for ( uint16_type i=0; i < numPointsInEntity[0] ; i++ )
244  {
245  if ( top_dim == 1 )
246  indices[0].push_back( RefConv.e2p( local_id, i ) );
247 
248  else if ( top_dim == 2 )
249  indices[0].push_back( RefConv.f2p( local_id, i ) );
250  }
251 
252  if ( ( top_dim == 2 ) && ( M_eid[1].size() != 0 ) )
253  {
254  for ( uint16_type i=0; i < numPointsInEntity[1] ; i++ )
255  indices[1].push_back( RefConv.f2e( local_id, i ) );
256  }
257 
258  if ( top_dim == 3 )
259  {
260  for ( uint16_type k=0; k<numPointsInEntity.size(); k++ )
261  for ( uint16_type i=0; i < numPointsInEntity[k]; i++ )
262  indices[k].push_back( i );
263  }
264 
265  if ( M_eid[top_dim].size() )
266  indices[top_dim].push_back( local_id );
267  }
268  }
269  }
270 
271  return indices;
272  }
273 
274 
275  points_type pointsBySubEntity( uint16_type top_dim, uint16_type local_id, bool boundary = 0 ) const
276  {
277  index_map_type index_list = entityToLocal( top_dim, local_id, boundary );
278 
279  uint16_type matrix_size = 0;
280 
281  if ( index_list[0].size() != 0 )
282  matrix_size +=index_list[0].size()*nbPtsPerVertex;
283 
284  if ( ( top_dim >= 1 ) && ( index_list[1].size() != 0 ) )
285  matrix_size +=index_list[1].size()*nbPtsPerEdge;
286 
287  if ( ( top_dim >= 2 ) && ( index_list[2].size() != 0 ) )
288  matrix_size +=index_list[2].size()*nbPtsPerFace;
289 
290  if ( ( top_dim == 3 ) && ( index_list[3].size() != 0 ) )
291  matrix_size +=nbPtsPerVolume;
292 
293  points_type G ( Dim, matrix_size );
294 
295  for ( uint16_type i=0, p=0; i < top_dim+1; i++ )
296  {
297  if ( index_list[i].size() )
298  {
299  for ( uint16_type j=0; j < index_list[i].size(); j++ )
300  {
301  points_type aux = interiorPointsById( i, index_list[i][j] );
302 
303  ublas::subrange( G, 0, Dim, p, p+aux.size2() ) = aux;
304 
305  p+=aux.size2();
306  }
307  }
308  }
309 
310  return G;
311  }
312 
313  index_map_type getEid ()
314  {
315  return M_eid;
316  }
317 
318  std::vector<range_type> getPtE()
319  {
320  return M_pt_to_entity;
321  }
322 
323  void setEid ( index_map_type eid )
324  {
325  M_eid = eid;
326  }
327 
328  void setPtE ( std::vector<range_type> pt_ent )
329  {
330  M_pt_to_entity = pt_ent;
331  }
332 
333  void addToEid ( uint16_type p, uint16_type q )
334  {
335  M_eid[p].push_back( q );
336  }
337 
338  void addToPtE ( uint16_type p, range_type q )
339  {
340  M_pt_to_entity[p] = q;
341  }
342 
344 
345 private:
346 
347  index_map_type M_eid;
348  std::vector<range_type> M_pt_to_entity;
349 
350  points_type makePoints( uint16_type topo_dim, uint16_type __id, int interior = 1 )
351  {
352  // vertices
353  if ( topo_dim == 0 )
354  {
355  points_type G( RefConv.vertices().size1(), 1 );
356  ublas::column( G, 0 ) = ublas::column( RefConv.vertices(), __id );
357  return G;
358  }
359 
360  // interior points of the convex
361  else if ( topo_dim == topological_dimension )
362  {
363  if ( __id == 0 )
364  return makeLattice<Shape>( interior );
365 
366  throw std::logic_error( "cannot make those points" );
367  return points_type();
368  }
369 
370  // all the other points
371  else
372  {
373  points_type G;
374  points_type Gret;
375 
376  if ( topo_dim == 1 )
377  {
378  G = makeLattice<SHAPE_LINE>( 1 );
379  Gret.resize( nRealDim, G.size2() );
380 
381  if ( is_simplex )
382  {
383  pt_to_entity_tetrahedron<Shape, 1> p_to_e( __id );
384 
385  for ( size_type i = 0; i < G.size2(); ++i )
386  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
387  }
388 
389  else
390  {
391  pt_to_entity_hexahedron<Shape, 1> p_to_e( __id );
392 
393  for ( size_type i = 0; i < G.size2(); ++i )
394  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
395  }
396 
397  return Gret;
398  }
399 
400  else if ( topo_dim == 2 )
401  {
402  if ( is_simplex )
403  {
404  G = makeLattice<SHAPE_TRIANGLE>( 1 );
405  Gret.resize( nRealDim, G.size2() );
406  pt_to_entity_tetrahedron<Shape, 2> p_to_e( __id );
407 
408  for ( size_type i = 0; i < G.size2(); ++i )
409  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
410  }
411 
412  else
413  {
414  G = makeLattice<SHAPE_QUAD>( 1 );
415  Gret.resize( nRealDim, G.size2() );
416  pt_to_entity_hexahedron<Shape, 2> p_to_e( __id );
417 
418  for ( size_type i = 0; i < G.size2(); ++i )
419  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
420  }
421 
422  return Gret;
423  }
424  }
425 
426  return points_type();
427  }
428 
429  template<size_type shape>
430  points_type makeLattice( uint16_type interior = 0 )
431  {
432  points_type G;
433 
434  if ( Order > 0 )
435  {
436  if ( shape == SHAPE_LINE )
437  G = make_line_points( interior );
438 
439  else if ( shape == SHAPE_TRIANGLE )
440  G = make_triangle_points( interior );
441 
442  else if ( shape == SHAPE_TETRA )
443  G = make_tetrahedron_points( interior );
444 
445  else if ( shape == SHAPE_QUAD )
446  return make_quad_points( interior );
447 
448  else if ( shape == SHAPE_HEXA )
449  return make_hexa_points( interior );
450  }
451 
452  else if ( Order == 0 )
453  G = glas::average( RefConv.vertices() );
454 
455  return G;
456  }
457 
458  //---------------------------------------------------------------------------------------------
459  int n_line_points( int interior = 0 )
460  {
461  return std::max( 0, int( Order )+1-2*interior );
462  }
463  int n_triangle_points( int interior = 0 )
464  {
465  if ( interior == 1 )
466  return std::max( 0, ( int( Order )+1-2*interior )*( int( Order )-2*interior )/2 );
467 
468  return ( Order+1 )*( Order+2 )/2;
469  }
470  int n_tetrahedron_points( int interior = 0 )
471  {
472  if ( interior == 1 )
473  return std::max( 0, ( int( Order )+1-2*interior )*( int( Order )-2*interior )*( int( Order )-1-2*interior )/6 );
474 
475  return ( Order+1 )*( Order+2 )*( Order+3 )/6;
476  }
477 
478  int n_quad_points( int interior = 0 ) const
479  {
480  if ( interior == 1 )
481  return std::max( 0, ( int( Order )+1-2*interior )*( int( Order )+1-2*interior ) );
482 
483  return ( Order+1 )*( Order+1 );
484  }
485 
486  int n_hexa_points( int interior = 0 ) const
487  {
488  if ( interior == 1 )
489  return std::max( 0, ( int( Order )+1-2*interior )*( int( Order )+1-2*interior )*( int( Order )+1-2*interior ) );
490 
491  return ( Order+1 )*( Order+1 )*( Order+1 );
492  }
493 
494  points_type
495  make_line_points( int interior = 0 )
496  {
497  points_type G;
498 
499  if ( Order > 0 )
500  {
501  ublas::vector<node_type> h ( 1 );
502  h( 0 ) = RefConv.vertex( 1 ) - RefConv.vertex( 0 );
503 
504  G.resize( Dim, n_line_points( interior ) );
505 
506  for ( int i = interior, indp = 0; i < int( Order )+1-interior; ++i, ++indp )
507  {
508  ublas::column( G, indp ) = RefConv.vertex( 0 ) + ( h( 0 ) * value_type( i ) )/value_type( Order );
509  }
510  }
511 
512  else
513  G = glas::average( RefConv.vertices() );
514 
515  return G;
516 
517  }
518 
519  points_type
520  make_triangle_points( int interior = 0 )
521  {
522  points_type G;
523 
524  if ( Order > 0 )
525  {
526  ublas::vector<node_type> h ( 2 );
527  h( 0 ) = RefConv.vertex( 1 ) - RefConv.vertex( 0 );
528  h( 1 ) = RefConv.vertex( 2 ) - RefConv.vertex( 0 );
529 
530  G.resize( Dim, n_triangle_points( interior ) );
531 
532  for ( int i = interior, p = 0; i < int( Order )+1-interior; ++i )
533  {
534  for ( int j = interior; j < int( Order ) + 1 - i-interior; ++j, ++p )
535  {
536  ublas::column( G, p ) = RefConv.vertex( 0 ) + ( value_type( i ) * h( 1 ) +
537  value_type( j ) * h( 0 ) )/ value_type( Order );
538  }
539  }
540  }
541 
542  else
543  G = glas::average( RefConv.vertices() );
544 
545  return G;
546  }
547 
548  points_type
549  make_tetrahedron_points( int interior = 0 )
550  {
551  points_type G;
552 
553  if ( Order > 0 )
554  {
555  ublas::vector<node_type> h ( 3 );
556  h( 0 ) = RefConv.vertex( 1 ) - RefConv.vertex( 0 );
557  h( 1 ) = RefConv.vertex( 2 ) - RefConv.vertex( 0 );
558  h( 2 ) = RefConv.vertex( 3 ) - RefConv.vertex( 0 );
559 
560  G.resize( Dim, n_tetrahedron_points( interior ) );
561 
562  for ( int i = interior, p = 0; i < int( Order )+1-interior; ++i )
563  {
564  for ( int j = interior; j < int( Order ) + 1 - i - interior; ++j )
565  {
566  for ( int k = interior; k < int( Order ) + 1 - i - j - interior; ++k, ++p )
567  {
568  ublas::column( G, p ) = RefConv.vertex( 0 ) + ( value_type( i ) * h( 2 ) +
569  value_type( j ) * h( 1 ) +
570  value_type( k ) * h( 0 ) ) / value_type( Order );
571 
572  }
573  }
574  }
575  }
576 
577  else
578  G = glas::average( RefConv.vertices() );
579 
580  return G;
581  }
582 
583  points_type
584  make_quad_points( int interior = 0 )
585  {
586  if ( Order > 0 )
587  {
588  ublas::vector<node_type> h ( 2 );
589  h( 0 ) = RefConv.vertex( 1 ) - RefConv.vertex( 0 );
590  h( 1 ) = RefConv.vertex( 3 ) - RefConv.vertex( 0 );
591 
592  DVLOG(2) << "n quad pts = " << n_quad_points( interior ) << "\n";
593  points_type G( Dim, n_quad_points( interior ) );
594 
595  for ( int i = interior, p = 0; i < int( Order )+1-interior; ++i )
596  {
597  for ( int j = interior; j < int( Order ) + 1 -interior; ++j, ++p )
598  {
599  ublas::column( G, p ) = RefConv.vertex( 0 ) + ( value_type( i ) * h( 0 ) +
600  value_type( j ) * h( 1 ) )/ value_type( Order );
601  }
602  }
603 
604  return G;
605  }
606 
607  else
608  return glas::average( RefConv.vertices() );
609  }
610 
611  points_type
612  make_hexa_points( int interior = 0 )
613  {
614  if ( Order > 0 )
615  {
616  ublas::vector<node_type> h ( 3 );
617  h( 0 ) = RefConv.vertex( 1 ) - RefConv.vertex( 0 );
618  h( 1 ) = RefConv.vertex( 3 ) - RefConv.vertex( 0 );
619  h( 2 ) = RefConv.vertex( 4 ) - RefConv.vertex( 0 );
620 
621  points_type G( Dim, n_hexa_points( interior ) );
622  DVLOG(2) << "n hexa pts = " << n_hexa_points( interior ) << "\n";
623 
624  for ( int i = interior, p = 0; i < int( Order )+1-interior; ++i )
625  {
626  for ( int j = interior; j < int( Order ) + 1 - interior; ++j )
627  {
628  for ( int k = interior; k < int( Order ) + 1 - interior; ++k, ++p )
629  {
630  ublas::column( G, p ) = RefConv.vertex( 0 ) + ( value_type( i ) * h( 0 ) +
631  value_type( j ) * h( 1 ) +
632  value_type( k ) * h( 2 ) ) / value_type( Order );
633 
634  }
635  }
636  }
637 
638  return G;
639  }
640 
641  else
642  return glas::average( RefConv.vertices() );
643  }
644 
645  template<size_type shape>
646  struct pt_to_edge
647  {
648  pt_to_edge( std::vector<uint16_type> vert_ids )
649  :
650  h( 1.0 ),
651  a( Entity<SHAPE_LINE, value_type>().vertex( 0 ) ),
652  b( Entity<SHAPE_LINE, value_type>().vertex( 1 ) ),
653  u( Entity<shape, value_type>().vertex( vert_ids[ 0 ] ) ),
654  v( Entity<shape, value_type>().vertex( vert_ids[ 1 ] ) ),
655  diff( v-u )
656  {
657  h = 1.0/( b[0]-a[0] );
658  }
659  node_type
660  operator()( node_type const& x ) const
661  {
662  return u + h * ( x[ 0 ] - a[ 0 ] ) * diff;
663  }
664  value_type h;
665  node_type a, b;
666  node_type u, v, diff;
667  };
668 
669  //
670  // pt_to_face tetrahedron
671  //
672  template<size_type shape>
673  struct pt_to_face_tetrahedron
674  {
675  pt_to_face_tetrahedron( std::vector<uint16_type> vert_ids )
676  :
677  u( Entity<shape, value_type>().vertex( vert_ids[ 0 ] ) ),
678  v( Entity<shape, value_type>().vertex( vert_ids[ 1 ] ) ),
679  w( Entity<shape, value_type>().vertex( vert_ids[ 2 ] ) ),
680  diff( 2 )
681  {
682  diff[0] = v-u;
683  diff[1] = w-u;
684  }
685  node_type
686  operator()( node_type const& x ) const
687  {
688  return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
689  }
690  node_type u, v, w;
691  ublas::vector<node_type> diff;
692  };
693 
694  //
695  // pt_to_face hexa
696  //
697  template<size_type shape>
698  struct pt_to_face_hexahedron
699  {
700  pt_to_face_hexahedron( std::vector<uint16_type> vert_ids )
701  :
702  u( Entity<shape, value_type>().vertex( vert_ids[ 0 ] ) ),
703  v( Entity<shape, value_type>().vertex( vert_ids[ 1 ] ) ),
704  w( Entity<shape, value_type>().vertex( vert_ids[ 3 ] ) ),
705  diff( 2 )
706  {
707  diff[0] = v-u;
708  diff[1] = w-u;
709  }
710  node_type
711  operator()( node_type const& x ) const
712  {
713  return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
714  }
715  node_type u, v, w;
716  ublas::vector<node_type> diff;
717  };
718 
719  template<size_type shape>
720  struct pt_to_element
721  {
722  pt_to_element() {}
723  pt_to_element( std::vector<uint16_type> const& ) {}
724  node_type operator()( node_type const& x ) const
725  {
726  return x;
727  }
728  };
729 
730  template<size_type shape, uint16_type topo_dim>
731  struct pt_to_entity_tetrahedron
732  {
733  typedef typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_LINE> >,
734  mpl::identity<mpl::vector<boost::none_t,pt_to_edge<shape>,pt_to_edge<shape> > >,
735  typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_TRIANGLE> >,
736  mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_element<shape> > >,
737  mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_face_tetrahedron<shape>, pt_to_element<shape> > >
738  >::type // 2
739  >::type::type _type;
740  typedef typename mpl::at<_type, mpl::int_<topo_dim> >::type mapping_type;
741  typedef mpl::vector<boost::none_t, edge_to_point_t, face_to_point_t> list_v;
742 
743  pt_to_entity_tetrahedron( uint16_type entity_id )
744  :
745  mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
746  {}
747 
748  node_type operator()( node_type const& x ) const
749  {
750  return mapping( x );
751  }
752  mapping_type mapping;
753  };
754 
755  template<size_type shape,uint16_type topo_dim>
756  struct pt_to_entity_hexahedron
757  {
758  typedef typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_LINE> >,
759  mpl::identity<mpl::vector<boost::none_t,pt_to_edge<shape>,pt_to_edge<shape> > >,
760  typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_QUAD> >,
761  mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_element<shape> > >,
762  mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_face_hexahedron<shape>, pt_to_element<shape> > >
763  >::type // 2
764  >::type::type _type;
765  typedef typename mpl::at<_type, mpl::int_<topo_dim> >::type mapping_type;
766  typedef mpl::vector<boost::none_t, edge_to_point_t, face_to_point_t> list_v;
767 
768  pt_to_entity_hexahedron( uint16_type entity_id )
769  :
770  mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
771  {}
772 
773  node_type operator()( node_type const& x ) const
774  {
775  return mapping( x );
776  }
777  mapping_type mapping;
778  };
779 
780 };
781 
782 } // Feel
783 #endif /* __PointSetEquiSpaced_H */

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