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
aitken2.hpp
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): Vincent Chabannes
6  Date: 24-03-2011
7 
8  Copyright (C) 2011 Universite Joseph Fourier (Grenoble I)
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 */
24 
25 
26 #ifndef __AitkenExtrapolation2
27 #define __AitkenExtrapolation2 1
28 
29 #include <feel/feeldiscr/functionspace.hpp>
30 #include <feel/feelalg/vector.hpp>
31 
32 
33 namespace Feel
34 {
64 template< typename fs_type >
65 class Aitken
66 {
67 
68 public:
69 
70  typedef fs_type functionspace_type;
71  typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
72 
73  typedef typename functionspace_type::element_type element_type;
74 
75  typedef typename functionspace_type::template Element<typename functionspace_type::value_type,
76  typename VectorUblas<typename functionspace_type::value_type>::range::type > element_range_type;
77 
78 
89  Aitken( functionspace_ptrtype _Xh, double _failsafeParameter = 0.1 )
90  :
91  Xh( _Xh ),
92  failsafeParameter( _failsafeParameter ),
93  previousParameter( _failsafeParameter ),
94  previousResidual( Xh, "previous residual" ),
95  previousElement( Xh, "previous element" ),
96  currentResidual( Xh, "current residual" ),
97  currentElement( Xh, "current element" )
98  {
99  }
100 
104  Aitken( Aitken const& tc )
105  :
106  Xh( tc.Xh ),
107  failsafeParameter( tc.failsafeParameter ),
108  previousParameter( tc.previousParameter ),
109  previousResidual( tc.previousResidual ),
110  previousElement( tc.previousElement ),
111  currentResidual( tc.currentResidual ),
112  currentElement( tc.currentElement )
113  {
114  }
115 
119  ~Aitken() {}
120 
126  void initialize( element_type const& residual, element_type const& elem )
127  {
128  previousResidual = residual;
129  previousElement = elem;
130  }
131 
132  void initialize( element_type const& residual, element_range_type const& elem )
133  {
134  previousResidual = residual;
135  previousElement.zero();
136  previousElement.add( 1.,elem );
137  /*previousElement = vf::project(previousElement.functionSpace(),
138  elements(previousElement.mesh()),
139  vf::idv(elem) );*/
140  }
141 
147  void setElement( element_type const& residual, element_type const& elem )
148  {
149  currentResidual = residual;
150  currentElement = elem;
151  }
152 
153  void setElement( element_type const& residual, element_range_type const& elem )
154  {
155  currentResidual = residual;
156  currentElement.zero();
157  currentElement.add( 1.,elem );
158  /*currentElement = vf::project(currentElement.functionSpace(),
159  elements(currentElement.mesh()),
160  vf::idv(elem) );*/
161  }
162 
163 
167  double calculateParameter();
168 
173  void relaxationStep( element_type& new_elem )
174  {
175  new_elem = currentResidual;
176  new_elem.scale( -previousParameter );
177 
178  new_elem += currentElement;
179  }
180 
185  void shiftRight()
186  {
187  previousResidual = currentResidual;
188  previousElement = currentElement;
189  }
190 
195  {
196  previousParameter = failsafeParameter;
197  }
198 
199 private:
200 
204  functionspace_ptrtype Xh;
205 
206  double failsafeParameter, previousParameter;
207 
208  element_type previousResidual, previousElement, currentResidual, currentElement;
209 
210 };
211 
212 
213 template< typename fs_type >
214 double
216 {
217  element_type aux( Xh, "aux" );
218 
219  aux = currentResidual;
220  aux -= previousResidual;
221 
222  double scalar = inner_product( aux, aux );
223 
224  aux.scale( 1.0/scalar );
225 
226  element_type aux2( Xh, "aux2" );
227 
228  aux2 = currentElement;
229  aux2 -= previousElement;
230 
231  scalar = inner_product( previousResidual , aux );
232  scalar = -previousParameter*scalar;
233 #if 1
234 
235  if ( scalar > 1 )
236  scalar = 1;//-failsafeParameter;
237 
238  if ( scalar < 0 )
239  scalar = failsafeParameter;
240 
241 #endif
242  previousParameter = scalar;
243 
244  return scalar;
245 }
246 
247 
248 } // End namespace Feel
249 
250 #endif

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