29 #ifndef __Quadmapped_H
30 #define __Quadmapped_H 1
38 #include <boost/numeric/ublas/matrix.hpp>
39 #include <boost/numeric/ublas/io.hpp>
41 #include <feel/feelcore/traits.hpp>
44 #include <feel/feelmesh/geoelement.hpp>
51 namespace ublas = boost::numeric::ublas;
58 template<
typename PsetType>
59 class QuadMapped :
public PsetType
61 typedef PsetType super;
63 typedef super pointset_type;
65 typedef typename super::value_type value_type;
67 typedef typename super::convex_type convex_type;
68 static const uint32_type Dim = convex_type::nDim;
69 static const uint32_type convexOrder = convex_type::nOrder;
71 static const bool is_simplex = convex_type::is_simplex;
73 typedef typename pointset_type::nodes_type nodes_type;
74 typedef typename matrix_node<value_type>::type points_type;
78 typedef ublas::vector<uint16_type> permutation_vector_type;
79 typedef ublas::mapped_matrix<uint16_type> permutation_matrix_type;
82 typedef convex_type conv_order_type;
85 static const uint32_type nbPtsPerVertex = conv_order_type::type::nbPtsPerVertex;
86 static const uint32_type nbPtsPerEdge = conv_order_type::type::nbPtsPerEdge;
87 static const uint32_type nbPtsPerFace = conv_order_type::type::nbPtsPerFace;
89 static const uint32_type nbPtsPerVertex = conv_order_type::nbPtsPerVertex;
90 static const uint32_type nbPtsPerEdge = conv_order_type::nbPtsPerEdge;
91 static const uint32_type nbPtsPerFace = conv_order_type::nbPtsPerFace;
93 typedef typename convex_type::vertex_permutation_type vertex_permutation_type;
94 typedef typename convex_type::edge_permutation_type edge_permutation_type;
95 typedef typename convex_type::face_permutation_type face_permutation_type;
96 typedef typename mpl::if_<mpl::equal_to<mpl::int_<Dim>, mpl::int_<2> >, mpl::identity<edge_permutation_type>, mpl::identity<face_permutation_type> >::type::type permutation_type;
97 typedef std::vector<std::map<permutation_type, points_type> > permutation_points_type;
107 permutation_points_type
108 operator()( pointset_type
const& )
110 permutation_points_type ppts( convex_type::numTopologicalFaces );
114 for (
int f = 0; f < convex_type::numTopologicalFaces; ++f )
116 __typeof__(
typename pointset_type::Face( *
this, f ) ) facequad = typename pointset_type::Face( *this, f );
118 generate( ppts, facequad, f, mpl::int_<Dim>() );
126 permutation_vector_type getVectorPermutation ( face_permutation_type P )
128 return vector_permutation[P];
131 permutation_matrix_type getMatrixPermutation ( face_permutation_type P )
133 return matrix_permutation[P];
137 generate( permutation_points_type& ppts,
typename pointset_type::Face
const& pts,
int f, mpl::int_<1> )
139 vertex_permutation_type pbegin( vertex_permutation_type::IDENTITY );
140 vertex_permutation_type pend( vertex_permutation_type::N_PERMUTATIONS );
142 for ( vertex_permutation_type p = pbegin; p < pend; ++p )
144 ppts[f].insert( std::make_pair( p, pts.points() ) );
145 DVLOG(2) <<
"[quadmapped] ppts[" << f <<
"]["
147 << ppts[f].find( p )->second <<
"\n";
151 generate( permutation_points_type& ppts,
typename pointset_type::Face
const& pts,
int f, mpl::int_<2> )
153 edge_permutation_type pbegin( edge_permutation_type::IDENTITY );
154 edge_permutation_type pend( edge_permutation_type::N_PERMUTATIONS );
156 for ( edge_permutation_type p = pbegin; p < pend; ++p )
158 ppts[f].insert( std::make_pair( p, permutatePoints( pts.points(), p ) ) );
159 DVLOG(2) <<
"[quadmapped] ppts[" << f <<
"][" << p.value() <<
"]=" << ppts[f].find( p )->second <<
"\n";
163 generate( permutation_points_type& ppts,
typename pointset_type::Face
const& pts,
int f, mpl::int_<3> )
165 uint16_type npts = pts.nPoints();
168 generateFacePermutations( npts, mpl::bool_<is_simplex>() );
170 face_permutation_type pbegin( face_permutation_type::IDENTITY );
171 face_permutation_type pend( face_permutation_type::N_PERMUTATIONS );
173 for ( face_permutation_type p = pbegin; p < pend; ++p )
175 ppts[f].insert( std::make_pair( p, permutatePoints( pts.points(), p ) ) );
176 DVLOG(2) <<
"[quadmapped] ppts[" << f <<
"][" << p.value() <<
"]=" << ppts[f].find( p )->second <<
"\n";
182 points_type permutatePoints ( points_type Gt,
183 edge_permutation_type edge_perm )
185 if ( edge_perm != edge_permutation_type( edge_permutation_type::IDENTITY ) )
188 for (
int i=0; i <= ( int( Gt.size2() )-1 )/2 - int( Gt.size2() )%2; i++ )
189 ublas::column( Gt, i ).swap( ublas::column( Gt, Gt.size2()-1-i ) );
197 points_type permutatePoints ( nodes_type
const& Gt,
198 face_permutation_type face_perm )
200 FEELPP_ASSERT( face_perm != face_permutation_type( face_permutation_type::NO_PERMUTATION ) )
201 ( face_perm ).error(
"invalid permutation" );
203 nodes_type res( Gt );
205 if ( face_perm != face_permutation_type( face_permutation_type::IDENTITY ) )
206 res = prod( Gt, getMatrixPermutation( face_perm ) );
208 FEELPP_ASSERT( res.size1() == Gt.size1() &&
209 res.size2() == Gt.size2() ) ( res )( Gt ).error (
"invalid permutation operation" );
213 permutation_matrix_type vectorToMatrixPermutation ( permutation_vector_type
const& v )
215 permutation_matrix_type P ( v.size(), v.size(), v.size()*v.size() );
218 for ( uint16_type i = 0; i < v.size(); ++i )
224 permutation_vector_type matrixToVectorPermutation ( permutation_matrix_type
const& P )
226 FEELPP_ASSERT( P.size1() == P.size2() ).error(
"invalid permutation" );
228 permutation_vector_type v ( P.size1() );
231 for ( uint16_type i = 0; i < v.size(); ++i )
232 for ( uint16_type j = 0; j < v.size(); ++j )
239 void setPermutation( face_permutation_type out, face_permutation_type first, face_permutation_type second )
241 matrix_permutation[out] = ublas::prod( matrix_permutation[first],
242 matrix_permutation[second] );
244 vector_permutation[out] = matrixToVectorPermutation( matrix_permutation[out] );
249 void generateFacePermutations( uint16_type npoints, mpl::bool_<true> )
253 uint16_type n_side_points = ( uint16_type )math::sqrt( (
double )npoints );
257 permutation_vector_type _vec( npoints );
261 for ( uint16_type i = 0; i < n_side_points; ++i )
263 for ( uint16_type j = 0; j <= i; j++ )
265 uint16_type _first = i + j*( j-1 )/2 + j*( n_side_points - j );
266 uint16_type _last = i + ( i-j )*( i-j-1 )/2 + ( i-j )*( n_side_points - ( i-j ) );
268 _vec( _first ) = _last;
272 vector_permutation[face_permutation_type( face_permutation_type::REVERSE_HYPOTENUSE )] = _vec;
273 matrix_permutation[face_permutation_type( face_permutation_type::REVERSE_HYPOTENUSE )] = vectorToMatrixPermutation( _vec );
276 for ( uint16_type i = 0; i < n_side_points; ++i )
278 uint16_type _begin = i*n_side_points - i*( i-1 )/2;
279 uint16_type _end = ( i+1 )*n_side_points - ( i+1 )*i/2 - 1;
281 for ( uint16_type j = 0; j <= n_side_points - i - 1; j++ )
282 _vec( _begin + j ) = _end - j;
285 vector_permutation[face_permutation_type( face_permutation_type::REVERSE_BASE )] = _vec;
286 matrix_permutation[face_permutation_type( face_permutation_type::REVERSE_BASE )] = vectorToMatrixPermutation( _vec );
288 setPermutation( ( face_permutation_type )face_permutation_type::ROTATION_ANTICLOCK,
289 ( face_permutation_type )face_permutation_type::REVERSE_BASE,
290 ( face_permutation_type )face_permutation_type::REVERSE_HYPOTENUSE );
292 setPermutation( ( face_permutation_type )face_permutation_type::ROTATION_CLOCKWISE,
293 ( face_permutation_type )face_permutation_type::REVERSE_HYPOTENUSE,
294 ( face_permutation_type )face_permutation_type::REVERSE_BASE );
296 setPermutation( ( face_permutation_type )face_permutation_type::REVERSE_HEIGHT,
297 ( face_permutation_type )face_permutation_type::REVERSE_BASE,
298 ( face_permutation_type )face_permutation_type::ROTATION_CLOCKWISE );
302 void generateFacePermutations( uint16_type npoints, mpl::bool_<false> )
305 uint16_type n_side_points = ( uint16_type )math::sqrt( (
double )npoints );
306 FEELPP_ASSERT( npoints == n_side_points*n_side_points )
307 ( npoints )( n_side_points ).error(
"invalid number of points" );
308 permutation_vector_type _vec( npoints );
313 for ( int16_type i = n_side_points-1; i >= 0; --i )
315 uint16_type _first = i*n_side_points;
317 for ( uint16_type j = 0; j <= n_side_points-1; j++ )
319 _vec( p ) = _first + j;
324 vector_permutation[( face_permutation_type )face_permutation_type::REVERSE_BASE] = _vec;
325 matrix_permutation[( face_permutation_type )face_permutation_type::REVERSE_BASE] = vectorToMatrixPermutation( _vec );
330 for ( int16_type i = n_side_points-1; i >= 0; --i )
332 for ( uint16_type j = 0; j <= n_side_points-1; j++ )
334 _vec( p ) = i + n_side_points*j;
339 vector_permutation[( face_permutation_type )face_permutation_type::ROTATION_ANTICLOCK] = _vec;
340 matrix_permutation[( face_permutation_type )face_permutation_type::ROTATION_ANTICLOCK] = vectorToMatrixPermutation( _vec );
342 setPermutation( ( face_permutation_type )face_permutation_type::SECOND_DIAGONAL,
343 ( face_permutation_type )face_permutation_type::REVERSE_BASE,
344 ( face_permutation_type )face_permutation_type::ROTATION_ANTICLOCK );
346 setPermutation( ( face_permutation_type )face_permutation_type::REVERSE_HEIGHT,
347 ( face_permutation_type )face_permutation_type::SECOND_DIAGONAL,
348 ( face_permutation_type )face_permutation_type::ROTATION_ANTICLOCK );
350 setPermutation( ( face_permutation_type )face_permutation_type::ROTATION_TWICE_CLOCKWISE,
351 ( face_permutation_type )face_permutation_type::REVERSE_HEIGHT,
352 ( face_permutation_type )face_permutation_type::REVERSE_BASE );
354 setPermutation( ( face_permutation_type )face_permutation_type::PRINCIPAL_DIAGONAL,
355 ( face_permutation_type )face_permutation_type::REVERSE_HEIGHT,
356 ( face_permutation_type )face_permutation_type::ROTATION_ANTICLOCK );
358 setPermutation( ( face_permutation_type )face_permutation_type::ROTATION_CLOCKWISE,
359 ( face_permutation_type )face_permutation_type::ROTATION_ANTICLOCK,
360 ( face_permutation_type )face_permutation_type::ROTATION_TWICE_CLOCKWISE );
365 std::map<face_permutation_type, permutation_vector_type> vector_permutation;
366 std::map<face_permutation_type, permutation_matrix_type> matrix_permutation;