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
warpblend.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-03
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 */
30 #ifndef __WarpBlend_H
31 #define __WarpBlend_H 1
32 
33 #include <feel/feelalg/glas.hpp>
35 #include <feel/feelpoly/jacobi.hpp>
36 #include <feel/feelpoly/geomap.hpp>
37 #include <feel/feelmesh/geond.hpp>
38 #include <feel/feelmesh/geo0d.hpp>
40 
42 
43 
44 namespace Feel
45 {
46 
47 template< class Convex,
48  uint16_type Order,
49  typename T = double >
50 class PointSetWarpBlend : public PointSetInterpolation<Convex::nDim, Order, T, Simplex>
51 {
52 
53 public :
54 
55  typedef PointSetEquiSpaced<Convex, Order, T> super;
56 
57  typedef T value_type;
58 
59  static const uint32_type Dim = Convex::nDim;
60  static const uint16_type nPoints1D = Order+1;
61 
62  typedef typename super::return_type return_type;
63 
64  typedef ublas::vector<value_type> vector_type;
65 
66  static const uint32_type topological_dimension = Convex::topological_dimension;
67  static const uint32_type numPoints = super::numPoints;
68 
69  typedef Reference<Convex, Dim, Convex::nOrder, /*Dim*/Convex::nRealDim, value_type> reference_convex_type;
70 
71  typedef typename reference_convex_type::points_type points_type;
72 
73  reference_convex_type RefConv;
74 
75  PointSetWarpBlend( int interior = 0 )
76  {
77  PointSetEquiSpaced<Convex, Order, T> G( interior );
78 
79  points_type final_pts = G.points();
80 
81  //Copies information about EquiSpaced points
82  this->setEid( G.getEid() );
83  this->setPtE( G.getPtE() );
84 
85  if ( Dim == 1 )
86  {
87  PointSetGaussLobatto<Hypercube<1,1>, Order, value_type> gausslobatto( interior );
88 
89  final_pts = gausslobatto.points();
90  }
91 
92  else if ( Order > 2 )
93  {
94  entities = std::make_pair( RefConv.vertices(), equiVertices() );
95 
96  if ( Dim == 2 )
97  final_pts = transformPoints<2>( final_pts );
98 
99  else
100  {
101  uint16_type max_entity_dim = 1;
102 
103  if ( Order > 3 )
104  max_entity_dim = 2;
105 
106  blendAxis();
107 
108  for ( uint16_type d = 1; d < max_entity_dim +1 ; ++d )
109  {
110  for ( int e = RefConv.entityRange( d ).begin();
111  e < RefConv.entityRange( d ).end();
112  ++e )
113  {
114  points_type pts = toEquilateral( G.pointsBySubEntity( d, e, 0 ), true );
115 
116  points_type coord_bar = toBarycentric( pts );
117 
118  pts += calculateFaceDeformation( coord_bar, entityMap( d,e ), mpl::int_<3>() );
119 
120  final_pts = putInPointset ( final_pts,
121  toEquilateral( pts, false ),
122  G.interiorRangeById( d, e ) );
123  }
124  }
125 
126  if ( Order > 3 )
127  {
128  final_pts = putInPointset ( final_pts,
129  transformPoints<3>( G.pointsBySubEntity( 3, 0, 0 ) ),
130  G.interiorRangeById( 3, 0 ) );
131  }
132  }
133  }
134 
135  this->setPoints( final_pts );
136 
137  this->setName( "warpblend", Order );
138  }
139 
140  ~PointSetWarpBlend() {}
141 
142 private :
143 
144  std::pair<points_type, points_type> entities;
145 
146  std::map<uint16_type, std::pair<vector_type, vector_type > > axis;
147 
148  points_type equiVertices ()
149  {
150  points_type V ( ublas::scalar_matrix<value_type>( Dim, Dim+1, value_type( 0 ) ) );
151 
152  for ( uint16_type i=0; i < 3; i++ )
153  {
154  value_type angle = M_PI*( value_type( 7 ) + value_type( 4 )*value_type( i ) )/value_type( 6 );
155 
156  V( 0,i ) = math::cos( angle );
157  V( 1,i ) = math::sin( angle );
158 
159  if ( Dim == 3 )
160  V( 2,i ) = - math::sqrt( value_type( 2 ) )/value_type( 4 );
161  }
162 
163  V *= value_type( 2 )/math::sqrt( value_type( 3 ) );
164 
165  if ( Dim == 3 )
166  V( 2,3 ) = math::sqrt( value_type( 6 ) )/value_type( 2 );
167 
168  return V;
169  }
170 
171  void blendAxis()
172  {
173  std::vector<uint16_type> indices;
174 
175  for ( uint16_type i=0; i<4; i++ )
176  {
177  for ( uint16_type j=0; j<4; j++ )
178  {
179  if ( j != ( 18 +i*( -6 + ( i-1 )*( -3 + 4*( i-2 ) ) ) )/6 )
180  indices.push_back( j );
181  }
182 
183  axis[i].first = getVertex( 1,indices[1] ) - getVertex( 1,indices[0] );
184  axis[i].second = getVertex( 1,indices[2] ) - ( getVertex( 1,indices[1] ) + getVertex( 1,indices[0] ) )/value_type( 2 );
185 
186  axis[i].first = axis[i].first/ublas::norm_2( axis[i].first );
187  axis[i].second = axis[i].second/ublas::norm_2( axis[i].second );
188 
189  indices.clear();
190  }
191  }
192 
193  //correspondence between the faces and edges of the reference and the equilateral simplex used
194  uint16_type entityMap ( uint16_type top_dim, uint16_type id )
195  {
196  if ( top_dim == 2 )
197  id = ( 12 + id*( 6 + ( id-1 )*( -9 + 4*( id-2 ) ) ) )/6;
198 
199  else
200  {
201  if ( id == 5 )
202  id = 2;
203 
204  else if ( id > 2 )
205  id = 1;
206 
207  else
208  id = 0;
209  }
210 
211  return id;
212  }
213 
214  vector_type getVertex ( uint16_type element, uint16_type i )
215  {
216  vector_type p;
217 
218  if ( element == 0 )
219  p = ublas::column( entities.first, i );
220 
221  else
222  p = ublas::column( entities.second, i );
223 
224  return p;
225  }
226 
227  points_type toEquilateral ( points_type pts, bool toEqui )
228  {
229  points_type coord_ref_elem ( Dim, Dim );
230  points_type coord_equi_elem ( Dim, Dim );
231 
232  for ( uint16_type i = 0; i < Dim; i++ )
233  {
234  ublas::column( coord_ref_elem, i ) = getVertex( 0,Dim ) - getVertex( 0,i );
235  ublas::column( coord_equi_elem, i ) = getVertex( 1,Dim ) - getVertex( 1,i );
236  }
237 
238  points_type A ( Dim, Dim );
239  points_type b ( Dim, 1 );
240 
241  LU< points_type > lu( coord_ref_elem );
242  lu.inverse( A );
243 
244  A = ublas::prod( coord_equi_elem , A );
245 
246  ublas::column( b,0 ) = getVertex( 1,0 ) - ublas::prod( A, getVertex( 0,0 ) );
247 
248  if ( toEqui )
249  {
250  pts = ublas::prod( A, pts );
251 
252  for ( uint16_type i=0; i < pts.size2(); i++ )
253  ublas::column( pts, i ) += ublas::column( b,0 );
254  }
255 
256  else
257  {
258  for ( uint16_type i=0; i < pts.size2(); i++ )
259  ublas::column( pts, i ) -= ublas::column( b,0 );
260 
261  LU< points_type > lu( A );
262  pts = lu.solve( pts ) ;
263  }
264 
265  return pts;
266  }
267 
268  points_type toBarycentric ( points_type pts )
269  {
270  points_type C( Dim, Dim );
271 
272  for ( uint16_type i = 0; i < Dim; i++ )
273  ublas::column( C, i ) = getVertex( 1, ( Dim - Dim%2 + i )%Dim ) - getVertex( 1,Dim );
274 
275  for ( uint16_type i=0; i < pts.size2(); i++ )
276  ublas::column( pts, i ) -= getVertex( 1,Dim );
277 
278  LU< points_type > lu( C );
279 
280  points_type solution = lu.solve( pts );
281 
282  pts.resize( entities.second.size2(), pts.size2() );
283 
284  ublas::subrange( pts, 1, Dim+1, 0, pts.size2() ) = solution;
285 
286  ublas::row( pts, 0 ) = ublas::scalar_vector<value_type>( pts.size2(), value_type( 1 ) );
287 
288  for ( uint16_type i = 1; i < Dim+1; i++ )
289  ublas::row( pts, 0 ) -= ublas::row( pts, i );
290 
291  return pts;
292  }
293 
294  vector_type gll()
295  {
296  vector_type gx( Order+1 );
297  vector_type gw( Order+1 );
298 
299  details::dyna::gausslobattojacobi<value_type, vector_type, vector_type>( Order+1, gw, gx );
300 
301  return gx;
302  }
303 
304  value_type alpha( uint16_type d ) const
305  {
306  value_type p;
307 
308  double __alpha_2d[ 16 ] = { 0,
309  0,
310  0,
311  1.4152,
312  0.1001,
313  0.2751,
314  0.9808,
315  1.0999,
316  1.2832,
317  1.3648,
318  1.4773,
319  1.4959,
320  1.5743,
321  1.5770,
322  1.6223,
323  1.6258
324  };
325 
326  double __alpha_3d[ 16 ] = { 0,
327  0,
328  0,
329  0,
330  0.1002,
331  1.1332,
332  1.5608,
333  1.3413,
334  1.2577,
335  1.1603,
336  1.0153,
337  0.6080,
338  0.4523,
339  0.8856,
340  0.8717,
341  0.9655
342  };
343 
344  if ( d == 2 )
345  p = ( Order<=15 )?__alpha_2d[ Order ]:( 5./3. );
346 
347  else
348  p = ( Order<=15 )?__alpha_3d[ Order ]:1.0;
349 
350  return p;
351  }
352 
353  template<typename AE>
354  vector_type warpFactor( vector_type const& x, ublas::vector_expression<AE> const& xout ) const
355  {
356  vector_type warp( xout().size() );
357  warp = ublas::scalar_vector<value_type>( xout().size(), value_type( 0 ) );
358  vector_type xeq( glas::linspace( value_type( -1 ), value_type( 1 ), nPoints1D, 0 ) );
359 
360  vector_type d( xout().size() );
361 
362  for ( int i = 0; i < nPoints1D; ++i )
363  {
364  d = ublas::scalar_vector<value_type>( xout().size(), x( i )-xeq( i ) );
365 
366  for ( int j = 1; j < nPoints1D-1; ++j )
367  if ( i != j )
368  {
369  d = ( ublas::element_prod( d, xout ) - d *xeq( j ) )/( xeq( i ) - xeq( j ) );
370  }
371 
372  // deflate end roots
373  if ( i != 0 )
374  d *= value_type( -1 )/( xeq( i ) - xeq( 0 ) );
375 
376  if ( i != nPoints1D-1 )
377  d *= value_type( 1 )/( xeq( i ) - xeq( nPoints1D-1 ) );
378 
379  warp += d;
380  }
381 
382  return warp;
383  }
384 
385  //retrieves the right barycentric coordinates for each warp function
386  points_type getCoordinates( points_type const& pts, uint16_type face_id )
387  {
388  std::map<uint16_type, std::vector<uint16_type> > faces;
389 
390  for ( uint16_type i=0; i<4; i++ )
391  for ( uint16_type j=0; j<4; j++ )
392  {
393  if ( j != i )
394  faces[i].push_back( j );
395  }
396 
397  for ( uint16_type i=2; i<4; i++ )
398  {
399  faces[i][1] = faces[i][2];
400  faces[i][2] = 1;
401  }
402 
403  points_type coord ( 3, pts.size2() );
404 
405  for ( uint16_type i=0; i<3; i++ )
406  ublas::row( coord, i ) = ublas::row( pts, faces[face_id][i] );
407 
408  return coord;
409  }
410 
411  points_type calculateFaceDeformation( points_type const& coord_bar, mpl::int_<2> )
412  {
413  points_type blend ( 3, coord_bar.size2() );
414 
415  for ( uint16_type i = 0; i < 3; i++ )
416  ublas::row( blend, i ) = ublas::element_prod( ublas::row( coord_bar, ( i+1 )%3 ), ublas::row( coord_bar, ( i+2 )%3 ) );
417 
418 
419  blend *= value_type( 4 );
420 
421  points_type warpfactor( 3, coord_bar.size2() );
422 
423  for ( uint16_type i = 0; i < 3; i++ )
424  ublas::row( warpfactor, i ) = warpFactor( gll(), ublas::row( coord_bar, ( i+2 )%3 ) - ublas::row( coord_bar, ( i+1 )%3 ) );
425 
426  points_type warp( 3, coord_bar.size2() );
427 
428  warp = ublas::element_prod( ublas::element_prod( blend, warpfactor ),
429  ublas::scalar_matrix<value_type>( 3, coord_bar.size2(), value_type( 1 ) ) +
430  alpha( 2 )*alpha( 2 )*ublas::element_prod( coord_bar, coord_bar ) );
431 
432  points_type deformation( ublas::scalar_matrix<value_type>( 2, warp.size2(), value_type( 0 ) ) );
433 
434  for ( uint16_type i=0; i < warp.size1(); i++ )
435  {
436  value_type angle = value_type( 2*i )*M_PI/value_type( 3 );
437 
438  ublas::row( deformation, 0 ) += math::cos( angle )*ublas::row( warp, i );
439  ublas::row( deformation, 1 ) += math::sin( angle )*ublas::row( warp, i );
440  }
441 
442  return deformation;
443  }
444 
445  //calculates the face deformation for one of the faces of the tetrahedra
446  points_type calculateFaceDeformation( points_type const& pts, uint16_type face_id, mpl::int_<3> )
447  {
448  points_type coord_bar;
449 
450  coord_bar = getCoordinates( pts, face_id );
451 
452  points_type def = calculateFaceDeformation( coord_bar, mpl::int_<2>() );
453 
454  points_type w ( 3, pts.size2() );
455 
456  for ( uint16_type i=0; i<3; i++ )
457  ublas::row( w, i ) = axis[face_id].first( i )*ublas::row( def,0 ) + axis[face_id].second( i )*ublas::row( def,1 );
458 
459  return w;
460  }
461 
462  points_type blendFunction( points_type const& pts )
463  {
464  points_type blend ( 4, pts.size2() );
465 
466  for ( uint16_type i=0; i<4; i++ )
467  {
468  vector_type aux1 ( ublas::scalar_vector<value_type>( pts.size2(), value_type( 1 ) ) );
469  vector_type aux2 = aux1;
470 
471  for ( uint16_type j=0; j<4; j++ )
472  {
473  if ( j != i )
474  {
475  aux1 = ublas::element_prod( aux1, ublas::row( pts, j ) );
476  aux2 = ublas::element_prod( aux2, value_type( 2 )*ublas::row( pts, j ) + ublas::row( pts, i ) );
477  }
478  }
479 
480  ublas::row( blend, i ) = ublas::element_div( aux1, aux2 );
481  }
482 
483  blend *= value_type( 8 );
484 
485  //add beta deformation
486  blend = ublas::element_prod( blend, ublas::scalar_matrix<value_type>( 4, pts.size2(), value_type( 1 ) ) +
487  alpha( 3 )*alpha( 3 )*ublas::element_prod( pts, pts ) );
488 
489  return blend;
490  }
491 
492  points_type calculateFaceDeformation( points_type const& coord_bar, mpl::int_<3> )
493  {
494  points_type blend = blendFunction( coord_bar );
495 
496  std::map<uint16_type, points_type > warp;
497 
498  for ( uint16_type face_id = 0; face_id < 4; face_id++ )
499  {
500  warp[face_id].resize( 3, coord_bar.size2() );
501 
502  for ( uint16_type i=0; i<3; i++ )
503  {
504  vector_type aux = ublas::row( calculateFaceDeformation( coord_bar,
505  entityMap( 2, face_id ),
506  mpl::int_<3>() ), i );
507 
508  ublas::row( warp[face_id], i ) = ublas::element_prod( ublas::row( blend, face_id ), aux );
509  }
510  }
511 
512  for ( uint16_type face_id = 1; face_id < 4; face_id++ )
513  warp[0] += warp[face_id];
514 
515  return warp[0];
516  }
517 
518  template<int N>
519  points_type transformPoints( points_type pts )
520  {
521  pts = toEquilateral( pts, true );
522 
523  points_type coord_bar = toBarycentric( pts );
524 
525  pts += calculateFaceDeformation( coord_bar, mpl::int_<N>() );
526 
527  return toEquilateral( pts, false );
528  }
529 
530  points_type putInPointset ( points_type final_pts, points_type pts, std::pair<uint16_type, uint16_type> position )
531  {
532  ublas::subrange( final_pts, 0, 3, position.first, position.second ) = pts;
533 
534  return final_pts;
535  }
536 };
537 
538 } // Feel
539 #endif /* __WarpBlend_H */

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