30 #ifndef __PointSetEquiSpaced_H
31 #define __PointSetEquiSpaced_H 1
39 #include <feel/feelcore/traits.hpp>
42 #include <boost/numeric/ublas/matrix.hpp>
43 #include <boost/numeric/ublas/io.hpp>
50 namespace ublas = boost::numeric::ublas;
52 template<
class Convex, u
int16_type Order,
typename T>
53 class PointSetEquiSpaced :
public PointSet< Convex, T>
58 typedef PointSet<Convex, T> super;
62 typedef typename super::return_type return_type;
63 typedef typename super::self_type self_type;
65 typedef typename super::node_type node_type;
67 typedef typename super::nodes_type nodes_type;
69 typedef typename matrix_node<value_type>::type points_type;
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;
76 static const size_type Shape = Convex::Shape;
78 static const bool is_simplex = Convex::is_simplex;
79 static const bool is_hypercube = Convex::is_hypercube;
81 typedef mpl::if_< mpl::bool_< is_simplex >,
82 Simplex<Dim, Order, Dim> ,
83 Hypercube<Dim, Order, Dim> > conv_order_type;
85 typedef Reference<Convex, Dim, convexOrder, Dim, value_type> RefElem;
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;
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;
97 typedef std::pair<uint16_type, uint16_type> range_type;
98 typedef std::vector< std::vector<size_type> > index_map_type;
102 PointSetEquiSpaced(
int interior = 0 )
104 super( numPoints, Dim ),
107 M_eid.resize( topological_dimension + 1 );
108 M_pt_to_entity.resize( numPoints );
110 nodes_type pts( Dim, numPoints );
112 if ( interior == 0 && Order > 0 )
116 for ( uint16_type d = 0, p = 0; d < topological_dimension+1; ++d )
120 for (
int e = RefConv.entityRange( d ).begin();
121 e < RefConv.entityRange( d ).end();
124 nodes_type Gt ( makePoints( d, e ) );
128 ublas::subrange( pts, 0, Dim, p, p+Gt.size2() ) = Gt;
130 for (
size_type j = 0; j < Gt.size2(); ++j )
133 addToPtE( p+j, std::make_pair( d, e ) );
144 else if ( interior == 1 && Order > 0 )
147 else if ( Order == 0 )
148 this->
setPoints( glas::average( RefConv.vertices() ) );
150 this->setName(
"equispaced", Order );
153 ~PointSetEquiSpaced() {}
155 ublas::matrix_range<nodes_type const> pointsByEntity( uint16_type e )
const
157 FEELPP_ASSERT( M_eid[e].size() )( e ).error(
"no points defined on this entity" );
160 ublas::range( 0,Dim ),
161 ublas::range( *M_eid[e].begin(), *M_eid[e].rbegin() + 1 ) );
164 std::pair<uint16_type, uint16_type> interiorRangeById( uint16_type e, uint16_type
id )
const
166 uint16_type numEntities = 1;
169 numEntities = Convex::numVertices;
172 numEntities = Convex::numEdges;
175 numEntities = Convex::numFaces;
177 uint16_type N = M_eid[e].size()/numEntities;
179 return std::make_pair( *M_eid[e].begin() +
id*N, *M_eid[e].begin() + (
id+1 )*N );
182 ublas::matrix_range<nodes_type const> interiorPointsById( uint16_type e, uint16_type
id )
const
184 std::pair<uint16_type, uint16_type> position = interiorRangeById( e,
id );
187 ublas::range( 0, Dim ),
188 ublas::range( position.first, position.second ) );
193 uint32_type entityIds(
int i,
int j )
const
198 uint32_type numEntities(
int i )
const
200 return M_eid[i].size();
203 std::pair<uint16_type, uint16_type>
const& pointToEntity(
int p )
const
205 return M_pt_to_entity[p];
209 index_map_type entityToLocal ( uint16_type top_dim, uint16_type local_id,
bool boundary = 0 )
const
211 index_map_type indices( top_dim+1 );
213 if ( top_dim == 0 && boundary )
215 range_type pair = interiorRangeById( top_dim, local_id );
217 indices[0].push_back( pair.first );
222 if ( !boundary && top_dim > 0 )
224 if ( numEntities( top_dim ) != 0 )
225 indices[top_dim].push_back( local_id );
237 std::map<uint16_type, uint16_type> numPointsInEntity;
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;
243 for ( uint16_type i=0; i < numPointsInEntity[0] ; i++ )
246 indices[0].push_back( RefConv.e2p( local_id, i ) );
248 else if ( top_dim == 2 )
249 indices[0].push_back( RefConv.f2p( local_id, i ) );
252 if ( ( top_dim == 2 ) && ( M_eid[1].size() != 0 ) )
254 for ( uint16_type i=0; i < numPointsInEntity[1] ; i++ )
255 indices[1].push_back( RefConv.f2e( local_id, i ) );
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 );
265 if ( M_eid[top_dim].size() )
266 indices[top_dim].push_back( local_id );
275 points_type pointsBySubEntity( uint16_type top_dim, uint16_type local_id,
bool boundary = 0 )
const
277 index_map_type index_list = entityToLocal( top_dim, local_id, boundary );
279 uint16_type matrix_size = 0;
281 if ( index_list[0].size() != 0 )
282 matrix_size +=index_list[0].size()*nbPtsPerVertex;
284 if ( ( top_dim >= 1 ) && ( index_list[1].size() != 0 ) )
285 matrix_size +=index_list[1].size()*nbPtsPerEdge;
287 if ( ( top_dim >= 2 ) && ( index_list[2].size() != 0 ) )
288 matrix_size +=index_list[2].size()*nbPtsPerFace;
290 if ( ( top_dim == 3 ) && ( index_list[3].size() != 0 ) )
291 matrix_size +=nbPtsPerVolume;
293 points_type G ( Dim, matrix_size );
295 for ( uint16_type i=0, p=0; i < top_dim+1; i++ )
297 if ( index_list[i].size() )
299 for ( uint16_type j=0; j < index_list[i].size(); j++ )
301 points_type aux = interiorPointsById( i, index_list[i][j] );
303 ublas::subrange( G, 0, Dim, p, p+aux.size2() ) = aux;
313 index_map_type getEid ()
318 std::vector<range_type> getPtE()
320 return M_pt_to_entity;
323 void setEid ( index_map_type eid )
328 void setPtE ( std::vector<range_type> pt_ent )
330 M_pt_to_entity = pt_ent;
333 void addToEid ( uint16_type p, uint16_type q )
335 M_eid[p].push_back( q );
338 void addToPtE ( uint16_type p, range_type q )
340 M_pt_to_entity[p] = q;
347 index_map_type M_eid;
348 std::vector<range_type> M_pt_to_entity;
350 points_type makePoints( uint16_type topo_dim, uint16_type __id,
int interior = 1 )
355 points_type G( RefConv.vertices().size1(), 1 );
356 ublas::column( G, 0 ) = ublas::column( RefConv.vertices(), __id );
361 else if ( topo_dim == topological_dimension )
364 return makeLattice<Shape>( interior );
366 throw std::logic_error(
"cannot make those points" );
367 return points_type();
378 G = makeLattice<SHAPE_LINE>( 1 );
379 Gret.resize( nRealDim, G.size2() );
383 pt_to_entity_tetrahedron<Shape, 1> p_to_e( __id );
385 for (
size_type i = 0; i < G.size2(); ++i )
386 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
391 pt_to_entity_hexahedron<Shape, 1> p_to_e( __id );
393 for (
size_type i = 0; i < G.size2(); ++i )
394 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
400 else if ( topo_dim == 2 )
404 G = makeLattice<SHAPE_TRIANGLE>( 1 );
405 Gret.resize( nRealDim, G.size2() );
406 pt_to_entity_tetrahedron<Shape, 2> p_to_e( __id );
408 for (
size_type i = 0; i < G.size2(); ++i )
409 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
414 G = makeLattice<SHAPE_QUAD>( 1 );
415 Gret.resize( nRealDim, G.size2() );
416 pt_to_entity_hexahedron<Shape, 2> p_to_e( __id );
418 for (
size_type i = 0; i < G.size2(); ++i )
419 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
426 return points_type();
429 template<
size_type shape>
430 points_type makeLattice( uint16_type interior = 0 )
436 if ( shape == SHAPE_LINE )
437 G = make_line_points( interior );
439 else if ( shape == SHAPE_TRIANGLE )
440 G = make_triangle_points( interior );
442 else if ( shape == SHAPE_TETRA )
443 G = make_tetrahedron_points( interior );
445 else if ( shape == SHAPE_QUAD )
446 return make_quad_points( interior );
448 else if ( shape == SHAPE_HEXA )
449 return make_hexa_points( interior );
452 else if ( Order == 0 )
453 G = glas::average( RefConv.vertices() );
459 int n_line_points(
int interior = 0 )
461 return std::max( 0,
int( Order )+1-2*interior );
463 int n_triangle_points(
int interior = 0 )
466 return std::max( 0, (
int( Order )+1-2*interior )*(
int( Order )-2*interior )/2 );
468 return ( Order+1 )*( Order+2 )/2;
470 int n_tetrahedron_points(
int interior = 0 )
473 return std::max( 0, (
int( Order )+1-2*interior )*(
int( Order )-2*interior )*(
int( Order )-1-2*interior )/6 );
475 return ( Order+1 )*( Order+2 )*( Order+3 )/6;
478 int n_quad_points(
int interior = 0 )
const
481 return std::max( 0, (
int( Order )+1-2*interior )*(
int( Order )+1-2*interior ) );
483 return ( Order+1 )*( Order+1 );
486 int n_hexa_points(
int interior = 0 )
const
489 return std::max( 0, (
int( Order )+1-2*interior )*(
int( Order )+1-2*interior )*(
int( Order )+1-2*interior ) );
491 return ( Order+1 )*( Order+1 )*( Order+1 );
495 make_line_points(
int interior = 0 )
501 ublas::vector<node_type> h ( 1 );
502 h( 0 ) = RefConv.vertex( 1 ) - RefConv.vertex( 0 );
504 G.resize( Dim, n_line_points( interior ) );
506 for (
int i = interior, indp = 0; i < int( Order )+1-interior; ++i, ++indp )
508 ublas::column( G, indp ) = RefConv.vertex( 0 ) + ( h( 0 ) * value_type( i ) )/value_type( Order );
513 G = glas::average( RefConv.vertices() );
520 make_triangle_points(
int interior = 0 )
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 );
530 G.resize( Dim, n_triangle_points( interior ) );
532 for (
int i = interior, p = 0; i < int( Order )+1-interior; ++i )
534 for (
int j = interior; j < int( Order ) + 1 - i-interior; ++j, ++p )
536 ublas::column( G, p ) = RefConv.vertex( 0 ) + ( value_type( i ) * h( 1 ) +
537 value_type( j ) * h( 0 ) )/ value_type( Order );
543 G = glas::average( RefConv.vertices() );
549 make_tetrahedron_points(
int interior = 0 )
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 );
560 G.resize( Dim, n_tetrahedron_points( interior ) );
562 for (
int i = interior, p = 0; i < int( Order )+1-interior; ++i )
564 for (
int j = interior; j < int( Order ) + 1 - i - interior; ++j )
566 for (
int k = interior; k < int( Order ) + 1 - i - j - interior; ++k, ++p )
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 );
578 G = glas::average( RefConv.vertices() );
584 make_quad_points(
int interior = 0 )
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 );
592 DVLOG(2) <<
"n quad pts = " << n_quad_points( interior ) <<
"\n";
593 points_type G( Dim, n_quad_points( interior ) );
595 for (
int i = interior, p = 0; i < int( Order )+1-interior; ++i )
597 for (
int j = interior; j < int( Order ) + 1 -interior; ++j, ++p )
599 ublas::column( G, p ) = RefConv.vertex( 0 ) + ( value_type( i ) * h( 0 ) +
600 value_type( j ) * h( 1 ) )/ value_type( Order );
608 return glas::average( RefConv.vertices() );
612 make_hexa_points(
int interior = 0 )
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 );
621 points_type G( Dim, n_hexa_points( interior ) );
622 DVLOG(2) <<
"n hexa pts = " << n_hexa_points( interior ) <<
"\n";
624 for (
int i = interior, p = 0; i < int( Order )+1-interior; ++i )
626 for (
int j = interior; j < int( Order ) + 1 - interior; ++j )
628 for (
int k = interior; k < int( Order ) + 1 - interior; ++k, ++p )
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 );
642 return glas::average( RefConv.vertices() );
645 template<
size_type shape>
648 pt_to_edge( std::vector<uint16_type> vert_ids )
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 ] ) ),
657 h = 1.0/( b[0]-a[0] );
660 operator()( node_type
const& x )
const
662 return u + h * ( x[ 0 ] - a[ 0 ] ) * diff;
666 node_type u, v, diff;
672 template<
size_type shape>
673 struct pt_to_face_tetrahedron
675 pt_to_face_tetrahedron( std::vector<uint16_type> vert_ids )
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 ] ) ),
686 operator()( node_type
const& x )
const
688 return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
691 ublas::vector<node_type> diff;
697 template<
size_type shape>
698 struct pt_to_face_hexahedron
700 pt_to_face_hexahedron( std::vector<uint16_type> vert_ids )
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 ] ) ),
711 operator()( node_type
const& x )
const
713 return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
716 ublas::vector<node_type> diff;
719 template<
size_type shape>
723 pt_to_element( std::vector<uint16_type>
const& ) {}
724 node_type operator()( node_type
const& x )
const
730 template<
size_type shape, u
int16_type topo_dim>
731 struct pt_to_entity_tetrahedron
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> > >
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;
743 pt_to_entity_tetrahedron( uint16_type entity_id )
745 mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
748 node_type operator()( node_type
const& x )
const
752 mapping_type mapping;
755 template<
size_type shape,u
int16_type topo_dim>
756 struct pt_to_entity_hexahedron
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> > >
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;
768 pt_to_entity_hexahedron( uint16_type entity_id )
770 mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
773 node_type operator()( node_type
const& x )
const
777 mapping_type mapping;
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
boost::tuple< mpl::size_t< MESH_POINTS >, typename MeshTraits< MeshType >::point_const_iterator, typename MeshTraits< MeshType >::point_const_iterator > points(MeshType const &mesh)
Definition: filters.hpp:1296
void setPoints(nodes_type const &pts)
Definition: pointset.hpp:267
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
#define FEELPP_DEFINE_VISITABLE()
Definition: visitor.hpp:319