1/*!
2 * \file   UAnderson.ixx
3 * \brief
4 * \author Étienne Castelier
5 * \date   10 oct. 2016
6 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights
7 * reserved.
8 * This project is publicly released under either the GNU GPL Licence
9 * or the CECILL-A licence. A copy of thoses licences are delivered
10 * with the sources of TFEL. CEA or EDF may also distribute this
11 * project under specific licensing conditions.
12 */
13
14#ifndef LIB_TFEL_MATH_UANDERSON_IXX
15#define LIB_TFEL_MATH_UANDERSON_IXX
16
17#include <utility>
18
19namespace tfel {
20
21  namespace math {
22
23    //! constructor
24    template <typename Field, typename real>
25    UAnderson<Field, real>::UAnderson(
26        const typename AndersonBase<Field, real>::Allocator a)
27        : AndersonBase<Field, real>(std::move(a)) {
28    }  // end of UAnderson<Field,real>::UAnderson
29
30    // Displacement fields for a new iteration
31    // uO,uN: Old and new displacement field
32    template <typename Field, typename real>
33    void UAnderson<Field, real>::newIter(Field*& uO, Field*& uN) {
34      using size_type = typename AndersonBase<Field, real>::size_type;
35      // Matrice upgrade
36      if (this->n < this->size()) {
37        ++(this->n);
38      }
39      *uO -= *uN;
40      this->reset();
41      // Projection
42      if (this->n < this->size()) {
43        uO = this->D[this->n];
44      } else {
45        uO = this->D[0];
46      }
47      if ((1 < this->n) && (this->alt == this->alMax)) {
48        this->cM.weightsGSchmidtD(this->w);
49        anderson::linear_combinaison(*uO, this->u, this->w, this->n);
50        this->alt = 1;
51      } else {
52        *uO = *uN;
53        if (this->alt < this->alMax) {
54          ++(this->alt);
55        }
56      }
57      // Shift
58      if (this->n == this->size()) {
59        this->cM.shift();
60        uN = this->u[0];
61        for (size_type i = 0; i < this->n - 1; ++i) {
62          this->u[i] = this->u[i + 1];
63        }
64        this->u[this->n - 1] = uN;
65        for (size_type i = 0; i < this->n - 1; ++i) {
66          this->D[i] = this->D[i + 1];
67        }
68        this->D[this->n - 1] = uO;
69      } else {
70        uN = this->u[this->n];
71      }
72    }
73
74    // First Iteration of a new time step
75    // uO,uN: Old and new displacement field
76    template <typename Field, typename real>
77    void UAnderson<Field, real>::restart(Field*& uO, Field*& uN) {
78      this->alloc();
79      uO = this->D[0];
80      uN = this->u[0];
81      this->alt = this->alMax;
82    }
83
84  }  // end of namespace math
85
86}  // end of namespace tfel
87
88#endif /* LIB_TFEL_MATH_UANDERSON_IXX */
89