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
mapped.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 __PointSetMapped_H
30 #define __PointSetMapped_H 1
31 
32 #include <stdexcept>
33 
34 #include <boost/numeric/ublas/matrix_sparse.hpp>
35 #include <boost/numeric/ublas/io.hpp>
36 
37 
38 #include <feel/feelcore/feel.hpp>
40 #include <feel/feelcore/traits.hpp>
41 
42 #include <feel/feelalg/glas.hpp>
43 
45 #include <feel/feelmesh/geoelement.hpp>
46 
48 
49 
50 
51 
52 
53 namespace Feel
54 {
55 namespace ublas = boost::numeric::ublas;
56 
57 template< typename element_type,
58  class Convex,
59  uint16_type Order,
60  typename T = double ,
61  template<class, uint16_type, class> class PointSetType = PointSetEquiSpaced>
62 class PointSetMapped : public PointSetType<Convex, Order, T>
63 {
64 
65 public :
66  typedef PointSetType<Convex, Order, T> pointset_type;
67 
68  typedef T value_type;
69 
70  static const uint32_type Dim = Convex::nDim;
71  static const uint32_type convexOrder = Convex::nOrder;
72 
73  static const bool is_simplex = Convex::is_simplex;
74 
75  typedef typename pointset_type::nodes_type nodes_type;
76  typedef typename matrix_node<value_type>::type points_type;
77 
78  typedef typename element_type::gm_type gm_type;
79  typedef typename element_type::gm_ptrtype gm_ptrtype;
80 
81  typedef typename element_type::edge_permutation_type edge_permutation_type;
82  typedef typename element_type::face_permutation_type face_permutation_type;
83 
84  typedef ublas::vector<uint16_type> permutation_vector_type;
85  typedef ublas::mapped_matrix<uint16_type> permutation_matrix_type;
86 
87  typedef mpl::if_< mpl::bool_< is_simplex >,
88  Simplex<Dim, Order, Dim> ,
89  Hypercube<Dim, Order, Dim> > conv_order_type;
90 
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;
94 
95  typedef Reference<Convex, Dim, convexOrder, Dim, value_type> RefElem;
96 
97  typedef typename pointset_type::range_type range_type;
98  typedef typename pointset_type::index_map_type index_map_type;
99 
100  RefElem RefConv;
101 
102  PointSetMapped( element_type const& _elt )
103  :
104  M_elt( _elt )
105  {
106  pointset_type pts;
107 
108  points_type Gt = updatePoints<2>( updatePoints<1>( pts.points(), pts, mpl::bool_< ( Order > 2 ) >() ),
109  pts,
110  mpl::bool_< ( Dim == 3 && Order > 2+uint16_type( is_simplex ) ) >() );
111 
112  this->setPoints( Gt );
113  }
114 
115  ~PointSetMapped() {}
116 
117  permutation_vector_type getVectorPermutation ( face_permutation_type P )
118  {
119  return vector_permutation[P];
120  }
121 
122  permutation_matrix_type getMatrixPermutation ( face_permutation_type P )
123  {
124  return matrix_permutation[P];
125  }
126 
127  points_type pointsBySubEntity( uint16_type top_dim, uint16_type local_id, bool boundary = 0, bool real = 0 )
128  {
129  index_map_type index_list = this->entityToLocal( top_dim, local_id, boundary );
130 
131  //Rearrange the order of the vertices, edges and faces
132 
133  uint16_type matrix_size = 0;
134 
135  if ( index_list[0].size() != 0 )
136  matrix_size +=index_list[0].size()*nbPtsPerVertex;
137 
138  if ( ( top_dim>=1 ) && ( index_list[1].size() != 0 ) )
139  matrix_size +=index_list[1].size()*nbPtsPerEdge;
140 
141  if ( ( top_dim>=2 ) && ( index_list[2].size() != 0 ) )
142  matrix_size +=index_list[2].size()*nbPtsPerFace;
143 
144  points_type G ( Dim, matrix_size );
145 
146  for ( uint16_type i=0, p=0; i < top_dim+1; i++ )
147  {
148  if ( index_list[i].size() )
149  {
150  for ( uint16_type j=0; j < index_list[i].size(); j++ )
151  {
152  //Something else needs to be done here: depending on the entity, and in the case
153  //we have boundary=1, the points in the boundary have to be reordered
154 
155  points_type aux = this->interiorPointsById( i, index_list[i][j] );
156 
157  ublas::subrange( G, 0, Dim, p, p+aux.size2() ) = aux;
158 
159  p+=aux.size2();
160  }
161  }
162  }
163 
164  if ( real )
165  G = transformToReal( G );
166 
167  return G;
168  }
169 
170 private:
171 
172  element_type M_elt;
173 
174  std::map<face_permutation_type, permutation_vector_type> vector_permutation;
175  std::map<face_permutation_type, permutation_matrix_type> matrix_permutation;
176 
177  template <int top_dim>
178  nodes_type updatePoints( nodes_type const& pts, pointset_type const& /*G*/, mpl::bool_<false> )
179  {
180  return pts;
181  }
182 
183  template <int top_dim>
184  nodes_type updatePoints( nodes_type const& thepts, pointset_type const& G, mpl::bool_<true> )
185  {
186  nodes_type pts( thepts );
187 
188  if ( top_dim == 2 )
189  {
190  generateFacePermutations( Order - 1-uint16_type( is_simplex ),
191  mpl::int_<Dim>(), mpl::bool_< is_simplex >() );
192  }
193 
194  for ( size_type e = RefConv.entityRange( top_dim ).begin();
195  e < RefConv.entityRange( top_dim ).end();
196  ++e )
197  {
198  points_type Gt = G.pointsBySubEntity( top_dim, e );
199 
200  Gt = permutatePoints ( Gt, e, mpl::int_<top_dim>() );
201 
202  uint16_type pos = G.interiorRangeById( top_dim, e ).first;
203 
204  ublas::subrange( pts, 0, Dim, pos, pos+Gt.size2() ) = Gt;
205  }
206 
207  return pts;
208  }
209 
210  // Permutates points in the reference element that, when mapped to the
211  // real one, match others generated in the same real edge
212  points_type permutatePoints ( nodes_type const& theGt, size_type entity_local_id, mpl::int_<1> )
213  {
214  nodes_type Gt( theGt );
215  edge_permutation_type permutation = M_elt.edgePermutation( entity_local_id );
216 
217  if ( permutation != edge_permutation_type( edge_permutation_type::IDENTITY ) )
218  {
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 ) );
221  }
222 
223  return Gt;
224  }
225 
226  // Permutates points in the reference element that, when mapped to the
227  // real one, match others generated in the same real face
228  points_type permutatePoints ( nodes_type const& Gt, size_type entity_local_id, mpl::int_<2> )
229  {
230  face_permutation_type permutation = M_elt.facePermutation( entity_local_id );
231  nodes_type res( Gt );
232 
233  if ( permutation != face_permutation_type( face_permutation_type::IDENTITY ) )
234  res = prod( Gt, getMatrixPermutation( permutation ) );
235 
236  return res;
237  }
238 
239  permutation_matrix_type vectorToMatrixPermutation ( permutation_vector_type const& v )
240  {
241  permutation_matrix_type P ( v.size(), v.size(), v.size()*v.size() );
242 
243  for ( uint16_type i = 0; i < v.size(); ++i )
244  P( i,v( i ) ) = 1;
245 
246  return P;
247  }
248 
249  permutation_vector_type matrixToVectorPermutation ( permutation_matrix_type const& P )
250  {
251  FEELPP_ASSERT( P.size1() == P.size2() ).error( "invalid permutation" );
252 
253  permutation_vector_type v ( P.size1() );
254 
255  for ( uint16_type i = 0; i < v.size(); ++i )
256  for ( uint16_type j = 0; j < v.size(); ++j )
257  if ( P( i,j ) == 1 )
258  v( i ) = j;
259 
260  return v;
261  }
262 
263  void setPermutation( int out, int first, int second )
264  {
265  matrix_permutation[face_permutation_type( out )] = ublas::prod( matrix_permutation[face_permutation_type( first )],
266  matrix_permutation[face_permutation_type( second )] );
267 
268  vector_permutation[face_permutation_type( out )] = matrixToVectorPermutation( matrix_permutation[face_permutation_type( out )] );
269  }
270 
271  void generateFacePermutations( uint16_type /*n_side_points*/, mpl::int_<2>, mpl::bool_<true> ) {}
272  void generateFacePermutations( uint16_type /*n_side_points*/, mpl::int_<2>, mpl::bool_<false> ) {}
273 
274  //Permutations for tetrahedra
275  void generateFacePermutations( uint16_type n_side_points, mpl::int_<3>, mpl::bool_<true> )
276  {
277  uint16_type npoints = n_side_points*( n_side_points+1 )/2;
278 
279  permutation_vector_type _vec( npoints );
280 
281  //define hypotenuse reverse permutation
282  for ( uint16_type i = 0; i < n_side_points; ++i )
283  {
284  for ( uint16_type j = 0; j <= i; j++ )
285  {
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 ) );
288 
289  _vec( _first ) = _last;
290  }
291  }
292 
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 );
295 
296  //define base reverse permutation
297  for ( uint16_type i = 0; i < n_side_points; ++i )
298  {
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;
301 
302  for ( uint16_type j = 0; j <= n_side_points - i - 1; j++ )
303  _vec( _begin + j ) = _end - j;
304  }
305 
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 );
308 
309  setPermutation( face_permutation_type::ROTATION_ANTICLOCK,
310  face_permutation_type::REVERSE_BASE,
311  face_permutation_type::REVERSE_HYPOTENUSE );
312 
313  setPermutation( face_permutation_type::ROTATION_CLOCKWISE,
314  face_permutation_type::REVERSE_HYPOTENUSE,
315  face_permutation_type::REVERSE_BASE );
316 
317  setPermutation( face_permutation_type::REVERSE_HEIGHT,
318  face_permutation_type::REVERSE_BASE,
319  face_permutation_type::ROTATION_CLOCKWISE );
320  }
321 
322  //Permutations for hexahedra
323  void generateFacePermutations( uint16_type n_side_points, mpl::int_<3>, mpl::bool_<false> )
324  {
325  uint16_type npoints = n_side_points*( n_side_points );
326 
327  permutation_vector_type _vec( npoints );
328 
329  //define base permutation (tau_2)
330  uint16_type p=0;
331 
332  for ( int16_type i = n_side_points-1; i >= 0; --i )
333  {
334  uint16_type _first = i*n_side_points;
335 
336  for ( uint16_type j = 0; j <= n_side_points-1; j++ )
337  {
338  _vec( p ) = _first + j;
339  p++;
340  }
341  }
342 
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 );
345 
346  //define once anticlockwise permutation (tau_3)
347  p=0;
348 
349  for ( int16_type i = n_side_points-1; i >= 0; --i )
350  {
351  for ( uint16_type j = 0; j <= n_side_points-1; j++ )
352  {
353  _vec( p ) = i + n_side_points*j;
354  p++;
355  }
356  }
357 
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 );
360 
361  setPermutation( face_permutation_type::SECOND_DIAGONAL,
362  face_permutation_type::REVERSE_BASE,
363  face_permutation_type::ROTATION_ANTICLOCK );
364 
365  setPermutation( face_permutation_type::REVERSE_HEIGHT,
366  face_permutation_type::SECOND_DIAGONAL,
367  face_permutation_type::ROTATION_ANTICLOCK );
368 
369  setPermutation( face_permutation_type::ROTATION_TWICE_CLOCKWISE,
370  face_permutation_type::REVERSE_HEIGHT,
371  face_permutation_type::REVERSE_BASE );
372 
373  setPermutation( face_permutation_type::PRINCIPAL_DIAGONAL,
374  face_permutation_type::REVERSE_HEIGHT,
375  face_permutation_type::ROTATION_ANTICLOCK );
376 
377  setPermutation( face_permutation_type::ROTATION_CLOCKWISE,
378  face_permutation_type::ROTATION_ANTICLOCK,
379  face_permutation_type::ROTATION_TWICE_CLOCKWISE );
380  }
381 
382  points_type transformToReal( points_type const& Gt )
383  {
384  gm_ptrtype _gm_ptr( new gm_type );
385 
386  typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( _gm_ptr, Gt ) );
387 
388  typename gm_type::template Context<vm::POINT, element_type> gmc( _gm_ptr, M_elt, __geopc );
389 
390  return gmc.xReal();
391  }
392 };
393 
394 } // Feel
395 #endif /* __PointSetMapped_H */

Generated on Sun Oct 20 2013 08:25:01 for Feel++ by doxygen 1.8.4