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
bdf.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-12-30
7 
8  Copyright (C) 2006 Universit� Joseph Fourier (Grenoble)
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 _BDF_H
30 #define _BDF_H
31 
32 #include <string>
33 #include <iostream>
34 #include <sstream>
35 #include <algorithm>
36 
37 #include <boost/shared_array.hpp>
38 #include <boost/lambda/lambda.hpp>
39 #include <boost/lambda/bind.hpp>
40 #include <boost/utility.hpp>
41 
42 #include <boost/numeric/ublas/vector.hpp>
43 #include <boost/numeric/ublas/vector_proxy.hpp>
44 
45 #include <feel/feelcore/feel.hpp>
46 #include <feel/feelalg/glas.hpp>
47 
48 namespace Feel
49 {
50 namespace ublas = boost::numeric::ublas;
51 
52 enum BDFTimeScheme { BDF_ORDER_ONE=1, BDF_ORDER_TWO, BDF_ORDER_THREE, BDF_ORDER_FOUR, BDF_MAX_ORDER = 4 };
53 
90 template<typename SpaceType>
91 class Bdf
92 {
93 public:
94 
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;
102  typedef typename node<value_type>::type node_type;
103 
110  Bdf( space_ptrtype const& space );
111 
112  ~Bdf();
113 
114 
121  BDFTimeScheme order() const
122  {
123  return M_order;
124  }
125 
130  void initialize( element_type const& u0 );
131 
136  void initialize( unknowns_type const& uv0 );
137 
143  template<typename container_type>
144  void shiftRight( typename space_type::template Element<value_type, container_type> const& u_curr );
145 
148  element_type derivate( int, scalar_type dt ) const;
149 
153  element_type derivate( int n ) const;
154 
157  element_type extrapolate( int n ) const;
158 
160  double derivateCoefficient( int n, size_type i, double dt = 1.0 ) const;
161 
163  double extrapolateCoefficient( int n, size_type i, double dt = 1.0 ) const;
164 
166  unknowns_type const& unknowns() const;
167 
169  element_type& unknown( int i );
170 
171  template<typename container_type>
172  void setUnknown( int i, typename space_type::template Element<value_type, container_type> const& e )
173  {
174  M_unknowns[i]->assign( e );
175  }
176 
177  void showMe( std::ostream& __out = std::cout ) const;
178 
179 private:
180 
182  space_ptrtype M_space;
183 
185  BDFTimeScheme M_order;
186 
188  size_type M_size;
189 
191  std::vector<ublas::vector<double> > M_alpha;
192 
194  std::vector<ublas::vector<double> > M_beta;
195 
197  unknowns_type M_unknowns;
198 };
199 
200 template <typename SpaceType>
201 Bdf<SpaceType>::Bdf( space_ptrtype const& __space )
202  :
203  M_space( __space ),
204  M_size( 0 ),
205  M_alpha( BDF_MAX_ORDER ),
206  M_beta( BDF_MAX_ORDER )
207 {
208  for ( int i = 0; i < BDF_MAX_ORDER; ++i )
209  {
210  M_alpha[ i ].resize( i+2 );
211  M_beta[ i ].resize( i+1 );
212  }
213 
214  for ( int i = 0; i < BDF_MAX_ORDER; ++i )
215  {
216  if ( i == 0 ) // BDF_ORDER_ONE:
217  {
218  M_alpha[i][ 0 ] = 1.; // Backward Euler
219  M_alpha[i][ 1 ] = 1.;
220  M_beta[i][ 0 ] = 1.; // u^{n+1} \approx u^n
221  }
222 
223  else if ( i == 1 ) // BDF_ORDER_TWO:
224  {
225  M_alpha[i][ 0 ] = 3. / 2.;
226  M_alpha[i][ 1 ] = 2.;
227  M_alpha[i][ 2 ] = -1. / 2.;
228  M_beta[i][ 0 ] = 2.;
229  M_beta[i][ 1 ] = -1.;
230  }
231 
232  else if ( i == 2 ) // BDF_ORDER_THREE:
233  {
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.;
238  M_beta[i][ 0 ] = 3.;
239  M_beta[i][ 1 ] = -3.;
240  M_beta[i][ 2 ] = 1.;
241  }
242 
243  else if ( i == 3 )
244  {
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.;
250  M_beta[i][ 0 ] = 4.;
251  M_beta[i][ 1 ] = -6.;
252  M_beta[i][ 2 ] = 4.;
253  M_beta[i][ 3 ] = -1.;
254  }
255  }
256 
257  M_unknowns.resize( BDF_MAX_ORDER );
258 
259  for ( uint8_type __i = 0; __i < ( uint8_type )BDF_MAX_ORDER; ++__i )
260  {
261  M_unknowns[__i] = unknown_type( new element_type( M_space ) );
262  M_unknowns[__i]->zero();
263  }
264 }
265 
266 
267 template <typename SpaceType>
269 {}
270 
271 
272 template <typename SpaceType>
273 void
274 Bdf<SpaceType>::initialize( element_type const& u0 )
275 {
276  std::for_each( M_unknowns.begin(), M_unknowns.end(), *boost::lambda::_1 = u0 );
277 }
278 
279 template <typename SpaceType>
280 void
281 Bdf<SpaceType>::initialize( unknowns_type const& uv0 )
282 {
283  // Check if uv0 has the right dimensions
284  //FEELPP_ASSERT( uv0.size() == uint16_type(M_order) ).error( "Initial data set are not enough for the selected BDF" );
285 
286  std::copy( uv0.begin(), uv0.end(), M_unknowns.begin() );
287 }
288 
289 template <typename SpaceType>
290 const
291 typename Bdf<SpaceType>::unknowns_type&
293 {
294  return M_unknowns;
295 }
296 
297 template <typename SpaceType>
298 typename Bdf<SpaceType>::element_type&
300 {
301  //VLOG(1) << "[Bdf::unknown] id: " << i << " l2norm = " << M_unknowns[i]->l2Norm() << "\n";
302  return *M_unknowns[i];
303 }
304 
305 
306 
307 
308 template <typename SpaceType>
309 double
310 Bdf<SpaceType>::derivateCoefficient( int n, size_type i, double dt ) const
311 {
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;
314 }
315 
316 template <typename SpaceType>
317 double
318 Bdf<SpaceType>::extrapolateCoefficient( int n, size_type i, double dt ) const
319 {
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 ];
322 }
323 
324 template <typename SpaceType>
325 void
326 Bdf<SpaceType>::showMe( std::ostream& __out ) const
327 {
328 #if 0
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";
333 #endif
334 }
335 
336 template <typename SpaceType>
337 template<typename container_type>
338 void
339 Bdf<SpaceType>::shiftRight( typename space_type::template Element<value_type, container_type> const& __new_unk )
340 {
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 ) ) );
345  // u(t^{n}) coefficient is in M_unknowns[0]
346  *M_unknowns[0] = __new_unk;
347 
348  /* int i = 0;
349  BOOST_FOREACH( boost::shared_ptr<element_type>& t, M_unknowns )
350  {
351  //VLOG(1) << "[Bdf::shiftright] id: " << i << " l2norm = " << t->l2Norm() << "\n";
352  ++i;
353  }
354  */
355 }
356 
357 template <typename SpaceType>
358 typename Bdf<SpaceType>::element_type
359 Bdf<SpaceType>::derivate( int n, scalar_type dt ) const
360 {
361  element_type __t( M_space );
362  __t.zero();
363 
364  __t.add( ( 1./dt ), derivate( n ) );
365 
366  return __t;
367 }
368 
369 template <typename SpaceType>
370 typename Bdf<SpaceType>::element_type
372 {
373  element_type __t( M_space );
374  __t.zero();
375 
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" );
378 
379  for ( uint8_type i = 0; i < n; ++i )
380  __t.add( M_alpha[n-1][ i+1 ], *M_unknowns[i] );
381 
382  return __t;
383 }
384 
385 template <typename SpaceType>
386 typename Bdf<SpaceType>::element_type
387 Bdf<SpaceType>::extrapolate( int n ) const
388 {
389  element_type __t( M_space );
390  __t.zero();
391 
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" );
394 
395  for ( uint8_type i = 0; i < n; ++i )
396  __t.add( M_beta[n-1][ i ], *M_unknowns[ i ] );
397 
398  return __t;
399 }
400 
401 }
402 #endif

Generated on Fri Oct 25 2013 14:24:04 for Feel++ by doxygen 1.8.4