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>
51 template<
typename tag,
typename T>
61 typedef ublas::matrix<value_type,ublas::row_major> matrix_type;
63 static const uint16_type nDim = PTraits::nDim;
64 static const uint16_type nOrder = PTraits::nOrder;
69 typedef PTraits traits_type;
71 typedef typename traits_type::value_type value_type;
73 typedef typename traits_type::convex_type convex_type;
74 typedef typename traits_type::reference_convex_type reference_convex_type;
76 typedef typename traits_type::diff_pointset_type diff_pointset_type;
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;
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
101 template<
typename PrimalBasis>
156 static matrix_type
const&
d( uint16_type i )
167 static matrix_type
const&
derivate( uint16_type i )
180 template<
typename PrimalBasis>
181 static void initDerivation( PrimalBasis
const& basis );
189 static bool _S_has_derivation;
196 static std::vector<matrix_type> _S_D;
199 template<
typename Tag,
typename T>
200 bool Basis<Tag, T>::_S_has_derivation =
false;
202 template<
typename Tag,
typename T>
203 std::vector<typename Basis<Tag, T>::matrix_type> Basis<Tag, T>::_S_D;
207 #if defined( FEELPP_HAS_QD_REAL)
208 template<
typename PrimalBasis>
209 static void initDerivation( PrimalBasis
const& basis )
211 initDerivation( basis, mpl::bool_<boost::is_same<value_type,double>::value>() );
214 template<
typename PrimalBasis>
215 static void initDerivation( PrimalBasis
const& basis, mpl::bool_<false> )
218 if ( _S_has_derivation ==
false )
220 _S_has_derivation =
true;
222 reference_convex_type refconvex;
226 diff_pointset_type diff_pts( 1 );
228 matrix_type A( sub_type::evaluate( diff_pts.points() ) );
230 matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
232 LU<matrix_type> lu( A );
234 matrix_type C = lu.solve( D );
236 vector_matrix_type d ( sub_type::derivate( diff_pts.points() ) );
237 _S_D.resize( d.size() );
239 for (
size_type i = 0; i < d.size(); ++i )
241 _S_D[i] = ublas::prod( d[i], C );
242 glas::clean( _S_D[i] );
246 ublas::permutation_matrix<value_type> perm( A.size1() );
247 ublas::lu_factorize( A );
248 _S_D = sub_type::derivate( diff_pts.points() );
250 for (
size_type i = 0; i < _S_D.size(); ++i )
252 ublas::lu_substitute( ublas::matrix_expression<matrix_type>( _S_D[i] ), A );
253 glas::clean( _S_D[i] );
261 template<
typename PrimalBasis>
262 static void initDerivation( PrimalBasis
const& basis, mpl::bool_<true> )
264 if ( _S_has_derivation ==
false )
266 _S_has_derivation =
true;
268 typename qd_basis_type::reference_convex_type refconvex;
270 qd_diff_pointset_type diff_pts( 1 );
272 typename qd_basis_type::matrix_type A( qd_basis_type::evaluate( diff_pts.points() ) );
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 );
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;
281 _qd_derivatives_matrix.resize( d.size() );
282 _S_D.resize( _qd_derivatives_matrix.size() );
284 for (
size_type i = 0; i < d.size(); ++i )
286 _qd_derivatives_matrix[i] = ublas::prod( d[i], C );
287 glas::clean( _qd_derivatives_matrix[i] );
289 _S_D[i].resize( _qd_derivatives_matrix[i].size1(),_qd_derivatives_matrix[i].size2() );
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 ) );
295 glas::clean( _S_D[i] );
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() ) );
303 _S_D.resize( d.size() );
305 for (
size_type i = 0; i < d.size(); ++i )
307 ublas::lu_substitute( ublas::matrix_expression<matrix_type>( d[i] ), A );
309 _S_D[i].resize( d[i].size1(),d[i].size2() );
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 ) );
315 glas::clean( _S_D[i] );
323 template<
typename Tag,
typename T>
324 template<
typename PrimalBasis>
326 Basis<Tag, T>::initDerivation( PrimalBasis
const& basis )
328 typedef typename PrimalBasis::traits_type traits_type;
329 typedef typename traits_type::value_type value_type;
331 typedef typename traits_type::convex_type convex_type;
332 typedef typename traits_type::reference_convex_type reference_convex_type;
334 typedef typename traits_type::diff_pointset_type diff_pointset_type;
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;
343 if ( _S_has_derivation ==
false )
345 _S_has_derivation =
true;
347 reference_convex_type refconvex;
350 diff_pointset_type diff_pts( 1 );
351 matrix_type A( basis.evaluate( diff_pts.points() ) );
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 );
358 vector_matrix_type d ( basis.derivate( diff_pts.points() ) );
359 _S_D.resize( d.size() );
361 for (
size_type i = 0; i < d.size(); ++i )
363 _S_D[i] = ublas::prod( d[i], C );
364 glas::clean( _S_D[i] );
368 ublas::permutation_matrix<double> perm( A.size1() );
369 ublas::lu_factorize( A );
370 _S_D = basis.derivate( diff_pts.points() );
372 for (
size_type i = 0; i < d.size(); ++i )
374 ublas::lu_substitute( _S_D[i], A );
375 glas::clean( _S_D[i] );
381 #endif // FEELPP_HAS_QD_REAL
Basis(PrimalBasis const &p)
Definition: basis.hpp:102
static matrix_type const & d(uint16_type i)
derivatives of Dubiner polynomials the derivatives are computed at the nodes of the lattice ...
Definition: basis.hpp:156
virtual ~Basis()
Definition: basis.hpp:117
static matrix_type const & derivate(uint16_type i)
derivatives of Dubiner polynomials the derivatives are computed at the nodes of the lattice ...
Definition: basis.hpp:167
Basis(Basis const &b)
Definition: basis.hpp:111
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Base class for basis.
Definition: basis.hpp:52