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
electrostatic.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 __electrostatic_H
31 #define __electrostatic_H 1
32 
34 
35 namespace Feel
36 {
37 template< class Convex,
38  uint16_type Order,
39  typename T = double >
40 class PointSetElectrostatic : public PointSetInterpolation<Convex::nDim, Order, T, Simplex>
41 {
42 
43 public :
44 
45  typedef PointSetWarpBlend<Convex, Order, T> super;
46 
47  typedef T value_type;
48 
49  static const uint32_type Dim = Convex::nDim;
50  static const uint16_type nPoints1D = Order+1;
51 
52  typedef typename super::return_type return_type;
53 
54  typedef ublas::vector<value_type> vector_type;
55 
56  static const uint32_type topological_dimension = Convex::topological_dimension;
57  static const uint32_type numPoints = super::numPoints;
58 
59  typedef Reference<Convex, Dim, Convex::nOrder, Dim, value_type> reference_convex_type;
60 
61  typedef typename reference_convex_type::points_type points_type;
62 
63  static const uint32_type nbPtsPerFace = super::nbPtsPerFace;
64 
65  typedef std::vector<uint16_type> orbits_type;
66 
67  reference_convex_type RefConv;
68 
69  PointSetElectrostatic( int interior = 0 )
70  {
71  PointSetWarpBlend<Convex, Order, T> G( interior );
72 
73  points_type final_pts = G.points();
74 
75  //Copies information about WarpBlend points
76  this->setEid( G.getEid() );
77  this->setPtE( G.getPtE() );
78 
79  if ( Order > 2 )
80  {
81  entities = std::make_pair( RefConv.vertices(), equiVertices() );
82 
83  defineOrbits();
84 
85  points_type Gt( Dim+1, nbPtsPerFace );
86 
87  uint16_type p=0;
88 
89  for ( uint16_type i = 0; i < getPrimalNumber(); i++ )
90  {
91  points_type orbit = generateOrbit( getOrbitId( i ),
92  ublas::subrange( getPrimalPts(), 0, 3, i, i+1 ) );
93 
94  ublas::subrange( Gt, 0, 3, p, p+orbit.size2() ) = orbit;
95 
96  p += orbit.size2();
97  }
98 
99  points_type pts = toEquilateral( toCartesian( Gt ) );
100 
101  pts = orderInteriorFace( pts );
102 
103  if ( interior )
104  final_pts = pts;
105 
106  else
107  final_pts = putInPointset ( final_pts, pts, G.interiorRangeById( 2, 0 ) );
108 
109  }
110 
111  this->setPoints( final_pts );
112 
113  this->setName( "fekete", Order );
114  }
115 
116  ~PointSetElectrostatic() {}
117 
118 private :
119 
120  std::pair<points_type, points_type> entities;
121 
122  orbits_type orbits;
123 
124  void addOrbit( uint16_type orbit_id, uint16_type num )
125  {
126  for ( uint16_type i=0; i<num; i++ )
127  orbits.push_back( orbit_id );
128  }
129 
130  //only define orbits for interior points; all the other we know that are the gausslobatto points in the edges
131  void defineOrbits()
132  {
133  if ( Order == 3 )
134  addOrbit( 1, 1 );
135 
136  else if ( Order == 4 )
137  addOrbit( 3, 1 );
138 
139  else if ( Order == 5 )
140  addOrbit( 3, 2 );
141 
142  else if ( Order == 6 )
143  {
144  addOrbit( 1, 1 );
145  addOrbit( 3, 1 );
146  addOrbit( 6, 1 );
147  }
148 
149  else if ( Order == 7 )
150  {
151  addOrbit( 3, 3 );
152  addOrbit( 6, 1 );
153  }
154 
155  else if ( Order == 8 )
156  {
157  addOrbit( 3, 3 );
158  addOrbit( 6, 2 );
159  }
160 
161  else if ( Order == 9 )
162  {
163  addOrbit( 1, 1 );
164  addOrbit( 3, 3 );
165  addOrbit( 6, 3 );
166  }
167 
168  else if ( Order == 10 )
169  {
170  addOrbit( 3, 4 );
171  addOrbit( 6, 4 );
172  }
173 
174  else if ( Order == 11 )
175  {
176  addOrbit( 3, 5 );
177  addOrbit( 6, 5 );
178  }
179 
180  else if ( Order == 12 )
181  {
182  addOrbit( 1, 1 );
183  addOrbit( 3, 4 );
184  addOrbit( 6, 7 );
185  }
186 
187  else if ( Order == 13 )
188  {
189  addOrbit( 3, 6 );
190  addOrbit( 6, 8 );
191  }
192  }
193 
194  uint16_type getOrbitId( uint16_type i )
195  {
196  return orbits[i];
197  }
198 
199  // returns the number of points defined in the triangle with which we calculate the orbits
200  uint16_type getPrimalNumber()
201  {
202  uint16_type __num[ 16 ] = { 1,
203  1,
204  1,
205  1,
206  2,
207  3,
208  4,
209  5
210  };
211 
212  return __num[Order-1];
213  }
214 
215  // defines the points with which we calculate the orbits
216  points_type getPrimalPts()
217  {
218  points_type bar_coord( Dim+1, getPrimalNumber() );
219 
220  if ( Order == 3 )
221  {
222  bar_coord( 0,0 ) = value_type( 1.0/3.0 );
223  bar_coord( 1,0 ) = value_type( 1.0/3.0 );
224  }
225 
226  else if ( Order == 4 )
227  {
228  bar_coord( 0,0 ) = value_type( 0.2371200168 );
229  bar_coord( 1,0 ) = value_type( 0.2371200168 );
230  }
231 
232  else if ( Order == 5 )
233  {
234  bar_coord( 0,0 ) = value_type( 0.410515151 );
235  bar_coord( 1,0 ) = value_type( 0.410515151 );
236 
237  bar_coord( 0,1 ) = value_type( 0.1575181512 );
238  bar_coord( 1,1 ) = value_type( 0.1575181512 );
239  }
240 
241  else if ( Order == 6 )
242  {
243  bar_coord( 0,0 ) = value_type( 1.0/3 );
244  bar_coord( 1,0 ) = value_type( 1.0/3 );
245 
246  bar_coord( 0,1 ) = value_type( 0.1061169285 );
247  bar_coord( 1,1 ) = value_type( 0.1061169285 );
248 
249  bar_coord( 0,2 ) = value_type( 0.3097982151 );
250  bar_coord( 1,2 ) = value_type( 0.5569099204 );
251  }
252 
253  else if ( Order == 7 )
254  {
255  bar_coord( 0,0 ) = value_type( 0.4477725053 );
256  bar_coord( 1,0 ) = value_type( 0.4477725053 );
257 
258  bar_coord( 0,1 ) = value_type( 0.2604038024 );
259  bar_coord( 1,1 ) = value_type( 0.2604038024 );
260 
261  bar_coord( 0,2 ) = value_type( 0.0660520784 );
262  bar_coord( 1,2 ) = value_type( 0.0660520784 );
263 
264  bar_coord( 0,3 ) = value_type( 0.2325524777 );
265  bar_coord( 1,3 ) = value_type( 0.6759625951 );
266  }
267 
268 #if 0
269 
270  else if ( Order == 8 )
271  {
272  bar_coord( 0,0 ) = value_type( 0.8444999609 );
273  bar_coord( 1,0 ) = value_type( 0.0777500195 );
274 
275  bar_coord( 0,1 ) = value_type( 0.4683305115 );
276  bar_coord( 1,1 ) = value_type( 0.4683305115 );
277 
278  bar_coord( 0,2 ) = value_type( 0.3853668203 );
279  bar_coord( 1,2 ) = value_type( 0.3853668204 );
280 
281  bar_coord( 0,3 ) = value_type( 0.7172965409 );
282  bar_coord( 1,3 ) = value_type( 0.2335581033 );
283 
284  bar_coord( 0,4 ) = value_type( 0.5853134902 );
285  bar_coord( 1,4 ) = value_type( 0.2667701010 );
286  }
287 
288  else if ( Order == 9 )
289  {
290  bar_coord( 0,0 ) = value_type( 1.0/3 );
291  bar_coord( 1,0 ) = value_type( 1.0/3 );
292 
293  bar_coord( 0,1 ) = value_type( 0.1704318201 );
294  bar_coord( 1,1 ) = value_type( 0.1704318201 );
295 
296  bar_coord( 0,2 ) = value_type( 0.0600824712 );
297  bar_coord( 1,2 ) = value_type( 0.4699587644 );
298 
299  bar_coord( 0,3 ) = value_type( 0.0489345696 );
300  bar_coord( 1,3 ) = value_type( 0.0489345696 );
301 
302  bar_coord( 0,4 ) = value_type( 0.1784337588 );
303  bar_coord( 1,4 ) = value_type( 0.3252434900 );
304 
305  bar_coord( 0,5 ) = value_type( 0.0588564879 );
306  bar_coord( 1,5 ) = value_type( 0.3010242110 );
307 
308  bar_coord( 0,6 ) = value_type( 0.0551758079 );
309  bar_coord( 1,6 ) = value_type( 0.1543901944 );
310  }
311 
312  else if ( Order == 10 )
313  {
314  bar_coord( 0,0 ) = value_type( 0.9145236987 );
315  bar_coord( 1,0 ) = value_type( 0.0427381507 );
316 
317  bar_coord( 0,1 ) = value_type( 0.5331019411 );
318  bar_coord( 1,1 ) = value_type( 0.2334490294 );
319 
320  bar_coord( 0,2 ) = value_type( 0.4814795342 );
321  bar_coord( 1,2 ) = value_type( 0.4814795342 );
322 
323  bar_coord( 0,3 ) = value_type( 0.3800851251 );
324  bar_coord( 1,3 ) = value_type( 0.3800851251 );
325 
326  bar_coord( 0,4 ) = value_type( 0.8150971991 );
327  bar_coord( 1,4 ) = value_type( 0.1351329831 );
328 
329  bar_coord( 0,5 ) = value_type( 0.6778669104 );
330  bar_coord( 1,5 ) = value_type( 0.2844305545 );
331 
332  bar_coord( 0,6 ) = value_type( 0.6759450113 );
333  bar_coord( 1,6 ) = value_type( 0.2079572403 );
334 
335  bar_coord( 0,7 ) = value_type( 0.5222323306 );
336  bar_coord( 1,7 ) = value_type( 0.3633472465 );
337  }
338 
339  else if ( Order == 11 )
340  {
341  bar_coord( 0,0 ) = value_type( 0.9201760661 );
342  bar_coord( 1,0 ) = value_type( 0.0399119670 );
343 
344  bar_coord( 0,1 ) = value_type( 0.8097416696 );
345  bar_coord( 1,1 ) = value_type( 0.0951291652 );
346 
347  bar_coord( 0,2 ) = value_type( 0.4216558161 );
348  bar_coord( 1,2 ) = value_type( 0.2891720920 );
349 
350  bar_coord( 0,3 ) = value_type( 0.4200100315 );
351  bar_coord( 1,3 ) = value_type( 0.4200100315 );
352 
353  bar_coord( 0,4 ) = value_type( 0.4832770031 );
354  bar_coord( 1,4 ) = value_type( 0.4832770031 );
355 
356  bar_coord( 0,5 ) = value_type( 0.8236881237 );
357  bar_coord( 1,5 ) = value_type( 0.1452587341 );
358 
359  bar_coord( 0,6 ) = value_type( 0.7030268141 );
360  bar_coord( 1,6 ) = value_type( 0.2021386640 );
361 
362  bar_coord( 0,7 ) = value_type( 0.6642752329 );
363  bar_coord( 1,7 ) = value_type( 0.3066778199 );
364 
365  bar_coord( 0,8 ) = value_type( 0.5605605456 );
366  bar_coord( 1,8 ) = value_type( 0.3510551601 );
367 
368  bar_coord( 0,9 ) = value_type( 0.5584153138 );
369  bar_coord( 1,9 ) = value_type( 0.2661283688 );
370 
371  }
372 
373  else if ( Order == 12 )
374  {
375  bar_coord( 0,0 ) = value_type( 1.0/3 );
376  bar_coord( 1,0 ) = value_type( 1.0/3 );
377 
378  bar_coord( 0,1 ) = value_type( 0.1988883477 );
379  bar_coord( 1,1 ) = value_type( 0.4005558262 );
380 
381  bar_coord( 0,2 ) = value_type( 0.2618405201 );
382  bar_coord( 1,2 ) = value_type( 0.2618405201 );
383 
384  bar_coord( 0,3 ) = value_type( 0.0807386775 );
385  bar_coord( 1,3 ) = value_type( 0.0807386775 );
386 
387  bar_coord( 0,4 ) = value_type( 0.0336975736 );
388  bar_coord( 1,4 ) = value_type( 0.0336975736 );
389 
390  bar_coord( 0,5 ) = value_type( 0.1089969290 );
391  bar_coord( 1,5 ) = value_type( 0.3837518758 );
392 
393  bar_coord( 0,6 ) = value_type( 0.1590834479 );
394  bar_coord( 1,6 ) = value_type( 0.2454317980 );
395 
396  bar_coord( 0,7 ) = value_type( 0.0887037176 );
397  bar_coord( 1,7 ) = value_type( 0.1697134458 );
398 
399  bar_coord( 0,8 ) = value_type( 0.0302317829 );
400  bar_coord( 1,8 ) = value_type( 0.4071849276 );
401 
402  bar_coord( 0,9 ) = value_type( 0.0748751152 );
403  bar_coord( 1,9 ) = value_type( 0.2874821712 );
404 
405  bar_coord( 0,10 ) = value_type( 0.0250122615 );
406  bar_coord( 1,10 ) = value_type( 0.2489279690 );
407 
408  bar_coord( 0,11 ) = value_type( 0.0262645218 );
409  bar_coord( 1,11 ) = value_type( 0.1206826354 );
410  }
411 
412 #endif
413 
414  for ( uint16_type i=0; i < getPrimalNumber(); i++ )
415  bar_coord( 2,i ) = 1 - bar_coord( 0,i ) - bar_coord( 1,i );
416 
417  return bar_coord;
418  }
419 
420  points_type cyclic( points_type pts, uint16_type permOrder )
421  {
422  if ( permOrder > 0 )
423  {
424  for ( uint16_type i=0; i<pts.size2(); i++ )
425  {
426  value_type a = pts( 0,i );
427  value_type b = pts( 1,i );
428  value_type c = pts( 2,i );
429 
430  pts( 0,i ) = c;
431  pts( 1,i ) = a;
432  pts( 2,i ) = b;
433  }
434  }
435 
436  if ( permOrder > 1 )
437  pts = cyclic( pts, permOrder-1 );
438 
439  return pts;
440  }
441 
442 
443  points_type symmetries( points_type pts, uint16_type permOrder )
444  {
445  if ( permOrder < 3 )
446  pts = cyclic( pts, permOrder );
447 
448  else
449  {
450  for ( uint16_type i=0; i<pts.size2(); i++ )
451  {
452  value_type a = pts( 0,i );
453  value_type b = pts( 1,i );
454  value_type c = pts( 2,i );
455 
456  pts( 0,i ) = a;
457  pts( 1,i ) = c;
458  pts( 2,i ) = b;
459 
460  pts = cyclic( pts, permOrder-3 );
461  }
462  }
463 
464  return pts;
465  }
466 
467 
468  points_type generateOrbit( uint16_type numOrbit, points_type orbitPoint )
469  {
470  points_type orbit( Dim+1, numOrbit );
471 
472  if ( numOrbit == 1 )
473  orbit = orbitPoint;
474 
475  else
476  {
477  for ( uint16_type i=0; i<numOrbit; i++ )
478  ublas::subrange( orbit, 0, 3, i, i+1 ) = symmetries( orbitPoint, i );
479  }
480 
481  return orbit;
482  }
483 
484 
485  points_type equiVertices ()
486  {
487  points_type V ( ublas::scalar_matrix<value_type>( Dim, Dim+1, value_type( 0 ) ) );
488 
489  for ( uint16_type i=0; i < 3; i++ )
490  {
491  value_type angle = M_PI*( value_type( 7 ) + value_type( 4 )*value_type( i ) )/value_type( 6 );
492 
493  V( 0,i ) = math::cos( angle );
494  V( 1,i ) = math::sin( angle );
495  }
496 
497  V *= value_type( 2 )/math::sqrt( value_type( 3 ) );
498 
499  return V;
500  }
501 
502  vector_type getVertex ( uint16_type element, uint16_type i )
503  {
504  vector_type p;
505 
506  if ( element == 0 )
507  p = ublas::column( entities.first, i );
508 
509  else
510  p = ublas::column( entities.second, i );
511 
512  return p;
513  }
514 
515  points_type toCartesian ( points_type pts )
516  {
517  points_type C( Dim, Dim );
518 
519  for ( uint16_type i = 0; i < Dim; i++ )
520  ublas::column( C, i ) = getVertex( 1, ( Dim - Dim%2 + i )%Dim ) - getVertex( 1,Dim );
521 
522  points_type bar_coord = ublas::subrange( pts, 1, Dim+1, 0, pts.size2() );
523 
524  points_type Gt = ublas::prod( C, bar_coord );
525 
526  for ( uint16_type i=0; i < pts.size2(); i++ )
527  ublas::column( Gt, i ) += getVertex( 1,Dim );
528 
529  return Gt;
530  }
531 
532  points_type toEquilateral ( points_type pts )
533  {
534  points_type coord_ref_elem ( Dim, Dim );
535  points_type coord_equi_elem ( Dim, Dim );
536 
537  for ( uint16_type i = 0; i < Dim; i++ )
538  {
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 );
541  }
542 
543  points_type A ( Dim, Dim );
544  points_type b ( Dim, 1 );
545 
546  LU< points_type > lu( coord_ref_elem );
547  lu.inverse( A );
548 
549  A = ublas::prod( coord_equi_elem , A );
550 
551  ublas::column( b,0 ) = getVertex( 1,0 ) - ublas::prod( A, getVertex( 0,0 ) );
552 
553  for ( uint16_type i=0; i < pts.size2(); i++ )
554  ublas::column( pts, i ) -= ublas::column( b,0 );
555 
556  LU< points_type > lu2( A );
557  pts = lu2.solve( pts ) ;
558 
559  return pts;
560  }
561 
562  points_type putInPointset ( points_type final_pts, points_type pts, std::pair<uint16_type, uint16_type> position )
563  {
564  ublas::subrange( final_pts, 0, 2, position.first, position.second ) = pts;
565 
566  return final_pts;
567  }
568 
569  template<int n>
570  static bool order( vector_type a, vector_type b )
571  {
572  return a( n ) < b( n );
573  }
574 
575  points_type orderInteriorFace( points_type pts )
576  {
577  std::vector<vector_type> aux( pts.size2() );
578 
579  for ( uint16_type i=0; i<pts.size2(); i++ )
580  aux[i] = ublas::column( pts, i );
581 
582  std::sort( aux.begin(), aux.end(), &order<1> );
583 
584  for ( uint16_type p=0, i=Order-3; i>0; i-- )
585  {
586  std::sort( aux.begin()+p, aux.begin()+p+i+1, &order<0> );
587  p += i+1;
588  }
589 
590  for ( uint16_type i=0; i<pts.size2(); i++ )
591  ublas::column( pts, i ) = aux[i];
592 
593  return pts;
594  }
595 };
596 
597 } // Feel
598 #endif /* __Electrostatic_H */
void setPoints(nodes_type const &pts)
Definition: pointset.hpp:267

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