29 #ifndef __PointSetInterpolation_H
30 #define __PointSetInterpolation_H 1
33 #include <boost/numeric/ublas/matrix.hpp>
34 #include <boost/numeric/ublas/io.hpp>
38 namespace ublas = boost::numeric::ublas;
48 template<uint16_type Dim,
51 template<u
int16_type,u
int16_type,u
int16_type>
class Convex = Simplex>
62 typedef typename super::return_type return_type;
66 typedef typename super::nodes_type nodes_type;
69 static const bool is_simplex = convex_type::is_simplex;
70 static const bool is_hypercube = convex_type::is_hypercube;
72 static const uint32_type convexOrder = convex_type::nOrder;
73 static const uint32_type topological_dimension = convex_type::topological_dimension;
75 typedef mpl::if_< mpl::bool_< is_simplex >,
77 Hypercube<Dim, Order, Dim> > conv_order_type;
81 static const uint32_type numPoints = conv_order_type::type::numPoints;
83 static const uint32_type nbPtsPerVertex = conv_order_type::type::nbPtsPerVertex;
84 static const uint32_type nbPtsPerEdge = conv_order_type::type::nbPtsPerEdge;
85 static const uint32_type nbPtsPerFace = conv_order_type::type::nbPtsPerFace;
86 static const uint32_type nbPtsPerVolume = conv_order_type::type::nbPtsPerVolume;
90 typedef std::pair<uint16_type, uint16_type> range_type;
91 typedef std::vector< std::vector<size_type> > index_map_type;
104 M_eid.resize( topological_dimension + 1 );
105 M_pt_to_entity.resize( numPoints );
111 M_eid.resize( topological_dimension + 1 );
112 M_pt_to_entity.resize( numPoints );
117 M_eid( psi.getEid() ),
118 M_pt_to_entity( psi.getPtE() )
124 ublas::matrix_range<nodes_type const> pointsByEntity( uint16_type e )
const
126 FEELPP_ASSERT( M_eid[e].size() )( e ).error(
"no points defined on this entity" );
129 ublas::range( 0,Dim ),
130 ublas::range( *M_eid[e].begin(), *M_eid[e].rbegin() + 1 ) );
133 range_type interiorRangeById( uint16_type e, uint16_type
id )
const
135 uint16_type numEntities = 1;
138 numEntities = convex_type::numVertices;
141 numEntities = convex_type::numEdges;
144 numEntities = convex_type::numFaces;
146 uint16_type N = M_eid[e].size()/numEntities;
148 return std::make_pair( *M_eid[e].begin() +
id*N, *M_eid[e].begin() + (
id+1 )*N );
151 ublas::matrix_range<nodes_type const> interiorPointsById( uint16_type e, uint16_type
id )
const
153 range_type position = interiorRangeById( e,
id );
155 ublas::matrix_range<nodes_type const> G =
ublas::project( this->points(),
156 ublas::range( 0, Dim ),
157 ublas::range( position.first, position.second ) );
162 uint32_type entityIds(
int i,
int j )
const
167 uint32_type numEntities(
int i )
const
169 return M_eid[i].size();
172 range_type
const& pointToEntity(
int p )
const
174 return M_pt_to_entity[p];
178 index_map_type entityToLocal ( uint16_type top_dim, uint16_type local_id,
bool boundary = 0 )
180 index_map_type indices( top_dim+1 );
182 if ( top_dim == 0 && boundary )
184 range_type pair = interiorRangeById( top_dim, local_id );
186 indices[0].push_back( pair.first );
191 if ( !boundary && top_dim > 0 )
193 if ( numEntities( top_dim ) != 0 )
194 indices[top_dim].push_back( local_id );
206 std::map<uint16_type, uint16_type> numPointsInEntity;
208 numPointsInEntity[0] = ( is_simplex )?( top_dim+1 ):( 2 + ( top_dim-1 )*( 2+3*( top_dim-2 ) ) );
209 numPointsInEntity[1] = ( 2 + ( top_dim-1 )*( 1+2*is_simplex + ( 5-4*is_simplex )*( top_dim-1 ) ) )/2;
210 numPointsInEntity[2] = 6-2*is_simplex;
212 for ( uint16_type i=0; i < numPointsInEntity[0] ; i++ )
215 indices[0].push_back( RefConv.e2p( local_id, i ) );
217 else if ( top_dim == 2 )
218 indices[0].push_back( RefConv.f2p( local_id, i ) );
221 if ( ( top_dim == 2 ) && ( M_eid[1].size() != 0 ) )
223 for ( uint16_type i=0; i < numPointsInEntity[1] ; i++ )
224 indices[1].push_back( RefConv.f2e( local_id, i ) );
229 for ( uint16_type k=0; k<numPointsInEntity.size(); k++ )
231 for ( uint16_type i=0; i < numPointsInEntity[k]; i++ )
233 indices[k].push_back( i );
238 if ( M_eid[top_dim].size() )
240 indices[top_dim].push_back( local_id );
250 points_type pointsBySubEntity( uint16_type top_dim, uint16_type local_id,
bool boundary = 0 )
252 index_map_type index_list = entityToLocal( top_dim, local_id, boundary );
254 uint16_type matrix_size = 0;
256 if ( index_list[0].size() != 0 )
257 matrix_size +=index_list[0].size()*nbPtsPerVertex;
259 if ( ( top_dim>=1 ) && ( index_list[1].size() != 0 ) )
260 matrix_size +=index_list[1].size()*nbPtsPerEdge;
262 if ( ( top_dim >=2 ) && ( index_list[2].size() != 0 ) )
263 matrix_size +=index_list[2].size()*nbPtsPerFace;
265 if ( ( top_dim ==3 ) && ( index_list[3].size() != 0 ) )
266 matrix_size +=nbPtsPerVolume;
268 points_type G ( Dim, matrix_size );
270 for ( uint16_type i=0, p=0; i < top_dim+1; i++ )
272 if ( index_list[i].size() )
274 for ( uint16_type j=0; j < index_list[i].size(); j++ )
276 points_type aux = interiorPointsById( i, index_list[i][j] );
278 ublas::subrange( G, 0, Dim, p, p+aux.size2() ) = aux;
288 index_map_type getEid ()
293 std::vector<range_type> getPtE()
295 return M_pt_to_entity;
298 void setEid ( index_map_type eid )
303 void setPtE ( std::vector<range_type> pt_ent )
305 M_pt_to_entity = pt_ent;
308 void addToEid ( uint16_type p, uint16_type q )
310 M_eid[p].push_back( q );
313 void addToPtE ( uint16_type p, range_type q )
315 M_pt_to_entity[p] = q;
318 FEELPP_DEFINE_VISITABLE();
322 index_map_type M_eid;
323 std::vector<range_type> M_pt_to_entity;