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
basis.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-09
7 
8  Copyright (C) 2006 EPFL
9  Copyright (C) 2007 Université Joseph Fourier Grenoble 1
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __Basis_H
31 #define __Basis_H 1
32 
33 #include <boost/numeric/ublas/vector.hpp>
34 #include <boost/numeric/ublas/vector_proxy.hpp>
35 #include <boost/numeric/ublas/fwd.hpp>
36 #include <boost/numeric/ublas/triangular.hpp>
37 #include <boost/numeric/ublas/matrix.hpp>
38 #include <boost/numeric/ublas/operation.hpp>
39 #include <boost/numeric/ublas/lu.hpp>
40 
41 namespace Feel
42 {
51 template<typename tag, typename T>
52 class Basis
53 {
54 public:
55 
56 
60  typedef T value_type;
61  typedef ublas::matrix<value_type,ublas::row_major> matrix_type;
62 #if 0
63  static const uint16_type nDim = PTraits::nDim;
64  static const uint16_type nOrder = PTraits::nOrder;
65 
66 
67 
68  typedef P sub_type;
69  typedef PTraits traits_type;
70 
71  typedef typename traits_type::value_type value_type;
72 
73  typedef typename traits_type::convex_type convex_type;
74  typedef typename traits_type::reference_convex_type reference_convex_type;
75 
76  typedef typename traits_type::diff_pointset_type diff_pointset_type;
77 
78  typedef typename traits_type::storage_policy storage_policy;
79  typedef typename traits_type::matrix_type matrix_type;
80  typedef typename traits_type::vector_matrix_type vector_matrix_type;
81  typedef typename traits_type::matrix_node_type matrix_node_type;
82  typedef typename traits_type::points_type points_type;
83  typedef typename traits_type::node_type node_type;
84 
85 
86 #if defined( FEELPP_HAS_QD_REAL )
87  typedef typename traits_type::template ChangeValueType<qd_real>::type qd_basis_type;
88  typedef typename traits_type::template ChangeValueType<qd_real>::traits_type::diff_pointset_type qd_diff_pointset_type;
89 #endif // FEELPP_HAS_QD_REAL
90 #endif
91 
92 
96 
101  template<typename PrimalBasis>
102  Basis( PrimalBasis const& p )
103  {
104  initDerivation( p );
105  }
106 
111  Basis( Basis const & b )
112  {}
113 
117  virtual ~Basis()
118  {}
119 
120 
122 
126 
127 
129 
133 
134 
136 
140 
141 
143 
147 
148 
149 
156  static matrix_type const& d( uint16_type i )
157  {
158  return _S_D[i];
159  }
160 
167  static matrix_type const& derivate( uint16_type i )
168  {
169  return _S_D[i];
170  }
171 
172 
173 
175 
176 
177 
178 protected:
179 
180  template<typename PrimalBasis>
181  static void initDerivation( PrimalBasis const& basis );
182 
183 private:
184 
189  static bool _S_has_derivation;
190 
195  //static std::vector<matrix_type> _S_D;
196  static std::vector<matrix_type> _S_D;
197 };
198 
199 template<typename Tag, typename T>
200 bool Basis<Tag, T>::_S_has_derivation = false;
201 
202 template<typename Tag, typename T>
203 std::vector<typename Basis<Tag, T>::matrix_type> Basis<Tag, T>::_S_D;
204 
205 
206 
207 #if defined( FEELPP_HAS_QD_REAL)
208 template<typename PrimalBasis>
209 static void initDerivation( PrimalBasis const& basis )
210 {
211  initDerivation( basis, mpl::bool_<boost::is_same<value_type,double>::value>() );
212 }
213 
214 template<typename PrimalBasis>
215 static void initDerivation( PrimalBasis const& basis, mpl::bool_<false> )
216 {
217 
218  if ( _S_has_derivation == false )
219  {
220  _S_has_derivation = true;
221 
222  reference_convex_type refconvex;
223 
224  // constructor pointset for differentiation only in
225  // the interior(1)
226  diff_pointset_type diff_pts( 1 );
227 
228  matrix_type A( sub_type::evaluate( diff_pts.points() ) );
229 #if 1
230  matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
231 
232  LU<matrix_type> lu( A );
233 
234  matrix_type C = lu.solve( D );
235 
236  vector_matrix_type d ( sub_type::derivate( diff_pts.points() ) );
237  _S_D.resize( d.size() );
238 
239  for ( size_type i = 0; i < d.size(); ++i )
240  {
241  _S_D[i] = ublas::prod( d[i], C );
242  glas::clean( _S_D[i] );
243  }
244 
245 #else
246  ublas::permutation_matrix<value_type> perm( A.size1() );
247  ublas::lu_factorize( A );
248  _S_D = sub_type::derivate( diff_pts.points() );
249 
250  for ( size_type i = 0; i < _S_D.size(); ++i )
251  {
252  ublas::lu_substitute( ublas::matrix_expression<matrix_type>( _S_D[i] ), A );
253  glas::clean( _S_D[i] );
254  }
255 
256 #endif
257 
258  }
259 }
260 
261 template<typename PrimalBasis>
262 static void initDerivation( PrimalBasis const& basis, mpl::bool_<true> )
263 {
264  if ( _S_has_derivation == false )
265  {
266  _S_has_derivation = true;
267 
268  typename qd_basis_type::reference_convex_type refconvex;
269 
270  qd_diff_pointset_type diff_pts( 1 );
271 
272  typename qd_basis_type::matrix_type A( qd_basis_type::evaluate( diff_pts.points() ) );
273 #if 1
274  typename qd_basis_type::matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
275  LU<typename qd_basis_type::matrix_type> lu( A );
276  typename qd_basis_type::matrix_type C = lu.solve( D );
277 
278  typename qd_basis_type::vector_matrix_type d ( qd_basis_type::derivate( diff_pts.points() ) );
279  std::vector<typename qd_basis_type::matrix_type> _qd_derivatives_matrix;
280 
281  _qd_derivatives_matrix.resize( d.size() );
282  _S_D.resize( _qd_derivatives_matrix.size() );
283 
284  for ( size_type i = 0; i < d.size(); ++i )
285  {
286  _qd_derivatives_matrix[i] = ublas::prod( d[i], C );
287  glas::clean( _qd_derivatives_matrix[i] );
288 
289  _S_D[i].resize( _qd_derivatives_matrix[i].size1(),_qd_derivatives_matrix[i].size2() );
290 
291  for ( size_type j = 0; j < _S_D[i].size1(); ++j )
292  for ( size_type k = 0; k < _S_D[i].size2(); ++k )
293  _S_D[i]( j,k ) = value_type( _qd_derivatives_matrix[i]( j,k ) );
294 
295  glas::clean( _S_D[i] );
296  }
297 
298 #else
299  ublas::permutation_matrix<typename qd_basis_type::value_type> perm( A.size1() );
300  ublas::lu_factorize( A );
301  typename qd_basis_type::vector_matrix_type d ( qd_basis_type::derivate( diff_pts.points() ) );
302 
303  _S_D.resize( d.size() );
304 
305  for ( size_type i = 0; i < d.size(); ++i )
306  {
307  ublas::lu_substitute( ublas::matrix_expression<matrix_type>( d[i] ), A );
308 
309  _S_D[i].resize( d[i].size1(),d[i].size2() );
310 
311  for ( size_type j = 0; j < _S_D[i].size1(); ++j )
312  for ( size_type k = 0; k < _S_D[i].size2(); ++k )
313  _S_D[i]( j,k ) = value_type( d[i]( j,k ) );
314 
315  glas::clean( _S_D[i] );
316  }
317 
318 #endif
319  }
320 }
321 
322 #else
323 template<typename Tag, typename T>
324 template<typename PrimalBasis>
325 void
326 Basis<Tag, T>::initDerivation( PrimalBasis const& basis )
327 {
328  typedef typename PrimalBasis::traits_type traits_type;
329  typedef typename traits_type::value_type value_type;
330 
331  typedef typename traits_type::convex_type convex_type;
332  typedef typename traits_type::reference_convex_type reference_convex_type;
333 
334  typedef typename traits_type::diff_pointset_type diff_pointset_type;
335 
336  typedef typename traits_type::storage_policy storage_policy;
337  typedef typename traits_type::matrix_type matrix_type;
338  typedef typename traits_type::vector_matrix_type vector_matrix_type;
339  typedef typename traits_type::matrix_node_type matrix_node_type;
340  typedef typename traits_type::points_type points_type;
341  typedef typename traits_type::node_type node_type;
342 
343  if ( _S_has_derivation == false )
344  {
345  _S_has_derivation = true;
346 
347  reference_convex_type refconvex;
348  // constructor pointset for differentiation only in
349  // the interior(1)
350  diff_pointset_type diff_pts( 1 );
351  matrix_type A( basis.evaluate( diff_pts.points() ) );
352 
353 #if 1
354  matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
355  LU<matrix_type> lu( A );
356  matrix_type C = lu.solve( D );
357 
358  vector_matrix_type d ( basis.derivate( diff_pts.points() ) );
359  _S_D.resize( d.size() );
360 
361  for ( size_type i = 0; i < d.size(); ++i )
362  {
363  _S_D[i] = ublas::prod( d[i], C );
364  glas::clean( _S_D[i] );
365  }
366 
367 #else
368  ublas::permutation_matrix<double> perm( A.size1() );
369  ublas::lu_factorize( A );
370  _S_D = basis.derivate( diff_pts.points() );
371 
372  for ( size_type i = 0; i < d.size(); ++i )
373  {
374  ublas::lu_substitute( _S_D[i], A );
375  glas::clean( _S_D[i] );
376  }
377 
378 #endif
379  }
380 }
381 #endif // FEELPP_HAS_QD_REAL
382 
383 
384 } // Feel
385 #endif /* __Basis_H */

Generated on Sun Oct 20 2013 08:24:54 for Feel++ by doxygen 1.8.4