29 #ifndef __PointSetMapped_H
30 #define __PointSetMapped_H 1
34 #include <boost/numeric/ublas/matrix_sparse.hpp>
35 #include <boost/numeric/ublas/io.hpp>
38 #include <feel/feelcore/feel.hpp>
40 #include <feel/feelcore/traits.hpp>
45 #include <feel/feelmesh/geoelement.hpp>
55 namespace ublas = boost::numeric::ublas;
57 template<
typename element_type,
61 template<
class, u
int16_type,
class>
class PointSetType = PointSetEquiSpaced>
62 class PointSetMapped :
public PointSetType<Convex, Order, T>
66 typedef PointSetType<Convex, Order, T> pointset_type;
70 static const uint32_type Dim = Convex::nDim;
71 static const uint32_type convexOrder = Convex::nOrder;
73 static const bool is_simplex = Convex::is_simplex;
75 typedef typename pointset_type::nodes_type nodes_type;
76 typedef typename matrix_node<value_type>::type points_type;
78 typedef typename element_type::gm_type gm_type;
79 typedef typename element_type::gm_ptrtype gm_ptrtype;
81 typedef typename element_type::edge_permutation_type edge_permutation_type;
82 typedef typename element_type::face_permutation_type face_permutation_type;
84 typedef ublas::vector<uint16_type> permutation_vector_type;
85 typedef ublas::mapped_matrix<uint16_type> permutation_matrix_type;
87 typedef mpl::if_< mpl::bool_< is_simplex >,
88 Simplex<Dim, Order, Dim> ,
89 Hypercube<Dim, Order, Dim> > conv_order_type;
91 static const uint32_type nbPtsPerVertex = conv_order_type::type::nbPtsPerVertex;
92 static const uint32_type nbPtsPerEdge = conv_order_type::type::nbPtsPerEdge;
93 static const uint32_type nbPtsPerFace = conv_order_type::type::nbPtsPerFace;
95 typedef Reference<Convex, Dim, convexOrder, Dim, value_type> RefElem;
97 typedef typename pointset_type::range_type range_type;
98 typedef typename pointset_type::index_map_type index_map_type;
102 PointSetMapped( element_type
const& _elt )
108 points_type Gt = updatePoints<2>( updatePoints<1>( pts.points(), pts, mpl::bool_< ( Order > 2 ) >() ),
110 mpl::bool_< ( Dim == 3 && Order > 2+uint16_type( is_simplex ) ) >() );
112 this->setPoints( Gt );
117 permutation_vector_type getVectorPermutation ( face_permutation_type P )
119 return vector_permutation[P];
122 permutation_matrix_type getMatrixPermutation ( face_permutation_type P )
124 return matrix_permutation[P];
127 points_type pointsBySubEntity( uint16_type top_dim, uint16_type local_id,
bool boundary = 0,
bool real = 0 )
129 index_map_type index_list = this->entityToLocal( top_dim, local_id, boundary );
133 uint16_type matrix_size = 0;
135 if ( index_list[0].size() != 0 )
136 matrix_size +=index_list[0].size()*nbPtsPerVertex;
138 if ( ( top_dim>=1 ) && ( index_list[1].size() != 0 ) )
139 matrix_size +=index_list[1].size()*nbPtsPerEdge;
141 if ( ( top_dim>=2 ) && ( index_list[2].size() != 0 ) )
142 matrix_size +=index_list[2].size()*nbPtsPerFace;
144 points_type G ( Dim, matrix_size );
146 for ( uint16_type i=0, p=0; i < top_dim+1; i++ )
148 if ( index_list[i].size() )
150 for ( uint16_type j=0; j < index_list[i].size(); j++ )
155 points_type aux = this->interiorPointsById( i, index_list[i][j] );
157 ublas::subrange( G, 0, Dim, p, p+aux.size2() ) = aux;
165 G = transformToReal( G );
174 std::map<face_permutation_type, permutation_vector_type> vector_permutation;
175 std::map<face_permutation_type, permutation_matrix_type> matrix_permutation;
177 template <
int top_dim>
178 nodes_type updatePoints( nodes_type
const& pts, pointset_type
const& , mpl::bool_<false> )
183 template <
int top_dim>
184 nodes_type updatePoints( nodes_type
const& thepts, pointset_type
const& G, mpl::bool_<true> )
186 nodes_type pts( thepts );
190 generateFacePermutations( Order - 1-uint16_type( is_simplex ),
191 mpl::int_<Dim>(), mpl::bool_< is_simplex >() );
194 for (
size_type e = RefConv.entityRange( top_dim ).begin();
195 e < RefConv.entityRange( top_dim ).end();
198 points_type Gt = G.pointsBySubEntity( top_dim, e );
200 Gt = permutatePoints ( Gt, e, mpl::int_<top_dim>() );
202 uint16_type pos = G.interiorRangeById( top_dim, e ).first;
204 ublas::subrange( pts, 0, Dim, pos, pos+Gt.size2() ) = Gt;
212 points_type permutatePoints ( nodes_type
const& theGt,
size_type entity_local_id, mpl::int_<1> )
214 nodes_type Gt( theGt );
215 edge_permutation_type permutation = M_elt.edgePermutation( entity_local_id );
217 if ( permutation != edge_permutation_type( edge_permutation_type::IDENTITY ) )
219 for ( uint16_type i=0; i <= ( Gt.size2()-1 )/2 - Gt.size2()%2; i++ )
220 ublas::column( Gt, i ).swap( ublas::column( Gt, Gt.size2()-1-i ) );
228 points_type permutatePoints ( nodes_type
const& Gt,
size_type entity_local_id, mpl::int_<2> )
230 face_permutation_type permutation = M_elt.facePermutation( entity_local_id );
231 nodes_type res( Gt );
233 if ( permutation != face_permutation_type( face_permutation_type::IDENTITY ) )
234 res = prod( Gt, getMatrixPermutation( permutation ) );
239 permutation_matrix_type vectorToMatrixPermutation ( permutation_vector_type
const& v )
241 permutation_matrix_type P ( v.size(), v.size(), v.size()*v.size() );
243 for ( uint16_type i = 0; i < v.size(); ++i )
249 permutation_vector_type matrixToVectorPermutation ( permutation_matrix_type
const& P )
251 FEELPP_ASSERT( P.size1() == P.size2() ).error(
"invalid permutation" );
253 permutation_vector_type v ( P.size1() );
255 for ( uint16_type i = 0; i < v.size(); ++i )
256 for ( uint16_type j = 0; j < v.size(); ++j )
263 void setPermutation(
int out,
int first,
int second )
265 matrix_permutation[face_permutation_type( out )] = ublas::prod( matrix_permutation[face_permutation_type( first )],
266 matrix_permutation[face_permutation_type( second )] );
268 vector_permutation[face_permutation_type( out )] = matrixToVectorPermutation( matrix_permutation[face_permutation_type( out )] );
271 void generateFacePermutations( uint16_type , mpl::int_<2>, mpl::bool_<true> ) {}
272 void generateFacePermutations( uint16_type , mpl::int_<2>, mpl::bool_<false> ) {}
275 void generateFacePermutations( uint16_type n_side_points, mpl::int_<3>, mpl::bool_<true> )
277 uint16_type npoints = n_side_points*( n_side_points+1 )/2;
279 permutation_vector_type _vec( npoints );
282 for ( uint16_type i = 0; i < n_side_points; ++i )
284 for ( uint16_type j = 0; j <= i; j++ )
286 uint16_type _first = i + j*( j-1 )/2 + j*( n_side_points - j );
287 uint16_type _last = i + ( i-j )*( i-j-1 )/2 + ( i-j )*( n_side_points - ( i-j ) );
289 _vec( _first ) = _last;
293 vector_permutation[face_permutation_type( face_permutation_type::REVERSE_HYPOTENUSE )] = _vec;
294 matrix_permutation[face_permutation_type( face_permutation_type::REVERSE_HYPOTENUSE )] = vectorToMatrixPermutation( _vec );
297 for ( uint16_type i = 0; i < n_side_points; ++i )
299 uint16_type _begin = i*n_side_points - i*( i-1 )/2;
300 uint16_type _end = ( i+1 )*n_side_points - ( i+1 )*i/2 - 1;
302 for ( uint16_type j = 0; j <= n_side_points - i - 1; j++ )
303 _vec( _begin + j ) = _end - j;
306 vector_permutation[face_permutation_type( face_permutation_type::REVERSE_BASE )] = _vec;
307 matrix_permutation[face_permutation_type( face_permutation_type::REVERSE_BASE )] = vectorToMatrixPermutation( _vec );
309 setPermutation( face_permutation_type::ROTATION_ANTICLOCK,
310 face_permutation_type::REVERSE_BASE,
311 face_permutation_type::REVERSE_HYPOTENUSE );
313 setPermutation( face_permutation_type::ROTATION_CLOCKWISE,
314 face_permutation_type::REVERSE_HYPOTENUSE,
315 face_permutation_type::REVERSE_BASE );
317 setPermutation( face_permutation_type::REVERSE_HEIGHT,
318 face_permutation_type::REVERSE_BASE,
319 face_permutation_type::ROTATION_CLOCKWISE );
323 void generateFacePermutations( uint16_type n_side_points, mpl::int_<3>, mpl::bool_<false> )
325 uint16_type npoints = n_side_points*( n_side_points );
327 permutation_vector_type _vec( npoints );
332 for ( int16_type i = n_side_points-1; i >= 0; --i )
334 uint16_type _first = i*n_side_points;
336 for ( uint16_type j = 0; j <= n_side_points-1; j++ )
338 _vec( p ) = _first + j;
343 vector_permutation[face_permutation_type( face_permutation_type::REVERSE_BASE )] = _vec;
344 matrix_permutation[face_permutation_type( face_permutation_type::REVERSE_BASE )] = vectorToMatrixPermutation( _vec );
349 for ( int16_type i = n_side_points-1; i >= 0; --i )
351 for ( uint16_type j = 0; j <= n_side_points-1; j++ )
353 _vec( p ) = i + n_side_points*j;
358 vector_permutation[face_permutation_type( face_permutation_type::ROTATION_ANTICLOCK )] = _vec;
359 matrix_permutation[face_permutation_type( face_permutation_type::ROTATION_ANTICLOCK )] = vectorToMatrixPermutation( _vec );
361 setPermutation( face_permutation_type::SECOND_DIAGONAL,
362 face_permutation_type::REVERSE_BASE,
363 face_permutation_type::ROTATION_ANTICLOCK );
365 setPermutation( face_permutation_type::REVERSE_HEIGHT,
366 face_permutation_type::SECOND_DIAGONAL,
367 face_permutation_type::ROTATION_ANTICLOCK );
369 setPermutation( face_permutation_type::ROTATION_TWICE_CLOCKWISE,
370 face_permutation_type::REVERSE_HEIGHT,
371 face_permutation_type::REVERSE_BASE );
373 setPermutation( face_permutation_type::PRINCIPAL_DIAGONAL,
374 face_permutation_type::REVERSE_HEIGHT,
375 face_permutation_type::ROTATION_ANTICLOCK );
377 setPermutation( face_permutation_type::ROTATION_CLOCKWISE,
378 face_permutation_type::ROTATION_ANTICLOCK,
379 face_permutation_type::ROTATION_TWICE_CLOCKWISE );
382 points_type transformToReal( points_type
const& Gt )
384 gm_ptrtype _gm_ptr(
new gm_type );
386 typename gm_type::precompute_ptrtype __geopc(
new typename gm_type::precompute_type( _gm_ptr, Gt ) );
388 typename gm_type::template Context<vm::POINT, element_type> gmc( _gm_ptr, M_elt, __geopc );