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
quadmapped.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 */
29 #ifndef __Quadmapped_H
30 #define __Quadmapped_H 1
31 
34 
36 
37 #include <stdexcept>
38 #include <boost/numeric/ublas/matrix.hpp>
39 #include <boost/numeric/ublas/io.hpp>
40 
41 #include <feel/feelcore/traits.hpp>
42 #include <feel/feelalg/glas.hpp>
43 
44 #include <feel/feelmesh/geoelement.hpp>
45 #include <feel/feelpoly/im.hpp>
46 
47 
48 
49 namespace Feel
50 {
51 namespace ublas = boost::numeric::ublas;
52 
53 /*template< class Convex,
54  uint16_type Order,
55  typename T = double ,
56  template<class, uint16_type, class> class PointSetType = Gauss>*/
57 
58 template<typename PsetType>
59 class QuadMapped : public PsetType
60 {
61  typedef PsetType super;
62 public :
63  typedef super pointset_type;
64 
65  typedef typename super::value_type value_type;
66 
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;
70 
71  static const bool is_simplex = convex_type::is_simplex;
72 
73  typedef typename pointset_type::nodes_type nodes_type;
74  typedef typename matrix_node<value_type>::type points_type;
75 
76 
77 
78  typedef ublas::vector<uint16_type> permutation_vector_type;
79  typedef ublas::mapped_matrix<uint16_type> permutation_matrix_type;
80 
81  //typedef mpl::if_< mpl::bool_< is_simplex >, Simplex<Dim, Order, Dim> , Hypercube<Dim, Order, Dim> > conv_order_type;
82  typedef convex_type conv_order_type;
83 
84 #if 0
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;
88 #else
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;
92 #endif
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;
98 
99  QuadMapped()
100  :
101  super()
102  {
103  }
104 
105  ~QuadMapped() {}
106 
107  permutation_points_type
108  operator()( pointset_type const& /*pts*/ )
109  {
110  permutation_points_type ppts( convex_type::numTopologicalFaces );
111 
112  // loop over all permutations on the reference convex and
113  // store all permutations of the pointset
114  for ( int f = 0; f < convex_type::numTopologicalFaces; ++f )
115  {
116  __typeof__( typename pointset_type::Face( *this, f ) ) facequad = typename pointset_type::Face( *this, f );
117 
118  generate( ppts, facequad, f, mpl::int_<Dim>() );
119  }
120 
121  return ppts;
122  }
123 
124 private:
125 
126  permutation_vector_type getVectorPermutation ( face_permutation_type P )
127  {
128  return vector_permutation[P];
129  }
130 
131  permutation_matrix_type getMatrixPermutation ( face_permutation_type P )
132  {
133  return matrix_permutation[P];
134  }
135 
136  void
137  generate( permutation_points_type& ppts, typename pointset_type::Face const& pts, int f, mpl::int_<1> )
138  {
139  vertex_permutation_type pbegin( vertex_permutation_type::IDENTITY );
140  vertex_permutation_type pend( vertex_permutation_type::N_PERMUTATIONS );
141 
142  for ( vertex_permutation_type p = pbegin; p < pend; ++p )
143  {
144  ppts[f].insert( std::make_pair( p, pts.points() ) );
145  DVLOG(2) << "[quadmapped] ppts[" << f << "]["
146  << p.value() << "]="
147  << ppts[f].find( p )->second << "\n";
148  }
149  }
150  void
151  generate( permutation_points_type& ppts, typename pointset_type::Face const& pts, int f, mpl::int_<2> )
152  {
153  edge_permutation_type pbegin( edge_permutation_type::IDENTITY );
154  edge_permutation_type pend( edge_permutation_type::N_PERMUTATIONS );
155 
156  for ( edge_permutation_type p = pbegin; p < pend; ++p )
157  {
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";
160  }
161  }
162  void
163  generate( permutation_points_type& ppts, typename pointset_type::Face const& pts, int f, mpl::int_<3> )
164  {
165  uint16_type npts = pts.nPoints();
166 
167  // generate all permutations
168  generateFacePermutations( npts, mpl::bool_<is_simplex>() );
169  //generateFacePermutations( npts, mpl::bool_<false>() );
170  face_permutation_type pbegin( face_permutation_type::IDENTITY );
171  face_permutation_type pend( face_permutation_type::N_PERMUTATIONS );
172 
173  for ( face_permutation_type p = pbegin; p < pend; ++p )
174  {
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";
177  }
178  }
179 
180  // Permutates points in the reference element that, when mapped to the
181  // real one, match others generated in the same real edge
182  points_type permutatePoints ( points_type Gt,
183  edge_permutation_type edge_perm )
184  {
185  if ( edge_perm != edge_permutation_type( edge_permutation_type::IDENTITY ) )
186  {
187  // this formula works with 1-point quadratures
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 ) );
190  }
191 
192  return Gt;
193  }
194 
195  // Permutates points in the reference element that, when mapped to the
196  // real one, match others generated in the same real face
197  points_type permutatePoints ( nodes_type const& Gt,
198  face_permutation_type face_perm )
199  {
200  FEELPP_ASSERT( face_perm != face_permutation_type( face_permutation_type::NO_PERMUTATION ) )
201  ( face_perm ).error( "invalid permutation" );
202 
203  nodes_type res( Gt );
204 
205  if ( face_perm != face_permutation_type( face_permutation_type::IDENTITY ) )
206  res = prod( Gt, getMatrixPermutation( face_perm ) );
207 
208  FEELPP_ASSERT( res.size1() == Gt.size1() &&
209  res.size2() == Gt.size2() ) ( res )( Gt ).error ( "invalid permutation operation" );
210  return res;
211  }
212 
213  permutation_matrix_type vectorToMatrixPermutation ( permutation_vector_type const& v )
214  {
215  permutation_matrix_type P ( v.size(), v.size(), v.size()*v.size() );
216  P.clear();
217 
218  for ( uint16_type i = 0; i < v.size(); ++i )
219  P( i,v( i ) ) = 1;
220 
221  return P;
222  }
223 
224  permutation_vector_type matrixToVectorPermutation ( permutation_matrix_type const& P )
225  {
226  FEELPP_ASSERT( P.size1() == P.size2() ).error( "invalid permutation" );
227 
228  permutation_vector_type v ( P.size1() );
229  v.clear();
230 
231  for ( uint16_type i = 0; i < v.size(); ++i )
232  for ( uint16_type j = 0; j < v.size(); ++j )
233  if ( P( i,j ) == 1 )
234  v( i ) = j;
235 
236  return v;
237  }
238 
239  void setPermutation( face_permutation_type out, face_permutation_type first, face_permutation_type second )
240  {
241  matrix_permutation[out] = ublas::prod( matrix_permutation[first],
242  matrix_permutation[second] );
243 
244  vector_permutation[out] = matrixToVectorPermutation( matrix_permutation[out] );
245  }
246 
247 
248  //Permutations for tetrahedra
249  void generateFacePermutations( uint16_type npoints, mpl::bool_<true> )
250  {
251  //uint16_type npoints = n_side_points*(n_side_points+1)/2;
252  //uint16_type n_side_points = (uint16_type) (-1 + math::sqrt( (double)8*npoints + 1 ))/2;
253  uint16_type n_side_points = ( uint16_type )math::sqrt( ( double )npoints );
254  //FEELPP_ASSERT( npoints == (n_side_points-1)*(n_side_points-2)/2 )
255  //FEELPP_ASSERT( npoints == (n_side_points)*(n_side_points))
256  //( npoints )( n_side_points ).error( "invalid number of points" );
257  permutation_vector_type _vec( npoints );
258  _vec.clear();
259 
260  //define hypotenuse reverse permutation
261  for ( uint16_type i = 0; i < n_side_points; ++i )
262  {
263  for ( uint16_type j = 0; j <= i; j++ )
264  {
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 ) );
267 
268  _vec( _first ) = _last;
269  }
270  }
271 
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 );
274 
275  //define base reverse permutation
276  for ( uint16_type i = 0; i < n_side_points; ++i )
277  {
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;
280 
281  for ( uint16_type j = 0; j <= n_side_points - i - 1; j++ )
282  _vec( _begin + j ) = _end - j;
283  }
284 
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 );
287 
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 );
291 
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 );
295 
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 );
299  }
300 
301  //Permutations for hexahedra
302  void generateFacePermutations( uint16_type npoints, mpl::bool_<false> )
303  {
304  //uint16_type npoints = n_side_points*(n_side_points);
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 );
309  _vec.clear();
310  //define base permutation (tau_2)
311  uint16_type p=0;
312 
313  for ( int16_type i = n_side_points-1; i >= 0; --i )
314  {
315  uint16_type _first = i*n_side_points;
316 
317  for ( uint16_type j = 0; j <= n_side_points-1; j++ )
318  {
319  _vec( p ) = _first + j;
320  p++;
321  }
322  }
323 
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 );
326 
327  //define once anticlockwise permutation (tau_3)
328  p=0;
329 
330  for ( int16_type i = n_side_points-1; i >= 0; --i )
331  {
332  for ( uint16_type j = 0; j <= n_side_points-1; j++ )
333  {
334  _vec( p ) = i + n_side_points*j;
335  p++;
336  }
337  }
338 
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 );
341 
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 );
345 
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 );
349 
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 );
353 
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 );
357 
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 );
361  }
362 
363 private:
364 
365  std::map<face_permutation_type, permutation_vector_type> vector_permutation;
366  std::map<face_permutation_type, permutation_matrix_type> matrix_permutation;
367 
368 
369 
370 };
371 
372 } // Feel
373 #endif /* __QuadMapped_H */

Generated on Sun Dec 22 2013 13:11:11 for Feel++ by doxygen 1.8.5