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
principal.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): Gilles Steiner <gilles.steiner@epfl.ch>
6  Date: 2005-12-14
7 
8  Copyright (C) 2005,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 __principal_H
30 #define __principal_H 1
31 
32 #include <vector>
33 
34 
35 #include <feel/feelalg/lu.hpp>
37 #include <feel/feelpoly/policy.hpp>
38 
39 namespace Feel
40 {
41 
57 template< uint16_type Degree, typename T = double,
58  template<class> class StoragePolicy = StorageUBlas>
59 class Principal
60 {
61 public:
62 
63  static const uint16_type nOrder = Degree;
64 
68 
70 
71  /*
72  * numerical type
73  */
74  typedef T value_type;
75 
76  /*
77  * storage policy
78  */
79  typedef StoragePolicy<value_type> storage_policy;
80  typedef typename storage_policy::matrix_type matrix_type;
81  typedef typename storage_policy::vector_matrix_type vector_matrix_type;
82  typedef typename storage_policy::matrix_node_type matrix_node_type;
83  typedef typename storage_policy::node_type node_type;
84  typedef typename storage_policy::vector_vector_matrix_type vector_vector_matrix_type;
85  typedef typename storage_policy::vector_type vector_type;
86 
88 
92 
93  Principal()
94  : M_a( value_type( 1.0 ) ), M_b( value_type( 1.0 ) )
95  {}
96 
97  Principal( value_type a,value_type b )
98  : M_a( a ), M_b( b )
99  {}
100 
101 
102  ~Principal()
103  {}
104 
106 
107 
108  self_type const& operator=( self_type const& d )
109  {
110  if ( this != &d )
111  {
112  M_a = d.M_a;
113  M_b = d.M_b;
114  }
115 
116  return *this;
117  }
118 
122 
123  uint16_type degree() const
124  {
125  return nOrder;
126  }
127 
129 
133 
134  value_type a() const
135  {
136  return M_a;
137  }
138  value_type b() const
139  {
140  return M_b;
141  }
142 
151  matrix_type evaluate_1( vector_type const& __pts ) const;
152 
155  vector_matrix_type evaluate_2( vector_type const& __pts ) const;
156 
157 
160  vector_vector_matrix_type evaluate_3( vector_type const& __pts ) const;
162 
163 
164  matrix_type derivate_1( vector_type const& __pts ) const;
165 
166  vector_matrix_type derivate_2( vector_type const& __pts ) const;
167 
168  vector_vector_matrix_type derivate_3( vector_type const& __pts ) const;
169 
170 private:
171 
172  value_type M_a;
173  value_type M_b;
174 
175 
176 };
177 
178 template<uint16_type Degree,
179  typename T,
180  template<class> class StoragePolicy>
181 typename Principal<Degree, T, StoragePolicy>::matrix_type
182 Principal<Degree, T, StoragePolicy>::evaluate_1( vector_type const& __pts ) const
183 {
184  // std::cout <<"[principal]1D evaluation ..."<< std::endl;
185  matrix_type J ( JacobiBatchEvaluation<nOrder-1,value_type>( M_a, M_b, __pts ) );
186  matrix_type D( nOrder+1,__pts.size() );
187 
188  vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
189 
190  ublas::row( D,0 ) = value_type( 0.5 )*( ones - __pts );
191  ublas::row( D,nOrder ) = value_type( 0.5 )*( ones + __pts );
192 
193  vector_type tmp( ublas::element_prod( ublas::row( D,0 ),ublas::row( D,D.size1()-1 ) ) );
194 
195  for ( uint16_type i = 1; i < D.size1()-1; ++i )
196  ublas::row( D, i ) = ublas::element_prod( tmp, ublas::row( J,i-1 ) );
197 
198  return D;
199 }
200 
201 
202 template<uint16_type Degree,
203  typename T,
204  template<class> class StoragePolicy>
205 typename Principal<Degree, T, StoragePolicy>::vector_matrix_type
206 Principal<Degree, T, StoragePolicy>::evaluate_2( vector_type const& __pts ) const
207 {
208  // std::cout <<"[principal]2D evaluation ..."<< std::endl;
209  matrix_type psi_1( evaluate_1( __pts ) );
210 
211  vector_matrix_type m( nOrder+1 );
212  vector_matrix_type J( nOrder-1 );
213 
214  m[0].resize( nOrder+1 ,__pts.size() );
215  m[0] = psi_1;
216  m[nOrder] = psi_1;
217 
218  vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
219 
220  vector_type tmp1 = ( ones - __pts )/value_type( 2.0 ); // (1-x)/2
221  vector_type tmp2 = ( ones + __pts )/value_type( 2.0 ); // (1+x)/2
222  vector_type tmp_i = tmp1;
223 
224  for ( uint16_type i= 1; i <= J.size() ; ++i )
225  {
226  m[i].resize( nOrder ,__pts.size() );
227  J[i-1] = JacobiBatchEvaluation<nOrder-1,value_type>( value_type( 2.0*i+1.0 ), M_b, __pts );
228 
229  tmp_i=ublas::element_prod( tmp_i,tmp1 );
230 
231  ublas::row( m[i],0 ) = tmp_i;
232  vector_type tmp3( ublas::element_prod( tmp2,tmp_i ) );
233 
234  for ( int16_type j=1; j < nOrder; ++j )
235  {
236  ublas::row( m[i],j ) = ublas::element_prod( tmp3, ublas::row( J[i-1],j-1 ) );
237  }
238  }
239 
240  return m;
241 }
242 
243 
244 template<uint16_type Degree,
245  typename T,
246  template<class> class StoragePolicy>
247 typename Principal<Degree, T, StoragePolicy>::vector_vector_matrix_type
248 Principal<Degree, T, StoragePolicy>::evaluate_3( vector_type const& __pts ) const
249 {
250  // std::cout <<"[principal]3D evaluation ..."<< std::endl;
251  vector_vector_matrix_type m( nOrder+1 );
252 
253  vector_matrix_type psi_2( evaluate_2( __pts ) );
254 
255  // i=0, 0 <= j <= nOrder , 0 <= k <= nOrder
256  m[0].resize( nOrder+1 );
257  m[0] = psi_2;
258 
259  // i=nOrder, 0 <= j <= nOrder , 0 <= k <= nOrder
260  m[nOrder].resize( nOrder+1 );
261  m[nOrder] = psi_2;
262 
263  vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
264  vector_type tmp1 = ( ones - __pts )/value_type( 2.0 ); // (1-x)/2
265  vector_type tmp2 = ( ones + __pts )/value_type( 2.0 ); // (1+x)/2
266  vector_type tmp_i = tmp1;
267 
268  vector_vector_matrix_type J( nOrder+1 );
269 
270  for ( int16_type i=1; i < nOrder; ++i )
271  {
272  m[i].resize( nOrder+1 );
273  J[i].resize( nOrder+1 );
274 
275  // 1 <= i <= nOrder-1, j = 0 , 0 <= k <= nOrder
276  m[i][0].resize( psi_2[i].size1(), psi_2[i].size2() );
277  m[i][0] = psi_2[i];
278 
279  // 1 <= i <= nOrder-1, j = nOrder , 0 <= k <= nOrder
280  m[i][nOrder].resize( psi_2[i].size1(),psi_2[i].size2() );
281  m[i][nOrder] = psi_2[i];
282 
283  tmp_i=ublas::element_prod( tmp_i,tmp1 );
284 
285  vector_type tmp_i_j( tmp_i );
286 
287  for ( int16_type j=1; j < nOrder; ++j )
288  {
289  // 1 <= i <= nOrder-1, 1 <= j <= nOrder-1 , 0 <= k <= nOrder-1
290  m[i][j].resize( nOrder,__pts.size() );
291  J[i][j] = JacobiBatchEvaluation<nOrder-1,value_type>( ( 2.0*i+2.0*j+1.0 ), 1.0, __pts );
292  tmp_i_j = ublas::element_prod( tmp_i_j,tmp1 );
293 
294  ublas::row( m[i][j],0 ) = tmp_i_j; // k=0
295 
296  for ( int16_type k=1; k < nOrder; ++k ) // 1 <= k <= nOrder-1
297  {
298  vector_type tmp3( ublas::element_prod( tmp2,tmp_i_j ) );
299  ublas::row( m[i][j],k ) = ublas::element_prod( tmp3, ublas::row( J[i][j],k-1 ) );
300  }
301  }
302  }
303 
304  return m;
305 }
306 
307 
308 template<uint16_type Degree,
309  typename T,
310  template<class> class StoragePolicy>
311 typename Principal<Degree, T, StoragePolicy>::matrix_type
313 {
314  // std::cout <<"[principal]1D derivation ..."<< std::endl;
315  matrix_type D( nOrder+1, __pts.size() );
316  matrix_type J ( JacobiBatchDerivation<nOrder-1,value_type>( M_a, M_b, __pts ) );
317  vector_type demi( ublas::scalar_vector<value_type>( __pts.size(), value_type( 0.5 ) ) );
318 
319  ublas::row( D,0 ) = -demi;
320  ublas::row( D,nOrder ) = demi;
321 
322  vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
323 
324 
325  matrix_type E ( JacobiBatchEvaluation<nOrder-1,value_type>( M_a, M_b, __pts ) );
326 
327  vector_type tmp( ublas::element_prod( 0.5*( ones - __pts ),0.5*( ones + __pts ) ) );
328 
329  for ( uint16_type i = 1; i < D.size1()-1; ++i )
330  ublas::row( D, i ) = ublas::element_prod( tmp, ublas::row( J,i-1 ) ) - 0.5*ublas::element_prod( __pts, ublas::row( E,i-1 ) ); // \f[ = \frac{1-x}{2}\frac{1+x}{2} \frac{d}{dx}P_{i-1}^{(1,1)}(x) - \frac{x}{2}P_{i-1}^{(1,1)}(x) \f]
331 
332  return D;
333 }
334 
335 template<uint16_type Degree,
336  typename T,
337  template<class> class StoragePolicy>
338 typename Principal<Degree, T, StoragePolicy>::vector_matrix_type
339 Principal<Degree, T, StoragePolicy>::derivate_2( vector_type const& __pts ) const
340 {
341  // std::cout <<"[principal]2D derivation ..."<< std::endl;
342  matrix_type dpsi_1( derivate_1( __pts ) );
343 
344  vector_matrix_type m( nOrder+1 );
345  vector_matrix_type J( nOrder-1 );
346  vector_matrix_type dJ( nOrder-1 );
347 
348  m[0].resize( nOrder+1 ,__pts.size() );
349  m[0] = dpsi_1;
350  m[nOrder] = dpsi_1;
351 
352  vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
353 
354  vector_type tmp1 = ( ones - __pts )/value_type( 2.0 ); // (1-x)/2
355  vector_type tmp2 = ( ones + __pts )/value_type( 2.0 ); // (1+x)/2
356 
357  vector_type tmp_prod = ublas::element_prod( tmp1,tmp2 );
358 
359  vector_type tmp_i = tmp1;
360 
361  for ( uint16_type i= 1; i <= J.size() ; ++i )
362  {
363  m[i].resize( nOrder ,__pts.size() );
364  J[i-1] = JacobiBatchEvaluation<nOrder-1,value_type>( ( 2.0*i+1.0 ), M_b, __pts );
365  dJ[i-1] = JacobiBatchDerivation<nOrder-1,value_type>( ( 2.0*i+1.0 ), M_b, __pts );
366 
367  ublas::row( m[i],0 ) = ( - value_type( i ) - 1.0 ) * tmp_i / 2.0;
368 
369  for ( int16_type j=1; j < nOrder; ++j )
370  {
371 #if 0
372  vector_type tmp3(
373  ublas::element_prod( tmp1, ublas::row( J[i-1],j-1 ) )
374  + 2.0*element_prod( tmp_prod, ublas::row( dJ[i-1],j-1 ) )
375  - ( value_type( i ) + 1.0 ) * ublas::element_prod( tmp2,ublas::row( J[i-1],j-1 ) )
376  );
377  ublas::row( m[i],j ) = ( ublas::element_prod( tmp_i, tmp3 ) ) / value_type( 2.0 ) ;
378 #else
379  vector_type A( ublas::element_prod( ublas::row( m[i],0 ), tmp2 ) );
380  A = ublas::element_prod( A, ublas::row( J[i-1],j-1 ) );
381 
382  vector_type B( ublas::element_prod( 0.5*tmp_i , ublas::row( J[i-1],j-1 ) ) );
383 
384  vector_type C( ublas::element_prod( tmp_i , tmp_prod ) );
385  C = ublas::element_prod( C ,ublas::row( dJ[i-1],j-1 ) );
386 
387  ublas::row( m[i],j ) = A + B + C;
388 #endif
389  }
390 
391  tmp_i=ublas::element_prod( tmp_i,tmp1 );
392  }
393 
394  return m;
395 }
396 
397 
398 template<uint16_type Degree,
399  typename T,
400  template<class> class StoragePolicy>
401 typename Principal<Degree, T, StoragePolicy>::vector_vector_matrix_type
402 Principal<Degree, T, StoragePolicy>::derivate_3( vector_type const& __pts ) const
403 {
404  // std::cout <<"[principal]3D derivation ..."<< std::endl;
405  vector_vector_matrix_type m( nOrder+1 );
406 
407  vector_matrix_type dpsi_2( derivate_2( __pts ) );
408 
409  // i=0, 0 <= j <= nOrder , 0 <= k <= nOrder
410  m[0].resize( nOrder+1 );
411  m[0] = dpsi_2;
412 
413  // i=nOrder, 0 <= j <= nOrder , 0 <= k <= nOrder
414  m[nOrder].resize( nOrder+1 );
415  m[nOrder] = dpsi_2;
416 
417  vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
418  vector_type tmp1 = ( ones - __pts )/value_type( 2.0 ); // (1-x)/2
419  vector_type tmp2 = ( ones + __pts )/value_type( 2.0 ); // (1+x)/2
420  vector_type tmp_i = ones;
421  vector_type tmp_prod = ublas::element_prod( tmp1,tmp2 );
422 
423  vector_vector_matrix_type J( nOrder+1 );
424  vector_vector_matrix_type dJ( nOrder+1 );
425 
426  for ( int16_type i=1; i < nOrder; ++i )
427  {
428  m[i].resize( nOrder+1 );
429  J[i].resize( nOrder+1 );
430  dJ[i].resize( nOrder+1 );
431 
432  // 1 <= i <= nOrder-1, j = 0 , 0 <= k <= nOrder
433  m[i][0].resize( dpsi_2[i].size1(), dpsi_2[i].size2() );
434  m[i][0] = dpsi_2[i];
435 
436  // 1 <= i <= nOrder-1, j = nOrder , 0 <= k <= nOrder
437  m[i][nOrder].resize( dpsi_2[i].size1(),dpsi_2[i].size2() );
438  m[i][nOrder] = dpsi_2[i];
439 
440  tmp_i=ublas::element_prod( tmp_i,tmp1 ); // [(1-x) / 2]^i
441 
442  vector_type tmp_i_j( tmp_i );
443 
444  for ( int16_type j=1; j < nOrder; ++j )
445  {
446  // 1 <= i <= nOrder-1, 1 <= j <= nOrder-1 , 0 <= k <= nOrder-1
447  m[i][j].resize( nOrder,__pts.size() );
448  J[i][j] = JacobiBatchEvaluation<nOrder-1,value_type>( ( 2.0*i+2.0*j+1.0 ), 1.0, __pts );
449  dJ[i][j] = JacobiBatchDerivation<nOrder-1,value_type>( ( 2.0*i+2.0*j+1.0 ), 1.0, __pts );
450 
451  tmp_i_j=ublas::element_prod( tmp_i_j,tmp1 );
452 
453  ublas::row( m[i][j],0 ) = -value_type( i+j+1.0 )*tmp_i_j / value_type( 2.0 ); // -(i+j+1)/2*((1-x)/2)^(i+j)
454 
455  for ( int16_type k=1; k < nOrder; ++k ) // 1 <= k <= nOrder-1
456  {
457  vector_type tmp3( ublas::element_prod( tmp1, ublas::row( J[i][j],k-1 ) )
458  + 2.0*element_prod( tmp_prod, ublas::row( dJ[i][j],k-1 ) )
459  - ( value_type( i+j+1.0 ) ) * ublas::element_prod( tmp2,ublas::row( J[i][j],k-1 ) ) );
460 
461  ublas::row( m[i][j],k ) = ( ublas::element_prod( tmp_i_j, tmp3 ) ) / 2.0 ;
462  }
463 
464  }
465  }
466 
467  return m;
468 }
469 
470 
471 } /* Feel */
472 
473 #endif /* __Principal_H */

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