GetFEM  5.4.2
getfem_continuation.cc
Go to the documentation of this file.
1 /*===========================================================================
2 
3  Copyright (C) 2011-2020 Tomas Ligursky, Yves Renard, Konstantinos Poulios
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 /** @file getfem_continuation.cc
23  @author Tomas Ligursky <tomas.ligursky@ugn.cas.cz>
24  @author Yves Renard <Yves.Renard@insa-lyon.fr>
25  @author Konstantinos Poulios <logari81@googlemail.com>
26  @date January 12, 2014.
27  @brief inexact Moore-Penrose continuation method.
28 */
29 
31 
32 namespace getfem {
33 
34  void cont_struct_getfem_model::set_variables
35  (const base_vector &x, double gamma) const {
36  md->set_real_variable(parameter_name)[0] = gamma;
37  if (!currentdata_name.empty()) {
38  gmm::add(gmm::scaled(md->real_variable(initdata_name), 1. - gamma),
39  gmm::scaled(md->real_variable(finaldata_name), gamma),
40  md->set_real_variable(currentdata_name));
41  }
42  md->to_variables(x);
43  }
44 
45  void cont_struct_getfem_model::update_matrix
46  (const base_vector &x, double gamma) const {
47  set_variables(x, gamma);
48  if (noisy() > 2) cout << "starting computing tangent matrix" << endl;
49  md->assembly(model::BUILD_MATRIX);
50  }
51 
52  // solve A * g = L
53  void cont_struct_getfem_model::solve
54  (const model_real_sparse_matrix &A, base_vector &g,
55  const base_vector &L) const {
56  if (noisy() > 2) cout << "starting linear solver" << endl;
57  gmm::iteration iter(maxres_solve, (noisy() >= 2) ? noisy() - 2 : 0,
58  40000);
59  (*lsolver)(A, g, L, iter);
60  if (noisy() > 2) cout << "linear solver done" << endl;
61  }
62 
63  // solve A * (g1|g2) = (L1|L2)
64  void cont_struct_getfem_model::solve
65  (const model_real_sparse_matrix &A, base_vector &g1, base_vector &g2,
66  const base_vector &L1, const base_vector &L2) const {
67  if (noisy() > 2) cout << "starting linear solver" << endl;
68  gmm::iteration iter(maxres_solve, (noisy() >= 2) ? noisy() - 2 : 0,
69  40000);
70  (*lsolver)(A, g1, L1, iter);
71  iter.init(); (*lsolver)(A, g2, L2, iter); // (can be optimised)
72  if (noisy() > 2) cout << "linear solver done" << endl;
73  }
74 
75  // F(x, gamma) --> f
76  void cont_struct_getfem_model::F
77  (const base_vector &x, double gamma, base_vector &f) const {
78  set_variables(x, gamma);
79  md->assembly(model::BUILD_RHS);
80  gmm::copy(gmm::scaled(md->real_rhs(), -1.), f);
81  }
82 
83  // (F(x, gamma + eps) - f0) / eps --> g
84  void cont_struct_getfem_model::F_gamma
85  (const base_vector &x, double gamma, const base_vector &f0,
86  base_vector &g) const {
87  const double eps = diffeps;
88  F(x, gamma + eps, g);
89  gmm::add(gmm::scaled(f0, -1.), g);
90  gmm::scale(g, 1./eps);
91  }
92 
93  // (F(x, gamma + eps) - F(x, gamma)) / eps --> g
94  void cont_struct_getfem_model::F_gamma
95  (const base_vector &x, double gamma, base_vector &g) const {
96  base_vector f0(x);
97  F(x, gamma, f0);
98  F_gamma(x, gamma, f0, g);
99  }
100 
101  // F_x(x, gamma) --> A
102  void cont_struct_getfem_model::F_x
103  (const base_vector &x, double gamma, model_real_sparse_matrix &A) const {
104  update_matrix(x, gamma);
105  size_type nbdof = md->nb_dof();
106  gmm::resize(A, nbdof, nbdof);
107  gmm::copy(md->real_tangent_matrix(), A);
108  }
109 
110  // solve F_x(x, gamma) * g = L
111  void cont_struct_getfem_model::solve_grad
112  (const base_vector &x, double gamma, base_vector &g,
113  const base_vector &L) const {
114  update_matrix(x, gamma);
115  solve(md->real_tangent_matrix(), g, L);
116 //int exp;
117 //cout<<"det="<<MUMPS_determinant(md->real_tangent_matrix(),exp)<<"*2^"<<exp<<endl;
118  }
119 
120  // solve F_x(x, gamma) * (g1|g2) = (L1|L2)
121  void cont_struct_getfem_model::solve_grad
122  (const base_vector &x, double gamma, base_vector &g1, base_vector &g2,
123  const base_vector &L1, const base_vector &L2) const {
124  update_matrix(x, gamma);
125  solve(md->real_tangent_matrix(), g1, g2, L1, L2);
126  }
127 
128  // F_x(x, gamma) * w --> y
129  void cont_struct_getfem_model::mult_grad
130  (const base_vector &x, double gamma,
131  const base_vector &w, base_vector &y) const {
132  update_matrix(x, gamma);
133  mult(md->real_tangent_matrix(), w, y);
134  }
135 
136  size_type cont_struct_getfem_model::estimated_memsize(void) {
137  return sizeof(cont_struct_getfem_model)
138  + virtual_cont_struct::estimated_memsize();
139  }
140 
141 } /* end of namespace getfem. */
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem_continuation.h
Inexact Moore-Penrose continuation method.