31 #define __WarpBlend_H 1
36 #include <feel/feelpoly/geomap.hpp>
38 #include <feel/feelmesh/geo0d.hpp>
47 template<
class Convex,
50 class PointSetWarpBlend :
public PointSetInterpolation<Convex::nDim, Order, T, Simplex>
55 typedef PointSetEquiSpaced<Convex, Order, T> super;
59 static const uint32_type Dim = Convex::nDim;
60 static const uint16_type nPoints1D = Order+1;
62 typedef typename super::return_type return_type;
64 typedef ublas::vector<value_type> vector_type;
66 static const uint32_type topological_dimension = Convex::topological_dimension;
67 static const uint32_type numPoints = super::numPoints;
69 typedef Reference<Convex, Dim, Convex::nOrder, Convex::nRealDim, value_type> reference_convex_type;
71 typedef typename reference_convex_type::points_type points_type;
73 reference_convex_type RefConv;
75 PointSetWarpBlend(
int interior = 0 )
77 PointSetEquiSpaced<Convex, Order, T> G( interior );
79 points_type final_pts = G.points();
82 this->setEid( G.getEid() );
83 this->setPtE( G.getPtE() );
87 PointSetGaussLobatto<Hypercube<1,1>, Order, value_type> gausslobatto( interior );
89 final_pts = gausslobatto.points();
94 entities = std::make_pair( RefConv.vertices(), equiVertices() );
97 final_pts = transformPoints<2>( final_pts );
101 uint16_type max_entity_dim = 1;
108 for ( uint16_type d = 1; d < max_entity_dim +1 ; ++d )
110 for (
int e = RefConv.entityRange( d ).begin();
111 e < RefConv.entityRange( d ).end();
114 points_type pts = toEquilateral( G.pointsBySubEntity( d, e, 0 ), true );
116 points_type coord_bar = toBarycentric( pts );
118 pts += calculateFaceDeformation( coord_bar, entityMap( d,e ), mpl::int_<3>() );
120 final_pts = putInPointset ( final_pts,
121 toEquilateral( pts,
false ),
122 G.interiorRangeById( d, e ) );
128 final_pts = putInPointset ( final_pts,
129 transformPoints<3>( G.pointsBySubEntity( 3, 0, 0 ) ),
130 G.interiorRangeById( 3, 0 ) );
137 this->setName(
"warpblend", Order );
140 ~PointSetWarpBlend() {}
144 std::pair<points_type, points_type> entities;
146 std::map<uint16_type, std::pair<vector_type, vector_type > > axis;
148 points_type equiVertices ()
150 points_type V ( ublas::scalar_matrix<value_type>( Dim, Dim+1, value_type( 0 ) ) );
152 for ( uint16_type i=0; i < 3; i++ )
154 value_type angle = M_PI*( value_type( 7 ) + value_type( 4 )*value_type( i ) )/value_type( 6 );
156 V( 0,i ) = math::cos( angle );
157 V( 1,i ) = math::sin( angle );
160 V( 2,i ) = - math::sqrt( value_type( 2 ) )/value_type( 4 );
163 V *= value_type( 2 )/math::sqrt( value_type( 3 ) );
166 V( 2,3 ) = math::sqrt( value_type( 6 ) )/value_type( 2 );
173 std::vector<uint16_type> indices;
175 for ( uint16_type i=0; i<4; i++ )
177 for ( uint16_type j=0; j<4; j++ )
179 if ( j != ( 18 +i*( -6 + ( i-1 )*( -3 + 4*( i-2 ) ) ) )/6 )
180 indices.push_back( j );
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 );
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 );
194 uint16_type entityMap ( uint16_type top_dim, uint16_type
id )
197 id = ( 12 +
id*( 6 + (
id-1 )*( -9 + 4*(
id-2 ) ) ) )/6;
214 vector_type getVertex ( uint16_type element, uint16_type i )
219 p = ublas::column( entities.first, i );
222 p = ublas::column( entities.second, i );
227 points_type toEquilateral ( points_type pts,
bool toEqui )
229 points_type coord_ref_elem ( Dim, Dim );
230 points_type coord_equi_elem ( Dim, Dim );
232 for ( uint16_type i = 0; i < Dim; i++ )
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 );
238 points_type A ( Dim, Dim );
239 points_type b ( Dim, 1 );
241 LU< points_type > lu( coord_ref_elem );
244 A = ublas::prod( coord_equi_elem , A );
246 ublas::column( b,0 ) = getVertex( 1,0 ) - ublas::prod( A, getVertex( 0,0 ) );
250 pts = ublas::prod( A, pts );
252 for ( uint16_type i=0; i < pts.size2(); i++ )
253 ublas::column( pts, i ) += ublas::column( b,0 );
258 for ( uint16_type i=0; i < pts.size2(); i++ )
259 ublas::column( pts, i ) -= ublas::column( b,0 );
261 LU< points_type > lu( A );
262 pts = lu.solve( pts ) ;
268 points_type toBarycentric ( points_type pts )
270 points_type C( Dim, Dim );
272 for ( uint16_type i = 0; i < Dim; i++ )
273 ublas::column( C, i ) = getVertex( 1, ( Dim - Dim%2 + i )%Dim ) - getVertex( 1,Dim );
275 for ( uint16_type i=0; i < pts.size2(); i++ )
276 ublas::column( pts, i ) -= getVertex( 1,Dim );
278 LU< points_type > lu( C );
280 points_type solution = lu.solve( pts );
282 pts.resize( entities.second.size2(), pts.size2() );
284 ublas::subrange( pts, 1, Dim+1, 0, pts.size2() ) = solution;
286 ublas::row( pts, 0 ) = ublas::scalar_vector<value_type>( pts.size2(), value_type( 1 ) );
288 for ( uint16_type i = 1; i < Dim+1; i++ )
289 ublas::row( pts, 0 ) -= ublas::row( pts, i );
296 vector_type gx( Order+1 );
297 vector_type gw( Order+1 );
299 details::dyna::gausslobattojacobi<value_type, vector_type, vector_type>( Order+1, gw, gx );
304 value_type alpha( uint16_type d )
const
308 double __alpha_2d[ 16 ] = { 0,
326 double __alpha_3d[ 16 ] = { 0,
345 p = ( Order<=15 )?__alpha_2d[ Order ]:( 5./3. );
348 p = ( Order<=15 )?__alpha_3d[ Order ]:1.0;
353 template<
typename AE>
354 vector_type warpFactor( vector_type
const& x, ublas::vector_expression<AE>
const& xout )
const
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 ) );
360 vector_type d( xout().size() );
362 for (
int i = 0; i < nPoints1D; ++i )
364 d = ublas::scalar_vector<value_type>( xout().size(), x( i )-xeq( i ) );
366 for (
int j = 1; j < nPoints1D-1; ++j )
369 d = ( ublas::element_prod( d, xout ) - d *xeq( j ) )/( xeq( i ) - xeq( j ) );
374 d *= value_type( -1 )/( xeq( i ) - xeq( 0 ) );
376 if ( i != nPoints1D-1 )
377 d *= value_type( 1 )/( xeq( i ) - xeq( nPoints1D-1 ) );
386 points_type getCoordinates( points_type
const& pts, uint16_type face_id )
388 std::map<uint16_type, std::vector<uint16_type> >
faces;
390 for ( uint16_type i=0; i<4; i++ )
391 for ( uint16_type j=0; j<4; j++ )
394 faces[i].push_back( j );
397 for ( uint16_type i=2; i<4; i++ )
399 faces[i][1] = faces[i][2];
403 points_type coord ( 3, pts.size2() );
405 for ( uint16_type i=0; i<3; i++ )
406 ublas::row( coord, i ) = ublas::row( pts, faces[face_id][i] );
411 points_type calculateFaceDeformation( points_type
const& coord_bar, mpl::int_<2> )
413 points_type blend ( 3, coord_bar.size2() );
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 ) );
419 blend *= value_type( 4 );
421 points_type warpfactor( 3, coord_bar.size2() );
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 ) );
426 points_type warp( 3, coord_bar.size2() );
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 ) );
432 points_type deformation( ublas::scalar_matrix<value_type>( 2, warp.size2(), value_type( 0 ) ) );
434 for ( uint16_type i=0; i < warp.size1(); i++ )
436 value_type angle = value_type( 2*i )*M_PI/value_type( 3 );
438 ublas::row( deformation, 0 ) += math::cos( angle )*ublas::row( warp, i );
439 ublas::row( deformation, 1 ) += math::sin( angle )*ublas::row( warp, i );
446 points_type calculateFaceDeformation( points_type
const& pts, uint16_type face_id, mpl::int_<3> )
448 points_type coord_bar;
450 coord_bar = getCoordinates( pts, face_id );
452 points_type def = calculateFaceDeformation( coord_bar, mpl::int_<2>() );
454 points_type w ( 3, pts.size2() );
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 );
462 points_type blendFunction( points_type
const& pts )
464 points_type blend ( 4, pts.size2() );
466 for ( uint16_type i=0; i<4; i++ )
468 vector_type aux1 ( ublas::scalar_vector<value_type>( pts.size2(), value_type( 1 ) ) );
469 vector_type aux2 = aux1;
471 for ( uint16_type j=0; j<4; j++ )
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 ) );
480 ublas::row( blend, i ) = ublas::element_div( aux1, aux2 );
483 blend *= value_type( 8 );
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 ) );
492 points_type calculateFaceDeformation( points_type
const& coord_bar, mpl::int_<3> )
494 points_type blend = blendFunction( coord_bar );
496 std::map<uint16_type, points_type > warp;
498 for ( uint16_type face_id = 0; face_id < 4; face_id++ )
500 warp[face_id].resize( 3, coord_bar.size2() );
502 for ( uint16_type i=0; i<3; i++ )
504 vector_type aux = ublas::row( calculateFaceDeformation( coord_bar,
505 entityMap( 2, face_id ),
506 mpl::int_<3>() ), i );
508 ublas::row( warp[face_id], i ) = ublas::element_prod( ublas::row( blend, face_id ), aux );
512 for ( uint16_type face_id = 1; face_id < 4; face_id++ )
513 warp[0] += warp[face_id];
519 points_type transformPoints( points_type pts )
521 pts = toEquilateral( pts,
true );
523 points_type coord_bar = toBarycentric( pts );
525 pts += calculateFaceDeformation( coord_bar, mpl::int_<N>() );
527 return toEquilateral( pts,
false );
530 points_type putInPointset ( points_type final_pts, points_type pts, std::pair<uint16_type, uint16_type> position )
532 ublas::subrange( final_pts, 0, 3, position.first, position.second ) = pts;