30 #ifndef __electrostatic_H
31 #define __electrostatic_H 1
37 template<
class Convex,
40 class PointSetElectrostatic :
public PointSetInterpolation<Convex::nDim, Order, T, Simplex>
45 typedef PointSetWarpBlend<Convex, Order, T> super;
49 static const uint32_type Dim = Convex::nDim;
50 static const uint16_type nPoints1D = Order+1;
52 typedef typename super::return_type return_type;
54 typedef ublas::vector<value_type> vector_type;
56 static const uint32_type topological_dimension = Convex::topological_dimension;
57 static const uint32_type numPoints = super::numPoints;
59 typedef Reference<Convex, Dim, Convex::nOrder, Dim, value_type> reference_convex_type;
61 typedef typename reference_convex_type::points_type points_type;
63 static const uint32_type nbPtsPerFace = super::nbPtsPerFace;
65 typedef std::vector<uint16_type> orbits_type;
67 reference_convex_type RefConv;
69 PointSetElectrostatic(
int interior = 0 )
71 PointSetWarpBlend<Convex, Order, T> G( interior );
73 points_type final_pts = G.points();
76 this->setEid( G.getEid() );
77 this->setPtE( G.getPtE() );
81 entities = std::make_pair( RefConv.vertices(), equiVertices() );
85 points_type Gt( Dim+1, nbPtsPerFace );
89 for ( uint16_type i = 0; i < getPrimalNumber(); i++ )
91 points_type orbit = generateOrbit( getOrbitId( i ),
92 ublas::subrange( getPrimalPts(), 0, 3, i, i+1 ) );
94 ublas::subrange( Gt, 0, 3, p, p+orbit.size2() ) = orbit;
99 points_type pts = toEquilateral( toCartesian( Gt ) );
101 pts = orderInteriorFace( pts );
107 final_pts = putInPointset ( final_pts, pts, G.interiorRangeById( 2, 0 ) );
113 this->setName(
"fekete", Order );
116 ~PointSetElectrostatic() {}
120 std::pair<points_type, points_type> entities;
124 void addOrbit( uint16_type orbit_id, uint16_type num )
126 for ( uint16_type i=0; i<num; i++ )
127 orbits.push_back( orbit_id );
136 else if ( Order == 4 )
139 else if ( Order == 5 )
142 else if ( Order == 6 )
149 else if ( Order == 7 )
155 else if ( Order == 8 )
161 else if ( Order == 9 )
168 else if ( Order == 10 )
174 else if ( Order == 11 )
180 else if ( Order == 12 )
187 else if ( Order == 13 )
194 uint16_type getOrbitId( uint16_type i )
200 uint16_type getPrimalNumber()
202 uint16_type __num[ 16 ] = { 1,
212 return __num[Order-1];
216 points_type getPrimalPts()
218 points_type bar_coord( Dim+1, getPrimalNumber() );
222 bar_coord( 0,0 ) = value_type( 1.0/3.0 );
223 bar_coord( 1,0 ) = value_type( 1.0/3.0 );
226 else if ( Order == 4 )
228 bar_coord( 0,0 ) = value_type( 0.2371200168 );
229 bar_coord( 1,0 ) = value_type( 0.2371200168 );
232 else if ( Order == 5 )
234 bar_coord( 0,0 ) = value_type( 0.410515151 );
235 bar_coord( 1,0 ) = value_type( 0.410515151 );
237 bar_coord( 0,1 ) = value_type( 0.1575181512 );
238 bar_coord( 1,1 ) = value_type( 0.1575181512 );
241 else if ( Order == 6 )
243 bar_coord( 0,0 ) = value_type( 1.0/3 );
244 bar_coord( 1,0 ) = value_type( 1.0/3 );
246 bar_coord( 0,1 ) = value_type( 0.1061169285 );
247 bar_coord( 1,1 ) = value_type( 0.1061169285 );
249 bar_coord( 0,2 ) = value_type( 0.3097982151 );
250 bar_coord( 1,2 ) = value_type( 0.5569099204 );
253 else if ( Order == 7 )
255 bar_coord( 0,0 ) = value_type( 0.4477725053 );
256 bar_coord( 1,0 ) = value_type( 0.4477725053 );
258 bar_coord( 0,1 ) = value_type( 0.2604038024 );
259 bar_coord( 1,1 ) = value_type( 0.2604038024 );
261 bar_coord( 0,2 ) = value_type( 0.0660520784 );
262 bar_coord( 1,2 ) = value_type( 0.0660520784 );
264 bar_coord( 0,3 ) = value_type( 0.2325524777 );
265 bar_coord( 1,3 ) = value_type( 0.6759625951 );
270 else if ( Order == 8 )
272 bar_coord( 0,0 ) = value_type( 0.8444999609 );
273 bar_coord( 1,0 ) = value_type( 0.0777500195 );
275 bar_coord( 0,1 ) = value_type( 0.4683305115 );
276 bar_coord( 1,1 ) = value_type( 0.4683305115 );
278 bar_coord( 0,2 ) = value_type( 0.3853668203 );
279 bar_coord( 1,2 ) = value_type( 0.3853668204 );
281 bar_coord( 0,3 ) = value_type( 0.7172965409 );
282 bar_coord( 1,3 ) = value_type( 0.2335581033 );
284 bar_coord( 0,4 ) = value_type( 0.5853134902 );
285 bar_coord( 1,4 ) = value_type( 0.2667701010 );
288 else if ( Order == 9 )
290 bar_coord( 0,0 ) = value_type( 1.0/3 );
291 bar_coord( 1,0 ) = value_type( 1.0/3 );
293 bar_coord( 0,1 ) = value_type( 0.1704318201 );
294 bar_coord( 1,1 ) = value_type( 0.1704318201 );
296 bar_coord( 0,2 ) = value_type( 0.0600824712 );
297 bar_coord( 1,2 ) = value_type( 0.4699587644 );
299 bar_coord( 0,3 ) = value_type( 0.0489345696 );
300 bar_coord( 1,3 ) = value_type( 0.0489345696 );
302 bar_coord( 0,4 ) = value_type( 0.1784337588 );
303 bar_coord( 1,4 ) = value_type( 0.3252434900 );
305 bar_coord( 0,5 ) = value_type( 0.0588564879 );
306 bar_coord( 1,5 ) = value_type( 0.3010242110 );
308 bar_coord( 0,6 ) = value_type( 0.0551758079 );
309 bar_coord( 1,6 ) = value_type( 0.1543901944 );
312 else if ( Order == 10 )
314 bar_coord( 0,0 ) = value_type( 0.9145236987 );
315 bar_coord( 1,0 ) = value_type( 0.0427381507 );
317 bar_coord( 0,1 ) = value_type( 0.5331019411 );
318 bar_coord( 1,1 ) = value_type( 0.2334490294 );
320 bar_coord( 0,2 ) = value_type( 0.4814795342 );
321 bar_coord( 1,2 ) = value_type( 0.4814795342 );
323 bar_coord( 0,3 ) = value_type( 0.3800851251 );
324 bar_coord( 1,3 ) = value_type( 0.3800851251 );
326 bar_coord( 0,4 ) = value_type( 0.8150971991 );
327 bar_coord( 1,4 ) = value_type( 0.1351329831 );
329 bar_coord( 0,5 ) = value_type( 0.6778669104 );
330 bar_coord( 1,5 ) = value_type( 0.2844305545 );
332 bar_coord( 0,6 ) = value_type( 0.6759450113 );
333 bar_coord( 1,6 ) = value_type( 0.2079572403 );
335 bar_coord( 0,7 ) = value_type( 0.5222323306 );
336 bar_coord( 1,7 ) = value_type( 0.3633472465 );
339 else if ( Order == 11 )
341 bar_coord( 0,0 ) = value_type( 0.9201760661 );
342 bar_coord( 1,0 ) = value_type( 0.0399119670 );
344 bar_coord( 0,1 ) = value_type( 0.8097416696 );
345 bar_coord( 1,1 ) = value_type( 0.0951291652 );
347 bar_coord( 0,2 ) = value_type( 0.4216558161 );
348 bar_coord( 1,2 ) = value_type( 0.2891720920 );
350 bar_coord( 0,3 ) = value_type( 0.4200100315 );
351 bar_coord( 1,3 ) = value_type( 0.4200100315 );
353 bar_coord( 0,4 ) = value_type( 0.4832770031 );
354 bar_coord( 1,4 ) = value_type( 0.4832770031 );
356 bar_coord( 0,5 ) = value_type( 0.8236881237 );
357 bar_coord( 1,5 ) = value_type( 0.1452587341 );
359 bar_coord( 0,6 ) = value_type( 0.7030268141 );
360 bar_coord( 1,6 ) = value_type( 0.2021386640 );
362 bar_coord( 0,7 ) = value_type( 0.6642752329 );
363 bar_coord( 1,7 ) = value_type( 0.3066778199 );
365 bar_coord( 0,8 ) = value_type( 0.5605605456 );
366 bar_coord( 1,8 ) = value_type( 0.3510551601 );
368 bar_coord( 0,9 ) = value_type( 0.5584153138 );
369 bar_coord( 1,9 ) = value_type( 0.2661283688 );
373 else if ( Order == 12 )
375 bar_coord( 0,0 ) = value_type( 1.0/3 );
376 bar_coord( 1,0 ) = value_type( 1.0/3 );
378 bar_coord( 0,1 ) = value_type( 0.1988883477 );
379 bar_coord( 1,1 ) = value_type( 0.4005558262 );
381 bar_coord( 0,2 ) = value_type( 0.2618405201 );
382 bar_coord( 1,2 ) = value_type( 0.2618405201 );
384 bar_coord( 0,3 ) = value_type( 0.0807386775 );
385 bar_coord( 1,3 ) = value_type( 0.0807386775 );
387 bar_coord( 0,4 ) = value_type( 0.0336975736 );
388 bar_coord( 1,4 ) = value_type( 0.0336975736 );
390 bar_coord( 0,5 ) = value_type( 0.1089969290 );
391 bar_coord( 1,5 ) = value_type( 0.3837518758 );
393 bar_coord( 0,6 ) = value_type( 0.1590834479 );
394 bar_coord( 1,6 ) = value_type( 0.2454317980 );
396 bar_coord( 0,7 ) = value_type( 0.0887037176 );
397 bar_coord( 1,7 ) = value_type( 0.1697134458 );
399 bar_coord( 0,8 ) = value_type( 0.0302317829 );
400 bar_coord( 1,8 ) = value_type( 0.4071849276 );
402 bar_coord( 0,9 ) = value_type( 0.0748751152 );
403 bar_coord( 1,9 ) = value_type( 0.2874821712 );
405 bar_coord( 0,10 ) = value_type( 0.0250122615 );
406 bar_coord( 1,10 ) = value_type( 0.2489279690 );
408 bar_coord( 0,11 ) = value_type( 0.0262645218 );
409 bar_coord( 1,11 ) = value_type( 0.1206826354 );
414 for ( uint16_type i=0; i < getPrimalNumber(); i++ )
415 bar_coord( 2,i ) = 1 - bar_coord( 0,i ) - bar_coord( 1,i );
420 points_type cyclic( points_type pts, uint16_type permOrder )
424 for ( uint16_type i=0; i<pts.size2(); i++ )
426 value_type a = pts( 0,i );
427 value_type b = pts( 1,i );
428 value_type c = pts( 2,i );
437 pts = cyclic( pts, permOrder-1 );
443 points_type symmetries( points_type pts, uint16_type permOrder )
446 pts = cyclic( pts, permOrder );
450 for ( uint16_type i=0; i<pts.size2(); i++ )
452 value_type a = pts( 0,i );
453 value_type b = pts( 1,i );
454 value_type c = pts( 2,i );
460 pts = cyclic( pts, permOrder-3 );
468 points_type generateOrbit( uint16_type numOrbit, points_type orbitPoint )
470 points_type orbit( Dim+1, numOrbit );
477 for ( uint16_type i=0; i<numOrbit; i++ )
478 ublas::subrange( orbit, 0, 3, i, i+1 ) = symmetries( orbitPoint, i );
485 points_type equiVertices ()
487 points_type V ( ublas::scalar_matrix<value_type>( Dim, Dim+1, value_type( 0 ) ) );
489 for ( uint16_type i=0; i < 3; i++ )
491 value_type angle = M_PI*( value_type( 7 ) + value_type( 4 )*value_type( i ) )/value_type( 6 );
493 V( 0,i ) = math::cos( angle );
494 V( 1,i ) = math::sin( angle );
497 V *= value_type( 2 )/math::sqrt( value_type( 3 ) );
502 vector_type getVertex ( uint16_type element, uint16_type i )
507 p = ublas::column( entities.first, i );
510 p = ublas::column( entities.second, i );
515 points_type toCartesian ( points_type pts )
517 points_type C( Dim, Dim );
519 for ( uint16_type i = 0; i < Dim; i++ )
520 ublas::column( C, i ) = getVertex( 1, ( Dim - Dim%2 + i )%Dim ) - getVertex( 1,Dim );
522 points_type bar_coord = ublas::subrange( pts, 1, Dim+1, 0, pts.size2() );
524 points_type Gt = ublas::prod( C, bar_coord );
526 for ( uint16_type i=0; i < pts.size2(); i++ )
527 ublas::column( Gt, i ) += getVertex( 1,Dim );
532 points_type toEquilateral ( points_type pts )
534 points_type coord_ref_elem ( Dim, Dim );
535 points_type coord_equi_elem ( Dim, Dim );
537 for ( uint16_type i = 0; i < Dim; i++ )
539 ublas::column( coord_ref_elem, i ) = getVertex( 0,Dim ) - getVertex( 0,i );
540 ublas::column( coord_equi_elem, i ) = getVertex( 1,Dim ) - getVertex( 1,i );
543 points_type A ( Dim, Dim );
544 points_type b ( Dim, 1 );
546 LU< points_type > lu( coord_ref_elem );
549 A = ublas::prod( coord_equi_elem , A );
551 ublas::column( b,0 ) = getVertex( 1,0 ) - ublas::prod( A, getVertex( 0,0 ) );
553 for ( uint16_type i=0; i < pts.size2(); i++ )
554 ublas::column( pts, i ) -= ublas::column( b,0 );
556 LU< points_type > lu2( A );
557 pts = lu2.solve( pts ) ;
562 points_type putInPointset ( points_type final_pts, points_type pts, std::pair<uint16_type, uint16_type> position )
564 ublas::subrange( final_pts, 0, 2, position.first, position.second ) = pts;
570 static bool order( vector_type a, vector_type b )
572 return a( n ) < b( n );
575 points_type orderInteriorFace( points_type pts )
577 std::vector<vector_type> aux( pts.size2() );
579 for ( uint16_type i=0; i<pts.size2(); i++ )
580 aux[i] = ublas::column( pts, i );
582 std::sort( aux.begin(), aux.end(), &order<1> );
584 for ( uint16_type p=0, i=Order-3; i>0; i-- )
586 std::sort( aux.begin()+p, aux.begin()+p+i+1, &order<0> );
590 for ( uint16_type i=0; i<pts.size2(); i++ )
591 ublas::column( pts, i ) = aux[i];
void setPoints(nodes_type const &pts)
Definition: pointset.hpp:267