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
operations.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: 2005-11-24
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 __operations_H
30 #define __operations_H 1
31 
32 
33 #include <feel/feelalg/svd.hpp>
34 #include <feel/feelpoly/im.hpp>
36 
37 namespace Feel
38 {
39 
46 template<
47 typename Poly,
48  template<uint16_type> class PsetType
49  >
50 typename Poly::value_type
52 {
53  return p.coeff()( 0, 0 );
54 }
55 
61 template<typename Poly, template<uint16_type> class Type>
62 Polynomial<Poly, Type>
64 {
65  return Polynomial<Poly, Type>( ublas::prod( p.coeff(), p.basis().d( 0 ) ), true );
66 }
72 template<typename Poly, template<uint16_type> class Type>
73 PolynomialSet<Poly, Type>
75 {
76  return PolynomialSet<Poly, Type>( ublas::prod( p.coeff(), p.basis().d( 0 ) ), true );
77 }
83 template<typename Poly, template<uint16_type> class Type>
84 Polynomial<Poly, Type>
86 {
87  BOOST_STATIC_ASSERT( ( Polynomial<Poly,Type>::nDim >= 2 ) );
88  return Polynomial<Poly, Type>( Poly(), ublas::prod( p.coeff(), p.basis().d( 1 ) ), true );
89 }
95 template<typename Poly, template<uint16_type> class Type>
96 PolynomialSet<Poly, Type>
98 {
99  BOOST_STATIC_ASSERT( ( Polynomial<Poly,Type>::nDim >= 2 ) );
100  return PolynomialSet<Poly, Type>( Poly(), ublas::prod( p.coeff(), p.basis().d( 1 ) ), true );
101 }
107 template<typename Poly, template<uint16_type> class Type>
108 Polynomial<Poly, Type>
110 {
111  BOOST_STATIC_ASSERT( ( Polynomial<Poly,Type>::nDim >= 3 ) );
112  return Polynomial<Poly, Type>( Poly(), ublas::prod( p.coeff(), p.basis().d( 2 ) ), true );
113 }
119 template<typename Poly, template<uint16_type> class Type>
120 PolynomialSet<Poly, Type>
122 {
123  BOOST_STATIC_ASSERT( ( PolynomialSet<Poly,Type>::nDim >= 3 ) );
124  return PolynomialSet<Poly, Type>( Poly(), ublas::prod( p.coeff(), p.basis().d( 2 ) ), true );
125 }
131 template<typename Poly>
132 Polynomial<Poly, Vectorial>
134 {
135  const int nDim = Polynomial<Poly, Scalar>::nDim;
136 
137  typedef typename Poly::value_type value_type;
138  const int ndof = p.coeff().size2();
139  ublas::matrix<value_type> c ( nDim*nDim, ndof );
140  c.clear();
141 
142  for ( int i = 0; i <nDim; ++i )
143  ublas::row( c, nDim*i+i ) = ublas::row( ublas::prod( p.coeff(), p.basis().d( i ) ), 0 );
144 
145  return Polynomial<Poly, Vectorial>( Poly(), c, true );
146 }
152 template<typename Poly>
153 Polynomial<Poly, Tensor2>
155 {
156  const int nDim = Polynomial<Poly, Vectorial>::nDim;
157  const int nComponents = Polynomial<Poly, Tensor2>::nComponents;
158 
159  typedef typename Poly::value_type value_type;
160  const int ndof = p.coeff().size2();
161  ublas::matrix<value_type> c ( nComponents, ndof );
162  c.clear();
163 
164  for ( int i = 0; i < nDim; ++i )
165  for ( int j = 0; j < nDim; ++j )
166  ublas::row( c, nComponents*i+nDim*j ) = ublas::row( ublas::prod( p[j].coeff(), p.basis().d( i ) ), 0 );
167 
168  return Polynomial<Poly, Tensor2>( Poly(), c, true );
169 }
175 template<typename Poly>
176 PolynomialSet<Poly, Vectorial>
178 {
179  const int nDim = PolynomialSet<Poly, Scalar>::nDim;
180 
181  typedef typename Poly::value_type value_type;
182  const int n1 = p.coeff().size1();
183  const int n2 = p.coeff().size2();
184  ublas::matrix<value_type> c ( nDim*nDim*n1, n2 );
185  c.clear();
186 
187  for ( int i = 0; i <nDim; ++i )
188  {
189  ublas::project( c,
190  ublas::slice( nDim*n1*i+i, nDim, n1 ),
191  ublas::slice( 0, 1, n2 ) ) = ublas::prod( p.coeff(), p.basis().d( i ) );
192  }
193 
194  return PolynomialSet<Poly, Vectorial>( Poly(), c, true );
195 }
201 template<typename Poly>
202 PolynomialSet<Poly, Tensor2>
204 {
205  const int nDim = PolynomialSet<Poly, Scalar>::nDim;
206  const int nComponents = PolynomialSet<Poly, Tensor2>::nComponents;
207 
208  typedef typename Poly::value_type value_type;
209  const int n1 = p.coeff().size1();
210  const int n2 = p.coeff().size2();
211  ublas::matrix<value_type> c ( nComponents*n1, n2 );
212  c.clear();
213 
214  for ( int i = 0; i <nDim; ++i )
215  for ( int j = 0; j < nDim; ++j )
216  {
217  ublas::project( c,
218  ublas::slice( nDim*n1*i+i, nDim, n1 ),
219  ublas::slice( 0, 1, n2 ) ) = ublas::prod( p.coeff(), p.basis().d( i ) );
220  }
221 
222  return PolynomialSet<Poly, Tensor2>( Poly(), c, true );
223 }
229 template<typename Poly>
230 PolynomialSet<Poly, Tensor2>
232 {
233  const int nDim = PolynomialSet<Poly, Scalar>::nDim;
234  const int nComponents = PolynomialSet<Poly, Tensor2>::nComponents;
235 
236  typedef typename Poly::value_type value_type;
237  const int n1 = p.coeff().size1();
238  const int n2 = p.coeff().size2();
239  ublas::matrix<value_type> c ( nComponents*nComponents*n1, n2 );
240  c.clear();
241 
242  for ( int i = 0; i <nDim; ++i )
243  for ( int j = 0; j < nDim; ++j )
244  {
245  ublas::project( c,
246  ublas::slice( nComponents*n1*( nDim*j+i ), nComponents, n1 ),
247  ublas::slice( 0, 1, n2 ) ) = ublas::prod( p.coeff(), p.d( i, j ) );
248  }
249 
250  return PolynomialSet<Poly, Tensor2>( Poly(), c, true );
251 }
252 
258 template<typename Poly>
259 Polynomial<Poly, Scalar>
261 {
262  const int nDim = Polynomial<Poly, Vectorial>::nDim;
263 
264  typedef typename Poly::value_type value_type;
265 
266  ublas::matrix<value_type> c ( ublas::zero_matrix<value_type>( 1, p.coeff().size2() ) );
267 
268  for ( int i = 0; i <nDim; ++i )
269  {
270  ublas::noalias( c ) += ublas::prod( p[i].coeff(), p.basis().d( i ) );
271  }
272 
273  return Polynomial<Poly, Scalar>( Poly(), c, true );
274 }
280 template<typename Poly>
281 PolynomialSet<Poly, Scalar>
283 {
284  const int nComponents = PolynomialSet<Poly, Vectorial>::nComponents;
285 
286  typedef typename Poly::value_type value_type;
287 
288  ublas::matrix<value_type> c ( ublas::zero_matrix<value_type>( p.coeff().size1()/nComponents,
289  p.coeff().size2() ) );
290 
291  for ( int i = 0; i <nComponents; ++i )
292  {
293  ublas::noalias( c ) += ublas::prod( p[i].coeff(), p.basis().d( i ) );
294  }
295 
296  return PolynomialSet<Poly, Scalar>( Poly(), c, true );
297 }
303 template<typename Poly>
304 Polynomial<Poly, Vectorial>
306 {
307  const int nDim = Polynomial<Poly, Vectorial>::nDim;
308  const int nComponents = Polynomial<Poly, Vectorial>::nComponents;
309  BOOST_STATIC_ASSERT( ( Polynomial<Poly, Vectorial>::nDim == 3 ) );
310 
311  typedef typename Poly::value_type value_type;
312 
313  ublas::matrix<value_type> c ( p.coeff().size1(), p.coeff().size2() );
314  c.clear();
315  const uint16_type dim_p = p.coeff().size1()/nComponents;
316  const uint16_type ndof = p.coeff().size2();
317 
318  ublas::project( c,
319  ublas::slice( 0, nComponents, dim_p/nComponents ),
320  ublas::slice( 0, 1, ndof ) ) = p[2].derivate( 1 ).coeff() - p[1].derivate( 2 ).coeff();
321  ublas::project( c,
322  ublas::slice( dim_p+1, nComponents, dim_p/nComponents ),
323  ublas::slice( 0, 1, ndof ) ) = p[0].derivate( 2 ).coeff() - p[2].derivate( 0 ).coeff();
324  ublas::project( c,
325  ublas::slice( 2*dim_p+2, nComponents, dim_p/nComponents ),
326  ublas::slice( 0, 1, ndof ) ) = p[1].derivate( 0 ).coeff() - p[0].derivate( 1 ).coeff();
327 
328  return Polynomial<Poly, Vectorial>( Poly(), c, true );
329 }
335 template<typename Poly>
336 PolynomialSet<Poly, Vectorial>
338 {
339  const int nDim = PolynomialSet<Poly, Vectorial>::nDim;
340  const int nComponents = PolynomialSet<Poly, Vectorial>::nComponents;
341  BOOST_STATIC_ASSERT( ( PolynomialSet<Poly, Vectorial>::nDim == 3 ) );
342 
343  typedef typename Poly::value_type value_type;
344 
345  ublas::matrix<value_type> c ( p.coeff().size1(), p.coeff().size2() );
346  c.clear();
347  const uint16_type dim_p = p.coeff().size1()/nComponents;
348  const uint16_type ndof = p.coeff().size2();
349 
350  ublas::project( c,
351  ublas::slice( 0, nComponents, dim_p/nComponents ),
352  ublas::slice( 0, 1, ndof ) ) = p[2].derivate( 1 ).coeff() - p[1].derivate( 2 ).coeff();
353  ublas::project( c,
354  ublas::slice( dim_p+1, nComponents, dim_p/nComponents ),
355  ublas::slice( 0, 1, ndof ) ) = p[0].derivate( 2 ).coeff() - p[2].derivate( 0 ).coeff();
356  ublas::project( c,
357  ublas::slice( 2*dim_p+2, nComponents, dim_p/nComponents ),
358  ublas::slice( 0, 1, ndof ) ) = p[1].derivate( 0 ).coeff() - p[0].derivate( 1 ).coeff();
359 
360  return PolynomialSet<Poly, Vectorial>( Poly(), c, true );
361 }
367 template<typename Poly>
368 Polynomial<Poly, Vectorial>
370 {
371  const int nDim = Polynomial<Poly, Vectorial>::nDim;
372  const int nComponents = Polynomial<Poly, Vectorial>::nComponents;
373  BOOST_STATIC_ASSERT( ( Polynomial<Poly, Vectorial>::nDim == 3 ) );
374 
375  typedef typename Poly::value_type value_type;
376 
377  ublas::matrix<value_type> c ( p.coeff().size1(), p.coeff().size2() );
378  const uint16_type dim_p = p.coeff().size1()/nComponents;
379  const uint16_type ndof = p.coeff().size2();
380 
381  ublas::project( c,
382  ublas::slice( 0, nComponents, dim_p/nComponents ),
383  ublas::slice( 0, 1, ndof ) ) = - p[2].derivate( 1 ).coeff() + p[1].derivate( 2 ).coeff();
384  ublas::project( c,
385  ublas::slice( dim_p+1, nComponents, dim_p/nComponents ),
386  ublas::slice( 0, 1, ndof ) ) = - p[0].derivate( 2 ).coeff() + p[2].derivate( 0 ).coeff();
387  ublas::project( c,
388  ublas::slice( 2*dim_p+2, nComponents, dim_p/nComponents ),
389  ublas::slice( 0, 1, ndof ) ) = - p[1].derivate( 0 ).coeff() + p[0].derivate( 1 ).coeff();
390 
391  return Polynomial<Poly, Vectorial>( Poly(), c, true );
392 }
398 template<typename Poly>
399 PolynomialSet<Poly, Vectorial>
401 {
402  const int nDim = PolynomialSet<Poly, Vectorial>::nDim;
403  const int nComponents = PolynomialSet<Poly, Vectorial>::nComponents;
404  BOOST_STATIC_ASSERT( ( PolynomialSet<Poly, Vectorial>::nDim == 3 ) );
405 
406  typedef typename Poly::value_type value_type;
407 
408  ublas::matrix<value_type> c ( p.coeff().size1(), p.coeff().size2() );
409  const uint16_type dim_p = p.coeff().size1()/nComponents;
410  const uint16_type ndof = p.coeff().size2();
411 
412  ublas::project( c,
413  ublas::slice( 0, nComponents, dim_p/nComponents ),
414  ublas::slice( 0, 1, ndof ) ) = - p[2].derivate( 1 ).coeff() + p[1].derivate( 2 ).coeff();
415  ublas::project( c,
416  ublas::slice( dim_p+1, nComponents, dim_p/nComponents ),
417  ublas::slice( 0, 1, ndof ) ) = - p[0].derivate( 2 ).coeff() + p[2].derivate( 0 ).coeff();
418  ublas::project( c,
419  ublas::slice( 2*dim_p+2, nComponents, dim_p/nComponents ),
420  ublas::slice( 0, 1, ndof ) ) = - p[1].derivate( 0 ).coeff() + p[0].derivate( 1 ).coeff();
421 
422  return PolynomialSet<Poly, Vectorial>( Poly(), c, true );
423 }
432 template<typename Pset,
433  typename Func,
434  typename IM>
435 Polynomial<Pset, Scalar>
436 project( Pset const& pset, Func const& f, IM const& im )
437 {
438  typedef typename Func::value_type value_type;
439 
440  // evaluate f at quad nodes
441  ublas::vector<value_type> fq( f( im.points() ) );
442 
443  // evaluate pset at quad nodes
444  ublas::matrix<value_type> psetq( pset.evaluate( im.points() ) );
445 
446  ublas::matrix<value_type> c( 1, psetq.size1() );
447  ublas::row( c, 0 ) = ublas::prod( psetq, ublas::element_prod( fq, im.weights() ) );
448  return Polynomial<Pset, Scalar>( pset, c, true );
449 }
450 
457 template<typename P, template<uint16_type> class Type,
458  template<class, template<uint16_type> class> class Poly1,
459  template<class, template<uint16_type> class> class Poly2 >
460 PolynomialSet<P, Type>
461 unite( Poly1<P, Type> const& pset1,
462  Poly2<P, Type> const& pset2 )
463 {
464  typedef PolynomialSet<P, Type> res_type;
465  typedef typename res_type::value_type value_type;
466 
467  FEELPP_ASSERT( pset1.coeff().size2() == pset2.coeff().size2() )
468  ( pset1.coeff().size2() )( pset2.coeff().size2() ).error( "incompatible size" );
469 
470  ublas::matrix<value_type> M( pset1.coeff().size1()+pset2.coeff().size1(), pset1.coeff().size2() );
471  project( M,
472  ublas::range( 0, pset1.coeff().size1() ),
473  ublas::range( 0,pset1.coeff().size2() ) ) = pset1.coeff();
474  project( M,
475  ublas::range( pset1.coeff().size1(), pset1.coeff().size1()+pset2.coeff().size1() ),
476  ublas::range( 0,pset1.coeff().size2() ) ) = pset2.coeff();
477 
478  M = res_type::polyset_type::toMatrix( M );
479 #if 0
480  std::cout << "before SVD M = " << M << "\n";
481  Eigen::MatrixXd A ( Eigen::MatrixXd::Zero( M.size1(), M.size2() ) );
482 
483  for ( int i = 0; i < M.size1(); ++i )
484  for ( int j = 0; j < M.size2(); ++j )
485  {
486  A( i,j ) = M( i,j );
487  }
488 
489  Eigen::SVD<Eigen::MatrixXd> svdOfA( A );
490  std::cout << "SVDofA.S() = " << svdOfA.singularValues() << "\n";
491  std::cout << "SVDofA.V() = " << svdOfA.matrixV().transpose() << "\n";
492 #endif
494 #if 0
495  //std::cout << "S() = " << svd.S() << "\n";
496  ublas::matrix<value_type> m ( svdOfA.singularValues().size(), M.size2() );
497 
498  for ( int i = 0; i < m.size1(); ++i )
499  for ( int j = 0; j < m.size2(); ++j )
500  {
501  m( i,j ) = svdOfA.matrixV()( j, i );
502  }
503 
504 #else
505  ublas::matrix<value_type> m ( ublas::subrange( svd.V(), 0, svd.S().size(), 0, M.size2() ) );
506 #endif
507  return res_type( P(), res_type::polyset_type::toType( m ), true );
508 }
509 
510 }
511 
512 #endif /* __operations_H */

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