37 #include <boost/shared_array.hpp>
38 #include <boost/lambda/lambda.hpp>
39 #include <boost/lambda/bind.hpp>
40 #include <boost/utility.hpp>
42 #include <boost/numeric/ublas/vector.hpp>
43 #include <boost/numeric/ublas/vector_proxy.hpp>
45 #include <feel/feelcore/feel.hpp>
50 namespace ublas = boost::numeric::ublas;
52 enum BDFTimeScheme { BDF_ORDER_ONE=1, BDF_ORDER_TWO, BDF_ORDER_THREE, BDF_ORDER_FOUR, BDF_MAX_ORDER = 4 };
90 template<
typename SpaceType>
95 typedef SpaceType space_type;
96 typedef boost::shared_ptr<space_type> space_ptrtype;
97 typedef typename space_type::element_type element_type;
98 typedef typename space_type::return_type return_type;
99 typedef typename element_type::value_type value_type;
100 typedef boost::shared_ptr<element_type> unknown_type;
101 typedef std::vector< unknown_type > unknowns_type;
110 Bdf( space_ptrtype
const& space );
143 template<
typename container_type>
144 void shiftRight(
typename space_type::template Element<value_type, container_type>
const& u_curr );
148 element_type
derivate(
int, scalar_type dt )
const;
153 element_type
derivate(
int n )
const;
166 unknowns_type
const&
unknowns()
const;
169 element_type&
unknown(
int i );
171 template<
typename container_type>
172 void setUnknown(
int i,
typename space_type::template Element<value_type, container_type>
const& e )
174 M_unknowns[i]->assign( e );
177 void showMe( std::ostream& __out = std::cout )
const;
182 space_ptrtype M_space;
185 BDFTimeScheme M_order;
191 std::vector<ublas::vector<double> > M_alpha;
194 std::vector<ublas::vector<double> > M_beta;
197 unknowns_type M_unknowns;
200 template <
typename SpaceType>
205 M_alpha( BDF_MAX_ORDER ),
206 M_beta( BDF_MAX_ORDER )
208 for (
int i = 0; i < BDF_MAX_ORDER; ++i )
210 M_alpha[ i ].resize( i+2 );
211 M_beta[ i ].resize( i+1 );
214 for (
int i = 0; i < BDF_MAX_ORDER; ++i )
218 M_alpha[i][ 0 ] = 1.;
219 M_alpha[i][ 1 ] = 1.;
225 M_alpha[i][ 0 ] = 3. / 2.;
226 M_alpha[i][ 1 ] = 2.;
227 M_alpha[i][ 2 ] = -1. / 2.;
229 M_beta[i][ 1 ] = -1.;
234 M_alpha[i][ 0 ] = 11. / 6.;
235 M_alpha[i][ 1 ] = 3.;
236 M_alpha[i][ 2 ] = -3. / 2.;
237 M_alpha[i][ 3 ] = 1. / 3.;
239 M_beta[i][ 1 ] = -3.;
245 M_alpha[i][ 0 ] = 25. / 12.;
246 M_alpha[i][ 1 ] = 4.;
247 M_alpha[i][ 2 ] = -3.;
248 M_alpha[i][ 3 ] = 4. / 3.;
249 M_alpha[i][ 4 ] = -1. / 4.;
251 M_beta[i][ 1 ] = -6.;
253 M_beta[i][ 3 ] = -1.;
257 M_unknowns.resize( BDF_MAX_ORDER );
259 for ( uint8_type __i = 0; __i < ( uint8_type )BDF_MAX_ORDER; ++__i )
261 M_unknowns[__i] = unknown_type(
new element_type( M_space ) );
262 M_unknowns[__i]->zero();
267 template <
typename SpaceType>
272 template <
typename SpaceType>
276 std::for_each( M_unknowns.begin(), M_unknowns.end(), *boost::lambda::_1 = u0 );
279 template <
typename SpaceType>
286 std::copy( uv0.begin(), uv0.end(), M_unknowns.begin() );
289 template <
typename SpaceType>
291 typename Bdf<SpaceType>::unknowns_type&
297 template <
typename SpaceType>
298 typename Bdf<SpaceType>::element_type&
302 return *M_unknowns[i];
308 template <
typename SpaceType>
312 FEELPP_ASSERT( i <
size_type( n + 1 ) ).error(
"Error in specification of the time derivative coefficient for the BDF formula (out of range error)" );
313 return M_alpha[n-1][ i ]/dt;
316 template <
typename SpaceType>
320 FEELPP_ASSERT( i < n ).error(
"Error in specification of the time derivative coefficient for the BDF formula (out of range error)" );
321 return M_beta[n-1][ i ];
324 template <
typename SpaceType>
326 Bdf<SpaceType>::showMe( std::ostream& __out )
const
329 __out <<
"*** BDF Time discretization of order " << M_order <<
" ***"
330 <<
" size : " << M_unknowns[0]->size() <<
"\n"
331 <<
" alpha : " << M_alpha <<
"\n"
332 <<
" beta : " << M_beta <<
"\n";
336 template <
typename SpaceType>
337 template<
typename container_type>
341 using namespace boost::lambda;
342 typename unknowns_type::reverse_iterator __it = boost::next( M_unknowns.rbegin() );
343 std::for_each( M_unknowns.rbegin(), boost::prior( M_unknowns.rend() ),
344 ( *lambda::_1 = *( *lambda::var( __it ) ), ++lambda::var( __it ) ) );
346 *M_unknowns[0] = __new_unk;
357 template <
typename SpaceType>
358 typename Bdf<SpaceType>::element_type
361 element_type __t( M_space );
364 __t.add( ( 1./dt ), derivate( n ) );
369 template <
typename SpaceType>
370 typename Bdf<SpaceType>::element_type
373 element_type __t( M_space );
376 FEELPP_ASSERT( __t.size() == M_space->nDof() )( __t.size() )( M_space->nDof() ).error(
"invalid space element size" );
377 FEELPP_ASSERT( __t.size() == M_unknowns[0]->size() )( __t.size() )( M_unknowns[0]->size() ).error(
"invalid space element size" );
379 for ( uint8_type i = 0; i < n; ++i )
380 __t.add( M_alpha[n-1][ i+1 ], *M_unknowns[i] );
385 template <
typename SpaceType>
386 typename Bdf<SpaceType>::element_type
389 element_type __t( M_space );
392 FEELPP_ASSERT( __t.size() == M_space->nDof() )( __t.size() )( M_space->nDof() ).error(
"invalid space element size" );
393 FEELPP_ASSERT( __t.size() == M_unknowns[0]->size() )( __t.size() )( M_unknowns[0]->size() ).error(
"invalid space element size" );
395 for ( uint8_type i = 0; i < n; ++i )
396 __t.add( M_beta[n-1][ i ], *M_unknowns[ i ] );