GetFEM  5.4.2
getfem_continuation.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2011-2020 Tomas Ligursky, Yves Renard, Konstantinos Poulios
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /** @file getfem_continuation.h
33  @author Tomas Ligursky <tomas.ligursky@ugn.cas.cz>
34  @author Yves Renard <Yves.Renard@insa-lyon.fr>
35  @author Konstantinos Poulios <logari81@googlemail.com>
36  @date October 17, 2011.
37  @brief Inexact Moore-Penrose continuation method.
38 */
39 #ifndef GETFEM_CONTINUATION_H__
40 #define GETFEM_CONTINUATION_H__
41 
42 #include <getfem/getfem_model_solvers.h>
43 
44 namespace getfem {
45 
46 
47  //=========================================================================
48  // Abstract Moore-Penrose continuation method
49  //=========================================================================
50 
51  template <typename VECT, typename MAT>
52  class virtual_cont_struct {
53 
54  protected:
55 #ifdef _MSC_VER
56  const double tau_bp_init = 1.e6;
57  const double diffeps = 1.e-8;
58 #else
59  static constexpr double tau_bp_init = 1.e6;
60  static constexpr double diffeps = 1.e-8;
61 #endif
62 
63  int singularities;
64 
65  private:
66 
67  bool non_smooth;
68  double scfac, h_init_, h_max_, h_min_, h_inc_, h_dec_;
69  size_type maxit_, thrit_;
70  double maxres_, maxdiff_, mincos_, delta_max_, delta_min_, thrvar_;
71  size_type nbdir_, nbspan_;
72  int noisy_;
73  double tau_lp, tau_bp_1, tau_bp_2;
74 
75  // stored singularities info
76  std::map<double, double> tau_bp_graph;
77  VECT alpha_hist, tau_bp_hist;
78  std::string sing_label;
79  VECT x_sing, x_next;
80  double gamma_sing, gamma_next;
81  std::vector<VECT> tx_sing, tx_predict;
82  std::vector<double> tgamma_sing, tgamma_predict;
83 
84  // randomized data
85  VECT bb_x_, cc_x_;
86  double bb_gamma, cc_gamma, dd;
87 
88  public:
89  /* Compute a unit tangent at (x, gamma) that is accute to the incoming
90  vector. */
91  void compute_tangent(const VECT &x, double gamma,
92  VECT &tx, double &tgamma) {
93  VECT g(x), y(x);
94  F_gamma(x, gamma, g); // g = F_gamma(x, gamma)
95  solve_grad(x, gamma, y, g); // y = F_x(x, gamma)^-1 * g
96  tgamma = 1. / (tgamma - w_sp(tx, y));
97  gmm::copy(gmm::scaled(y, -tgamma), tx); // tx = -tgamma * y
98 
99  scale(tx, tgamma, 1./w_norm(tx, tgamma)); // [tx,tgamma] /= w_norm(tx,tgamma)
100 
101  mult_grad(x, gamma, tx, y); // y = F_x(x, gamma) * tx
102  gmm::add(gmm::scaled(g, tgamma), y); // y += tgamma * g
103  double r = norm(y);
104  if (r > 1.e-10)
105  GMM_WARNING2("Tangent computed with the residual " << r);
106  }
107 
108  private:
109 
110  /* Calculate a tangent vector at (x, gamma) + h * (tX, tGamma) and test
111  whether it is close to (tX, tGamma). Informatively, compare it with
112  (tx, tgamma), as well. */
113  bool test_tangent(const VECT &x, double gamma,
114  const VECT &tX, double tGamma,
115  const VECT &tx, double tgamma, double h) {
116  bool res = false;
117  double Gamma1, tGamma1(tgamma);
118  VECT X1(x), tX1(tx);
119 
120  scaled_add(x, gamma, tX, tGamma, h, X1, Gamma1); // [X1,Gamma1] = [x,gamma] + h * [tX,tGamma]
121  compute_tangent(X1, Gamma1, tX1, tGamma1);
122 
123  double cang = cosang(tX1, tX, tGamma1, tGamma);
124  if (noisy() > 1)
125  cout << "cos of the angle with the tested tangent " << cang << endl;
126  if (cang >= mincos())
127  res = true;
128  else {
129  cang = cosang(tX1, tx, tGamma1, tGamma);
130  if (noisy() > 1)
131  cout << "cos of the angle with the initial tangent " << cang << endl;
132  }
133  return res;
134  }
135 
136  /* Simple tangent switch. */
137  bool switch_tangent(const VECT &x, double gamma,
138  VECT &tx, double &tgamma, double &h) {
139  double Gamma, tGamma(tgamma);
140  VECT X(x), tX(tx);
141 
142  if (noisy() > 0) cout << "Trying a simple tangent switch" << endl;
143  if (noisy() > 1) cout << "Computing a new tangent" << endl;
144  h *= 1.5;
145  scaled_add(x, gamma, tx, tgamma, h, X, Gamma); // [X,Gamma] = [x,gamma] + h * [tx,tgamma]
146  compute_tangent(X, Gamma, tX, tGamma);
147  // One can test the cosine of the angle between (tX, tGamma) and
148  // (tx, tgamma), for sure, and increase h_min if it were greater or
149  // equal to mincos(). However, this seems to be superfluous.
150 
151  if (noisy() > 1)
152  cout << "Starting testing the computed tangent" << endl;
153  double h_test = -0.9 * h_min();
154  bool accepted(false);
155  while (!accepted && (h_test > -h_max())) {
156  h_test = -h_test
157  + pow(10., floor(log10(-h_test / h_min()))) * h_min();
158  accepted = test_tangent(x, gamma, tX, tGamma, tx, tgamma, h_test);
159  if (!accepted) {
160  h_test *= -1.;
161  accepted = test_tangent(x, gamma, tX, tGamma, tx, tgamma, h_test);
162  }
163  }
164 
165  if (accepted) {
166  if (h_test < 0) {
167  gmm::scale(tX, -1.);
168  tGamma *= -1.;
169  h_test *= -1.;
170  }
171  if (noisy() > 0)
172  cout << "Tangent direction switched, "
173  << "starting computing a suitable step size" << endl;
174  h = h_init();
175  bool h_adapted = false;
176  while (!h_adapted && (h > h_test)) {
177  h_adapted = test_tangent(x, gamma, tX, tGamma, tx, tgamma, h);
178  h *= h_dec();
179  }
180  h = h_adapted ? h / h_dec() : h_test;
181  copy(tX, tGamma, tx, tgamma);
182  } else
183  if (noisy() > 0) cout << "Simple tangent switch has failed!" << endl;
184 
185  return accepted;
186  }
187 
188  /* Test for limit points (also called folds or turning points). */
189  bool test_limit_point(double tgamma) {
190  double tau_lp_old = get_tau_lp();
191  set_tau_lp(tgamma);
192  return (tgamma * tau_lp_old < 0);
193  }
194 
195  void init_test_functions(const VECT &x, double gamma,
196  const VECT &tx, double tgamma) {
197  set_tau_lp(tgamma);
198  if (this->singularities > 1) {
199  if (noisy() > 1) cout << "Computing an initial value of the "
200  << "test function for bifurcations" << endl;
201  set_tau_bp_2(test_function_bp(x, gamma, tx, tgamma));
202  }
203  }
204 
205  /* Test function for bifurcation points for a given matrix. The first part
206  of the solution of the augmented system is passed in
207  (v_x, v_gamma). */
208  double test_function_bp(const MAT &A, const VECT &g,
209  const VECT &tx, double tgamma,
210  VECT &v_x, double &v_gamma) {
211  VECT y(g), z(g);
212  size_type nn = gmm::vect_size(g);
213 
214  solve(A, y, z, g, bb_x(nn)); // [y,z] = A^-1 * [g,bb_x]
215  v_gamma = (bb_gamma - sp(tx, z)) / (tgamma - sp(tx, y));
216  scaled_add(z, y, -v_gamma, v_x); // v_x = y - v_gamma*z
217  double tau = 1. / (dd - sp(cc_x(nn), v_x) - cc_gamma * v_gamma);
218  scale(v_x, v_gamma, -tau); // [v_x,v_gamma] *= -tau
219 
220  // control of the norm of the residual
221  mult(A, v_x, y);
222  gmm::add(gmm::scaled(g, v_gamma), y); // y += v_gamma*g
223  gmm::add(gmm::scaled(bb_x(nn), tau), y); // y += bb_x*tau
224  double r = sp(tx, v_x) + tgamma * v_gamma + bb_gamma * tau;
225  double q = sp(cc_x(nn), v_x) + cc_gamma * v_gamma + dd * tau - 1.;
226  r = sqrt(sp(y, y) + r * r + q * q);
227  if (r > 1.e-10)
228  GMM_WARNING2("Test function evaluated with the residual " << r);
229 
230  return tau;
231  }
232 
233  double test_function_bp(const MAT &A, const VECT &g,
234  const VECT &tx, double tgamma) {
235  VECT v_x(g); double v_gamma;
236  return test_function_bp(A, g, tx, tgamma, v_x, v_gamma);
237  }
238 
239  /* Test function for bifurcation points for the gradient computed at
240  (x, gamma). */
241  double test_function_bp(const VECT &x, double gamma,
242  const VECT &tx, double tgamma,
243  VECT &v_x, double &v_gamma) {
244  MAT A; VECT g(x);
245  F_x(x, gamma, A);
246  F_gamma(x, gamma, g);
247  return test_function_bp(A, g, tx, tgamma, v_x, v_gamma);
248  }
249 
250  double test_function_bp(const VECT &x, double gamma,
251  const VECT &tx, double tgamma) {
252  VECT v_x(x); double v_gamma;
253  return test_function_bp(x, gamma, tx, tgamma, v_x, v_gamma);
254  }
255 
256  public:
257  /* Test for smooth bifurcation points. */
258  bool test_smooth_bifurcation(const VECT &x, double gamma,
259  const VECT &tx, double tgamma) {
260  double tau_bp_1_old = get_tau_bp_1();
261  set_tau_bp_1(get_tau_bp_2());
262  set_tau_bp_2(test_function_bp(x, gamma, tx, tgamma));
263  return (get_tau_bp_2() * get_tau_bp_1() < 0) &&
264  (gmm::abs(get_tau_bp_1()) < gmm::abs(tau_bp_1_old));
265  }
266 
267  /* Test for non-smooth bifurcation points. */
268  bool test_nonsmooth_bifurcation(const VECT &x1, double gamma1,
269  const VECT &tx1, double tgamma1,
270  const VECT &x2, double gamma2,
271  const VECT &tx2, double tgamma2) {
272  VECT g1(x1), g2(x1), g(x1), tx(x1);
273 
274  // compute gradients at the two given points
275  MAT A1, A2;
276  F_x(x2, gamma2, A2);
277  F_gamma(x2, gamma2, g2);
278  F_x(x1, gamma1, A1);
279  F_gamma(x1, gamma1, g1);
280  double tau1 = test_function_bp(A1, g1, tx1, tgamma1);
281  double tau2 = test_function_bp(A2, g2, tx2, tgamma2);
282  double tau_var_ref = std::max(gmm::abs(tau2 - tau1), 1.e-8);
283  set_tau_bp_2(tau1);
284  init_tau_bp_graph();
285  MAT A(A2);
286 
287  // monitor the sign changes of the test function on the convex
288  // combination
289  size_type nb_changes = 0;
290  double delta = delta_min(), tau0 = tau_bp_init, tgamma;
291  for (double alpha=0.; alpha < 1.; ) {
292  alpha = std::min(alpha + delta, 1.);
293  scaled_add(A1, 1.-alpha, A2, alpha, A); // A = (1-alpha)*A1 + alpha*A2
294  scaled_add(g1, 1.-alpha, g2, alpha, g); // g = (1-alpha)*g1 + alpha*g2
295  scaled_add(tx1, tgamma1, 1.-alpha, tx2, tgamma2, alpha, tx, tgamma);
296  //[tx,tgamma] = (1-alpha)*[tx1,tgamma1] + alpha*[tx2,tgamma2]
297 
298  tau2 = test_function_bp(A, g, tx, tgamma);
299  if ((tau2 * tau1 < 0) && (gmm::abs(tau1) < gmm::abs(tau0)))
300  ++nb_changes;
301  insert_tau_bp_graph(alpha, tau2);
302 
303  if (gmm::abs(tau2 - tau1) < 0.5 * thrvar() * tau_var_ref)
304  delta = std::min(2 * delta, delta_max());
305  else if (gmm::abs(tau2 - tau1) > thrvar() * tau_var_ref)
306  delta = std::max(0.1 * delta, delta_min());
307  tau0 = tau1;
308  tau1 = tau2;
309  }
310 
311  set_tau_bp_1(tau_bp_init);
312  set_tau_bp_2(tau2);
313  return nb_changes % 2;
314  }
315 
316  private:
317  /* Newton-type corrections for the couple ((X, Gamma), (tX, tGamma)).
318  The current direction of (tX, tGamma) is informatively compared with
319  (tx, tgamma). */
320  bool newton_corr(VECT &X, double &Gamma, VECT &tX,
321  double &tGamma, const VECT &tx, double tgamma,
322  size_type &it) {
323  bool converged = false;
324  double Delta_Gamma, res(0), diff;
325  VECT f(X), g(X), Delta_X(X), y(X);
326 
327  if (noisy() > 1) cout << "Starting correction" << endl;
328  F(X, Gamma, f); // f = F(X, Gamma) = -rhs(X, Gamma)
329 //CHANGE 1: line search
330 //double res0 = norm(f);
331 
332  for (it=0; it < maxit() && res < 1.e8; ++it) {
333  F_gamma(X, Gamma, f, g); // g = F_gamma(X, Gamma)
334  solve_grad(X, Gamma, Delta_X, y, f, g); // y = F_x(X, Gamma)^-1 * g
335  // Delta_X = F_x(X, Gamma)^-1 * f
336  Delta_Gamma = sp(tX, Delta_X) / (sp(tX, y) - tGamma); // Delta_Gamma = tX.Delta_X / (tX.y - tGamma)
337  if (isnan(Delta_Gamma)) {
338  if (noisy() > 1) cout << "Newton correction failed with NaN" << endl;
339  return false;
340  }
341  gmm::add(gmm::scaled(y, -Delta_Gamma), Delta_X); // Delta_X -= Delta_Gamma * y
342  scaled_add(X, Gamma, Delta_X, Delta_Gamma, -1,
343  X, Gamma); // [X,Gamma] -= [Delta_X,Delta_Gamma]
344  F(X, Gamma, f); // f = F(X, gamma) = -rhs(X, gamma)
345  res = norm(f);
346 
347 //CHANGE 1: line search
348 //for (size_type ii=0; ii < 4 && (isnan(res) || res > res0); ++ii) { // some basic linesearch
349 // scale(Delta_X, Delta_Gamma, 0.5);
350 // scaled_add(X, Gamma, Delta_X, Delta_Gamma, 1, X, Gamma); // [X,Gamma] += [Delta_X,Delta_Gamma]
351 // F(X, Gamma, f); // f = F(X, gamma) = -rhs(X, gamma)
352 // res = norm(f);
353 //}
354 
355  tGamma = 1. / (tGamma - w_sp(tX, y)); // tGamma = 1 / (tGamma - k*tX.y)
356  gmm::copy(gmm::scaled(y, -tGamma), tX); // tX = -tGamma * y
357  scale(tX, tGamma, 1./w_norm(tX, tGamma)); // [tX,tGamma] /= w_norm(tX,tGamma)
358 
359  diff = w_norm(Delta_X, Delta_Gamma);
360  if (noisy() > 1)
361  cout << " Correction " << std::setw(3) << it << ":"
362  << " Gamma = " << std::fixed << std::setprecision(6) << Gamma
363  << " residual = " << std::scientific << std::setprecision(3) << res
364  << " difference = " << std::scientific << std::setprecision(3) << diff
365  << " cosang = " << std::fixed << std::setprecision(6)
366  << cosang(tX, tx, tGamma, tgamma) << endl;
367 
368  if (res <= maxres() && diff <= maxdiff()) {
369  converged = true;
370  // recalculate the final tangent, for sure
371  compute_tangent(X, Gamma, tX, tGamma);
372  break;
373  }
374  }
375  if (noisy() > 1) cout << "Correction finished with Gamma = "
376  << Gamma << endl;
377  return converged;
378  }
379 
380  bool newton_corr(VECT &X, double &Gamma, VECT &tX,
381  double &tGamma, const VECT &tx, double tgamma) {
382  size_type it;
383  return newton_corr(X, Gamma, tX, tGamma, tx, tgamma, it);
384  }
385 
386  /* Try to perform one predictor-corrector step starting from the couple
387  ((x, gamma), (tx, tgamma)). Return the resulting couple in the case of
388  convergence. */
389  bool test_predict_dir(VECT &x, double &gamma,
390  VECT &tx, double &tgamma) {
391  bool converged = false;
392  double h = h_init(), Gamma, tGamma;
393  VECT X(x), tX(x);
394  while (!converged) { //step control
395  // prediction
396  scaled_add(x, gamma, tx, tgamma, h, X, Gamma); // [X,Gamma] = [x,gamma] + h * [tx,tgamma]
397  if (noisy() > 1)
398  cout << "(TPD) Prediction : Gamma = " << Gamma
399  << " (for h = " << h << ", tgamma = " << tgamma << ")" << endl;
400  copy(tx, tgamma, tX, tGamma);
401  //correction
402  converged = newton_corr(X, Gamma, tX, tGamma, tx, tgamma);
403 
404  if (h > h_min())
405  h = std::max(0.199 * h_dec() * h, h_min());
406  else
407  break;
408  }
409  if (converged) {
410  // check the direction of the tangent found
411  scaled_add(X, Gamma, x, gamma, -1., tx, tgamma); // [tx,tgamma] = [X,Gamma] - [x,gamma]
412  if (sp(tX, tx, tGamma, tgamma) < 0)
413  scale(tX, tGamma, -1.); // [tX,tGamma] *= -1
414  copy(X, Gamma, x, gamma);
415  copy(tX, tGamma, tx, tgamma);
416  }
417  return converged;
418  }
419 
420  /* A tool for approximating a smooth bifurcation point close to (x, gamma)
421  and locating the two branches emanating from there. */
422  void treat_smooth_bif_point(const VECT &x, double gamma,
423  const VECT &tx, double tgamma, double h) {
424  double tau0(get_tau_bp_1()), tau1(get_tau_bp_2());
425  double gamma0(gamma), Gamma,
426  tgamma0(tgamma), tGamma(tgamma), v_gamma;
427  VECT x0(x), X(x), tx0(tx), tX(tx), v_x(tx);
428 
429  if (noisy() > 0)
430  cout << "Starting locating a bifurcation point" << endl;
431 
432  // predictor-corrector steps with a secant-type step-length adaptation
433  h *= tau1 / (tau0 - tau1);
434  for (size_type i=0; i < 10 && (gmm::abs(h) >= h_min()); ++i) {
435  scaled_add(x0, gamma0, tx0, tgamma0, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h * [tx0,tgamma0]
436  if (noisy() > 1)
437  cout << "(TSBP) Prediction : Gamma = " << Gamma
438  << " (for h = " << h << ", tgamma = " << tgamma << ")" << endl;
439  if (newton_corr(X, Gamma, tX, tGamma, tx0, tgamma0)) {
440  copy(X, Gamma, x0, gamma0);
441  if (cosang(tX, tx0, tGamma, tgamma0) >= mincos())
442  copy(tX, tGamma, tx0, tgamma0);
443  tau0 = tau1;
444  tau1 = test_function_bp(X, Gamma, tx0, tgamma0, v_x, v_gamma);
445  h *= tau1 / (tau0 - tau1);
446  } else {
447  scaled_add(x0, gamma0, tx0, tgamma0, h, x0, gamma0); // [x0,gamma0] += h*[tx0,tgamma0]
448  test_function_bp(x0, gamma0, tx0, tgamma0, v_x, v_gamma);
449  break;
450  }
451  }
452  if (noisy() > 0)
453  cout << "Bifurcation point located" << endl;
454  set_sing_point(x0, gamma0);
455  insert_tangent_sing(tx0, tgamma0);
456 
457  if (noisy() > 0)
458  cout << "Starting searching for the second branch" << endl;
459  double no = w_norm(v_x, v_gamma);
460  scale(v_x, v_gamma, 1./no); // [v_x,v_gamma] /= no
461  if (test_predict_dir(x0, gamma0, v_x, v_gamma)
462  && insert_tangent_sing(v_x, v_gamma))
463  { if (noisy() > 0) cout << "Second branch found" << endl; }
464  else if (noisy() > 0) cout << "Second branch not found!" << endl;
465  }
466 
467  public:
468 
469  /* A tool for approximating a non-smooth point close to (x, gamma) and
470  consequent locating one-sided smooth solution branches emanating from
471  there. It is supposed that (x, gamma) is a point on a smooth solution
472  branch within the distance of h_min() from the end point of this
473  branch and (tx, tgamma) is the corresponding tangent directed towards
474  the end point. The boolean set_next determines whether the first new
475  branch found (if any) is to be chosen for further continuation. */
476  void treat_nonsmooth_point(const VECT &x, double gamma,
477  const VECT &tx, double tgamma, bool set_next) {
478  double gamma0(gamma), Gamma(gamma);
479  double tgamma0(tgamma), tGamma(tgamma);
480  double h = h_min(), mcos = mincos();
481  VECT x0(x), X(x), tx0(tx), tX(tx);
482 
483  // approximate the end point more precisely by a bisection-like algorithm
484  if (noisy() > 0)
485  cout << "Starting locating a non-smooth point" << endl;
486 
487  scaled_add(x0, gamma0, tx0, tgamma0, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h*[tx0,tgamma0]
488  if (newton_corr(X, Gamma, tX, tGamma, tx0, tgamma0)) { // --> X, Gamma, tX, tGamma
489  double cang = cosang(tX, tx0, tGamma, tgamma0);
490  if (cang >= mcos) mcos = (cang + 1.) / 2.;
491  }
492 
493  copy(tx0, tgamma0, tX, tGamma);
494  h /= 2.;
495  for (size_type i = 0; i < 15; i++) {
496  scaled_add(x0, gamma0, tx0, tgamma0, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h*[tx0,tgamma0]
497  if (noisy() > 1)
498  cout << "(TNSBP) Prediction : Gamma = " << Gamma
499  << " (for h = " << h << ", tgamma = " << tgamma << ")" << endl;
500  if (newton_corr(X, Gamma, tX, tGamma, tx0, tgamma0)
501  && (cosang(tX, tx0, tGamma, tgamma0) >= mcos)) {
502  copy(X, Gamma, x0, gamma0);
503  copy(tX, tGamma, tx0, tgamma0);
504  } else {
505  copy(tx0, tgamma0, tX, tGamma);
506  }
507  h /= 2.;
508  }
509  if (noisy() > 0)
510  cout << "Non-smooth point located" << endl;
511  set_sing_point(x0, gamma0);
512 
513  // take two reference vectors to span a subspace of directions emanating
514  // from the end point
515  if (noisy() > 0)
516  cout << "Starting a thorough search for different branches" << endl;
517  double tgamma1 = tgamma0, tgamma2 = tgamma0;
518  VECT tx1(tx0), tx2(tx0);
519  scale(tx1, tgamma1, -1.); // [tx1,tgamma1] *= -1
520  insert_tangent_sing(tx1, tgamma1);
521  h = h_min();
522  scaled_add(x0, gamma0, tx0, tgamma0, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h*[tx0,tgamma0]
523  compute_tangent(X, Gamma, tx2, tgamma2);
524 
525  // try systematically the directions of linear combinations of the couple
526  // of the reference vectors for finding new possible tangent predictions
527  // emanating from the end point
528  size_type i1 = 0, i2 = 0, nspan = 0;
529  double a, a1, a2, no;
530 
531  do {
532  for (size_type i = 0; i < nbdir(); i++) {
533  a = (2 * M_PI * double(i)) / double(nbdir());
534  a1 = sin(a);
535  a2 = cos(a);
536  scaled_add(tx1, tgamma1, a1, tx2, tgamma2, a2, tX, tGamma); // [tX,tGamma] = a1*[tx1,tgamma1] + a2*[tx2,tgamma2]
537  no = w_norm(tX, tGamma);
538  scaled_add(x0, gamma0, tX, tGamma, h/no, X, Gamma); // [X,Gamma] = [x0,gamma0] + h/no * [tX,tGamma]
539  compute_tangent(X, Gamma, tX, tGamma);
540 
541  if (gmm::abs(cosang(tX, tx0, tGamma, tgamma0)) < mincos()
542  || (i == 0 && nspan == 0)) {
543  copy(tX, tGamma, tx0, tgamma0);
544  if (insert_tangent_predict(tX, tGamma)) {
545  if (noisy() > 1)
546  cout << "New potential tangent vector found, "
547  << "trying one predictor-corrector step" << endl;
548  copy(x0, gamma0, X, Gamma);
549 
550  if (test_predict_dir(X, Gamma, tX, tGamma)) {
551  if (insert_tangent_sing(tX, tGamma)) {
552  if ((i == 0) && (nspan == 0)
553  // => (tX, tGamma) = (tx2, tgamma2)
554  && (gmm::abs(cosang(tX, tx0, tGamma, tgamma0))
555  >= mincos())) { i2 = 1; }
556  if (noisy() > 0) cout << "A new branch located (for nspan = "
557  << nspan << ")" << endl;
558  if (set_next) set_next_point(X, Gamma);
559 
560  }
561  copy(x0, gamma0, X, Gamma);
562  copy(tx0, tgamma0, tX, tGamma);
563  }
564 
565  scale(tX, tGamma, -1.); // [tX,tGamma] *= -1
566  if (test_predict_dir(X, Gamma, tX, tGamma)
567  && insert_tangent_sing(tX, tGamma)) {
568  if (noisy() > 0) cout << "A new branch located (for nspan = "
569  << nspan << ")" << endl;
570  if (set_next) set_next_point(X, Gamma);
571  }
572  }
573  }
574  }
575 
576  // heuristics for varying the reference vectors
577  bool perturb = true;
578  if (i1 + 1 < i2) { ++i1; perturb = false; }
579  else if(i2 + 1 < nb_tangent_sing())
580  { ++i2; i1 = 0; perturb = false; }
581  if (!perturb) {
582  copy(get_tx_sing(i1), get_tgamma_sing(i1), tx1, tgamma1);
583  copy(get_tx_sing(i2), get_tgamma_sing(i2), tx2, tgamma2);
584  } else {
585  gmm::fill_random(tX);
586  tGamma = gmm::random(1.);
587  no = w_norm(tX, tGamma);
588  scaled_add(tx2, tgamma2, tX, tGamma, 0.1/no, tx2, tgamma2);
589  // [tx2,tgamma2] += 0.1/no * [tX,tGamma]
590  scaled_add(x0, gamma0, tx2, tgamma2, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h*[tx2,tgamma2]
591  compute_tangent(X, Gamma, tx2, tgamma2);
592  }
593  } while (++nspan < nbspan());
594 
595  if (noisy() > 0)
596  cout << "Number of branches emanating from the non-smooth point "
597  << nb_tangent_sing() << endl;
598  }
599 
600 
601  void init_Moore_Penrose_continuation(const VECT &x,
602  double gamma, VECT &tx,
603  double &tgamma, double &h) {
604  gmm::clear(tx);
605  tgamma = (tgamma >= 0) ? 1. : -1.;
606  if (noisy() > 1)
607  cout << "Computing an initial tangent" << endl;
608  compute_tangent(x, gamma, tx, tgamma);
609  h = h_init();
610  if (this->singularities > 0)
611  init_test_functions(x, gamma, tx, tgamma);
612  }
613 
614 
615  /* Perform one step of the (non-smooth) Moore-Penrose continuation.
616  NOTE: The new point need not to be saved in the model in the end! */
617  void Moore_Penrose_continuation(VECT &x, double &gamma,
618  VECT &tx, double &tgamma,
619  double &h, double &h0) {
620  bool converged, new_point = false, tangent_switched = false;
621  size_type it, step_dec = 0;
622  double tgamma0 = tgamma, Gamma, tGamma;
623  VECT tx0(tx), X(x), tX(x);
624 
625  clear_tau_bp_currentstep();
626  clear_sing_data();
627 
628  do {
629  h0 = h;
630  // prediction
631  scaled_add(x, gamma, tx, tgamma, h, X, Gamma); // [X,Gamma] = [x,gamma] + h*[tx,tgamma]
632  if (noisy() > 1)
633  cout << " Prediction : Gamma = " << Gamma
634  << " (for h = " << std::scientific << std::setprecision(3) << h
635  << ", tgamma = " << tgamma << ")" << endl;
636  copy(tx, tgamma, tX, tGamma);
637 
638  // correction
639  converged = newton_corr(X, Gamma, tX, tGamma, tx, tgamma, it);
640  double cang(converged ? cosang(tX, tx, tGamma, tgamma) : 0);
641 
642  if (converged && cang >= mincos()) {
643  new_point = true;
644  if (this->singularities > 0) {
645  if (test_limit_point(tGamma)) {
646  set_sing_label("limit point");
647  if (noisy() > 0) cout << "Limit point detected!" << endl;
648  }
649  if (this->singularities > 1) { // Treat bifurcations
650  if (noisy() > 1)
651  cout << "New point found, computing a test function "
652  << "for bifurcations" << endl;
653  if (!tangent_switched) {
654  if (test_smooth_bifurcation(X, Gamma, tX, tGamma)) {
655  set_sing_label("smooth bifurcation point");
656  if (noisy() > 0)
657  cout << "Smooth bifurcation point detected!" << endl;
658  treat_smooth_bif_point(X, Gamma, tX, tGamma, h);
659  }
660  } else if (test_nonsmooth_bifurcation(x, gamma, tx0,
661  tgamma0, X, Gamma, tX,
662  tGamma)) {
663  set_sing_label("non-smooth bifurcation point");
664  if (noisy() > 0)
665  cout << "Non-smooth bifurcation point detected!" << endl;
666  treat_nonsmooth_point(x, gamma, tx0, tgamma0, false);
667  }
668  }
669  }
670 
671 //CHANGE 2: avoid false step increases
672 //if (step_dec == 0 && it < thrit() && h_inc()*(1-cang) < (1-mincos()))
673  if (step_dec == 0 && it < thrit())
674  h = std::min(h_inc() * h, h_max());
675  } else if (h > h_min()) {
676  h = std::max(h_dec() * h, h_min());
677  step_dec++;
678  } else if (this->non_smooth && !tangent_switched) {
679  if (noisy() > 0)
680  cout << "Classical continuation has failed!" << endl;
681  if (switch_tangent(x, gamma, tx, tgamma, h)) {
682  tangent_switched = true;
683  step_dec = (h >= h_init()) ? 0 : 1;
684  if (noisy() > 0)
685  cout << "Restarting the classical continuation" << endl;
686  } else break;
687  } else break;
688  } while (!new_point);
689 
690  if (new_point) {
691  copy(X, Gamma, x, gamma);
692  copy(tX, tGamma, tx, tgamma);
693  } else if (this->non_smooth && this->singularities > 1) {
694  h0 = h_min();
695  treat_nonsmooth_point(x, gamma, tx0, tgamma0, true);
696  if (gmm::vect_size(get_x_next()) > 0) {
697  if (test_limit_point(tGamma)) {
698  set_sing_label("limit point");
699  if (noisy() > 0) cout << "Limit point detected!" << endl;
700  }
701  if (noisy() > 1)
702  cout << "Computing a test function for bifurcations"
703  << endl;
704  bool bifurcation_detected = (nb_tangent_sing() > 2);
705  if (bifurcation_detected) {
706  // update the stored values of the test function only
707  set_tau_bp_1(tau_bp_init);
708  set_tau_bp_2(test_function_bp(get_x_next(),
709  get_gamma_next(),
710  get_tx_sing(1),
711  get_tgamma_sing(1)));
712  } else
713  bifurcation_detected
714  = test_nonsmooth_bifurcation(x, gamma, tx, tgamma,
715  get_x_next(),
716  get_gamma_next(),
717  get_tx_sing(1),
718  get_tgamma_sing(1));
719  if (bifurcation_detected) {
720  set_sing_label("non-smooth bifurcation point");
721  if (noisy() > 0)
722  cout << "Non-smooth bifurcation point detected!" << endl;
723  }
724 
725  copy(get_x_next(), get_gamma_next(), x, gamma);
726  copy(get_tx_sing(1), get_tgamma_sing(1), tx, tgamma);
727  h = h_init();
728  new_point = true;
729  }
730  }
731 
732  if (!new_point) {
733  cout << "Continuation has failed!" << endl;
734  h0 = h = 0;
735  }
736  }
737 
738  void Moore_Penrose_continuation(VECT &x, double &gamma,
739  VECT &tx, double &tgamma, double &h) {
740  double h0;
741  Moore_Penrose_continuation(x, gamma, tx, tgamma, h, h0);
742  }
743 
744  protected:
745  // Linear algebra functions
746  void copy(const VECT &v1, const double &a1, VECT &v, double &a) const
747  { gmm::copy(v1, v); a = a1; }
748  void scale(VECT &v, double &a, double c) const { gmm::scale(v, c); a *= c; }
749  void scaled_add(const VECT &v1, const VECT &v2, double c2, VECT &v) const
750  { gmm::add(v1, gmm::scaled(v2, c2), v); }
751  void scaled_add(const VECT &v1, double c1,
752  const VECT &v2, double c2, VECT &v) const
753  { gmm::add(gmm::scaled(v1, c1), gmm::scaled(v2, c2), v); }
754  void scaled_add(const VECT &v1, const double &a1,
755  const VECT &v2, const double &a2, double c2,
756  VECT &v, double &a) const
757  { gmm::add(v1, gmm::scaled(v2, c2), v); a = a1 + c2*a2; }
758  void scaled_add(const VECT &v1, const double &a1, double c1,
759  const VECT &v2, const double &a2, double c2,
760  VECT &v, double &a) const {
761  gmm::add(gmm::scaled(v1, c1), gmm::scaled(v2, c2), v);
762  a = c1*a1 + c2*a2;
763  }
764  void scaled_add(const MAT &m1, double c1,
765  const MAT &m2, double c2, MAT &m) const
766  { gmm::add(gmm::scaled(m1, c1), gmm::scaled(m2, c2), m); }
767  void mult(const MAT &A, const VECT &v1, VECT &v) const
768  { gmm::mult(A, v1, v); }
769 
770  double norm(const VECT &v) const
771  { return gmm::vect_norm2(v); }
772 
773  double sp(const VECT &v1, const VECT &v2) const
774  { return gmm::vect_sp(v1, v2); }
775  double sp(const VECT &v1, const VECT &v2, double w1, double w2) const
776  { return sp(v1, v2) + w1 * w2; }
777 
778  virtual double intrv_sp(const VECT &v1, const VECT &v2) const = 0;
779 
780  double w_sp(const VECT &v1, const VECT &v2) const
781  { return scfac * intrv_sp(v1, v2); }
782  double w_sp(const VECT &v1, const VECT &v2, double w1, double w2) const
783  { return w_sp(v1, v2) + w1 * w2; }
784  double w_norm(const VECT &v, double w) const
785  { return sqrt(w_sp(v, v) + w * w); }
786 
787  double cosang(const VECT &v1, const VECT &v2) const {
788  double no = sqrt(intrv_sp(v1, v1) * intrv_sp(v2, v2));
789  return (no == 0) ? 0: intrv_sp(v1, v2) / no;
790  }
791  double cosang(const VECT &v1, const VECT &v2, double w1, double w2) const {
792 //CHANGE 3: new definition of cosang
793 //double wgamma(0.1);
794 //double no = w_norm(v1, wgamma*w1) * w_norm(v2, wgamma*w2);
795 //return (no == 0) ? 0 : w_sp(v1, v2, wgamma*w1, wgamma*w2) / no;
796  double no = sqrt((intrv_sp(v1, v1) + w1*w1)*
797  (intrv_sp(v2, v2) + w2*w2));
798  return (no == 0) ? 0 : (intrv_sp(v1, v2) + w1*w2) / no;
799  }
800 
801  public:
802 
803  // Misc. for accessing private data
804  int noisy(void) const { return noisy_; }
805  double h_init(void) const { return h_init_; }
806  double h_min(void) const { return h_min_; }
807  double h_max(void) const { return h_max_; }
808  double h_dec(void) const { return h_dec_; }
809  double h_inc(void) const { return h_inc_; }
810  size_type maxit(void) const { return maxit_; }
811  size_type thrit(void) const { return thrit_; }
812  double maxres(void) const { return maxres_; }
813  double maxdiff(void) const { return maxdiff_; }
814  double mincos(void) const { return mincos_; }
815  double delta_max(void) const { return delta_max_; }
816  double delta_min(void) const { return delta_min_; }
817  double thrvar(void) const { return thrvar_; }
818  size_type nbdir(void) const { return nbdir_; }
819  size_type nbspan(void) const { return nbspan_; }
820 
821  void set_tau_lp(double tau) { tau_lp = tau; }
822  double get_tau_lp(void) const { return tau_lp; }
823  void set_tau_bp_1(double tau) { tau_bp_1 = tau; }
824  double get_tau_bp_1(void) const { return tau_bp_1; }
825  void set_tau_bp_2(double tau) { tau_bp_2 = tau; }
826  double get_tau_bp_2(void) const { return tau_bp_2; }
827  void clear_tau_bp_currentstep(void) {
828  tau_bp_graph.clear();
829  }
830  void init_tau_bp_graph(void) { tau_bp_graph[0.] = tau_bp_2; }
831  void insert_tau_bp_graph(double alpha, double tau) {
832  tau_bp_graph[alpha] = tau;
833  gmm::resize(alpha_hist, 0);
834  gmm::resize(tau_bp_hist, 0);
835  }
836  const VECT &get_alpha_hist(void) {
837  size_type i = 0;
838  gmm::resize(alpha_hist, tau_bp_graph.size());
839  for (std::map<double, double>::iterator it = tau_bp_graph.begin();
840  it != tau_bp_graph.end(); it++) {
841  alpha_hist[i] = (*it).first; i++;
842  }
843  return alpha_hist;
844  }
845  const VECT &get_tau_bp_hist(void) {
846  size_type i = 0;
847  gmm::resize(tau_bp_hist, tau_bp_graph.size());
848  for (std::map<double, double>::iterator it = tau_bp_graph.begin();
849  it != tau_bp_graph.end(); it++) {
850  tau_bp_hist[i] = (*it).second; i++;
851  }
852  return tau_bp_hist;
853  }
854 
855  void clear_sing_data(void) {
856  sing_label = "";
857  gmm::resize(x_sing, 0);
858  gmm::resize(x_next, 0);
859  tx_sing.clear();
860  tgamma_sing.clear();
861  tx_predict.clear();
862  tgamma_predict.clear();
863  }
864  void set_sing_label(std::string label) { sing_label = label; }
865  const std::string get_sing_label(void) const { return sing_label; }
866  void set_sing_point(const VECT &x, double gamma) {
867  gmm::resize(x_sing, gmm::vect_size(x));
868  copy(x, gamma, x_sing, gamma_sing);
869  }
870  const VECT &get_x_sing(void) const { return x_sing; }
871  double get_gamma_sing(void) const { return gamma_sing; }
872  size_type nb_tangent_sing(void) const { return tx_sing.size(); }
873  bool insert_tangent_sing(const VECT &tx, double tgamma){
874  bool is_included = false;
875  for (size_type i = 0; (i < tx_sing.size()) && (!is_included); ++i) {
876  double cang = cosang(tx_sing[i], tx, tgamma_sing[i], tgamma);
877  is_included = (cang >= mincos_);
878  }
879  if (!is_included) {
880  tx_sing.push_back(tx);
881  tgamma_sing.push_back(tgamma);
882  }
883  return !is_included;
884  }
885  const VECT &get_tx_sing(size_type i) const { return tx_sing[i]; }
886  double get_tgamma_sing(size_type i) const { return tgamma_sing[i]; }
887  const std::vector<VECT> &get_tx_sing(void) const { return tx_sing; }
888  const std::vector<double> &get_tgamma_sing(void) const { return tgamma_sing; }
889 
890  void set_next_point(const VECT &x, double gamma) {
891  if (gmm::vect_size(x_next) == 0) {
892  gmm::resize(x_next, gmm::vect_size(x));
893  copy(x, gamma, x_next, gamma_next);
894  }
895  }
896  const VECT &get_x_next(void) const { return x_next; }
897  double get_gamma_next(void) const { return gamma_next; }
898 
899  bool insert_tangent_predict(const VECT &tx, double tgamma) {
900  bool is_included = false;
901  for (size_type i = 0; (i < tx_predict.size()) && (!is_included); ++i) {
902  double cang = gmm::abs(cosang(tx_predict[i], tx, tgamma_predict[i], tgamma));
903  is_included = (cang >= mincos_);
904  }
905  if (!is_included) {
906  tx_predict.push_back(tx);
907  tgamma_predict.push_back(tgamma);
908  }
909  return !is_included;
910  }
911 
912  void init_border(size_type nbdof) {
913  srand(unsigned(time(NULL)));
914  gmm::resize(bb_x_, nbdof); gmm::fill_random(bb_x_);
915  gmm::resize(cc_x_, nbdof); gmm::fill_random(cc_x_);
916  bb_gamma = gmm::random(1.)/scalar_type(nbdof);
917  cc_gamma = gmm::random(1.)/scalar_type(nbdof);
918  dd = gmm::random(1.)/scalar_type(nbdof);
919  gmm::scale(bb_x_, scalar_type(1)/scalar_type(nbdof));
920  gmm::scale(cc_x_, scalar_type(1)/scalar_type(nbdof));
921  }
922 
923  protected:
924 
925  const VECT &bb_x(size_type nbdof)
926  { if (gmm::vect_size(bb_x_) != nbdof) init_border(nbdof); return bb_x_; }
927  const VECT &cc_x(size_type nbdof)
928  { if (gmm::vect_size(cc_x_) != nbdof) init_border(nbdof); return cc_x_; }
929 
930  size_type estimated_memsize(void) {
931  size_type szd = sizeof(double);
932  return (this->singularities == 0) ? 0
933  : (2 * gmm::vect_size(bb_x_) * szd
934  + 4 * gmm::vect_size(get_tau_bp_hist()) * szd
935  + (1 + nb_tangent_sing()) * gmm::vect_size(get_x_sing()) * szd);
936  }
937 
938  // virtual methods
939 
940  // solve A * g = L
941  virtual void solve(const MAT &A, VECT &g, const VECT &L) const = 0;
942  // solve A * (g1|g2) = (L1|L2)
943  virtual void solve(const MAT &A, VECT &g1, VECT &g2,
944  const VECT &L1, const VECT &L2) const = 0;
945  // F(x, gamma) --> f
946  virtual void F(const VECT &x, double gamma, VECT &f) const = 0;
947  // (F(x, gamma + eps) - f0) / eps --> g
948  virtual void F_gamma(const VECT &x, double gamma, const VECT &f0,
949  VECT &g) const = 0;
950  // (F(x, gamma + eps) - F(x, gamma)) / eps --> g
951  virtual void F_gamma(const VECT &x, double gamma, VECT &g) const = 0;
952  // F_x(x, gamma) --> A
953  virtual void F_x(const VECT &x, double gamma, MAT &A) const = 0;
954  // solve F_x(x, gamma) * g = L
955  virtual void solve_grad(const VECT &x, double gamma,
956  VECT &g, const VECT &L) const = 0;
957  // solve F_x(x, gamma) * (g1|g2) = (L1|L2)
958  virtual void solve_grad(const VECT &x, double gamma, VECT &g1,
959  VECT &g2, const VECT &L1, const VECT &L2) const = 0;
960  // F_x(x, gamma) * w --> y
961  virtual void mult_grad(const VECT &x, double gamma,
962  const VECT &w, VECT &y) const = 0;
963 
964  public:
965 
966  virtual_cont_struct
967  (int sing = 0, bool nonsm = false, double sfac=0.,
968  double hin = 1.e-2, double hmax = 1.e-1, double hmin = 1.e-5,
969  double hinc = 1.3, double hdec = 0.5,
970  size_type mit = 10, size_type tit = 4,
971  double mres = 1.e-6, double mdiff = 1.e-6, double mcos = 0.9,
972  double dmax = 0.005, double dmin = 0.00012, double tvar = 0.02,
973  size_type ndir = 40, size_type nspan = 1, int noi = 0)
974  : singularities(sing), non_smooth(nonsm), scfac(sfac),
975  h_init_(hin), h_max_(hmax), h_min_(hmin), h_inc_(hinc), h_dec_(hdec),
976  maxit_(mit), thrit_(tit), maxres_(mres), maxdiff_(mdiff), mincos_(mcos),
977  delta_max_(dmax), delta_min_(dmin), thrvar_(tvar),
978  nbdir_(ndir), nbspan_(nspan), noisy_(noi),
979  tau_lp(0.), tau_bp_1(tau_bp_init), tau_bp_2(tau_bp_init),
980  gamma_sing(0.), gamma_next(0.)
981  {}
982  virtual ~virtual_cont_struct() {}
983 
984  };
985 
986 
987 
988 
989 
990 
991  //=========================================================================
992  // Moore-Penrose continuation method for Getfem models
993  //=========================================================================
994 
995 
996 #ifdef GETFEM_MODELS_H__
997 
998  class cont_struct_getfem_model
999  : public virtual_cont_struct<base_vector, model_real_sparse_matrix>,
1000  virtual public dal::static_stored_object {
1001 
1002  private:
1003  mutable model *md;
1004  std::string parameter_name;
1005  std::string initdata_name, finaldata_name, currentdata_name;
1006  gmm::sub_interval I; // for continuation based on a subset of model variables
1007  rmodel_plsolver_type lsolver;
1008  double maxres_solve;
1009 
1010  void set_variables(const base_vector &x, double gamma) const;
1011  void update_matrix(const base_vector &x, double gamma) const;
1012 
1013  // implemented virtual methods
1014 
1015  double intrv_sp(const base_vector &v1, const base_vector &v2) const {
1016  return (I.size() > 0) ? gmm::vect_sp(gmm::sub_vector(v1,I),
1017  gmm::sub_vector(v2,I))
1018  : gmm::vect_sp(v1, v2);
1019  }
1020 
1021  // solve A * g = L
1022  void solve(const model_real_sparse_matrix &A, base_vector &g, const base_vector &L) const;
1023  // solve A * (g1|g2) = (L1|L2)
1024  void solve(const model_real_sparse_matrix &A, base_vector &g1, base_vector &g2,
1025  const base_vector &L1, const base_vector &L2) const;
1026  // F(x, gamma) --> f
1027  void F(const base_vector &x, double gamma, base_vector &f) const;
1028  // (F(x, gamma + eps) - f0) / eps --> g
1029  void F_gamma(const base_vector &x, double gamma, const base_vector &f0,
1030  base_vector &g) const;
1031  // (F(x, gamma + eps) - F(x, gamma)) / eps --> g
1032  void F_gamma(const base_vector &x, double gamma, base_vector &g) const;
1033 
1034  // F_x(x, gamma) --> A
1035  void F_x(const base_vector &x, double gamma, model_real_sparse_matrix &A) const;
1036  // solve F_x(x, gamma) * g = L
1037  void solve_grad(const base_vector &x, double gamma,
1038  base_vector &g, const base_vector &L) const;
1039  // solve F_x(x, gamma) * (g1|g2) = (L1|L2)
1040  void solve_grad(const base_vector &x, double gamma, base_vector &g1,
1041  base_vector &g2,
1042  const base_vector &L1, const base_vector &L2) const;
1043  // F_x(x, gamma) * w --> y
1044  void mult_grad(const base_vector &x, double gamma,
1045  const base_vector &w, base_vector &y) const;
1046 
1047  public:
1048  size_type estimated_memsize(void);
1049  const model &linked_model(void) { return *md; }
1050 
1051  void set_parametrised_data_names
1052  (const std::string &in, const std::string &fn, const std::string &cn) {
1053  initdata_name = in;
1054  finaldata_name = fn;
1055  currentdata_name = cn;
1056  }
1057 
1058  void set_interval_from_variable_name(const std::string &varname) {
1059  if (varname == "") I = gmm::sub_interval(0,0);
1060  else I = md->interval_of_variable(varname);
1061  }
1062 
1063  cont_struct_getfem_model
1064  (model &md_, const std::string &pn, double sfac, rmodel_plsolver_type ls,
1065  double hin = 1.e-2, double hmax = 1.e-1, double hmin = 1.e-5,
1066  double hinc = 1.3, double hdec = 0.5, size_type mit = 10,
1067  size_type tit = 4, double mres = 1.e-6, double mdiff = 1.e-6,
1068  double mcos = 0.9, double mress = 1.e-8, int noi = 0, int sing = 0,
1069  bool nonsm = false, double dmax = 0.005, double dmin = 0.00012,
1070  double tvar = 0.02, size_type ndir = 40, size_type nspan = 1)
1071  : virtual_cont_struct(sing, nonsm, sfac, hin, hmax, hmin, hinc, hdec,
1072  mit, tit, mres, mdiff, mcos, dmax, dmin, tvar,
1073  ndir, nspan, noi),
1074  md(&md_), parameter_name(pn),
1075  initdata_name(""), finaldata_name(""), currentdata_name(""),
1076  I(0,0), lsolver(ls), maxres_solve(mress)
1077  {
1078  GMM_ASSERT1(!md->is_complex(),
1079  "Continuation has only a real version, sorry.");
1080  }
1081 
1082  };
1083 
1084 #endif
1085 
1086 
1087 } /* end of namespace getfem. */
1088 
1089 
1090 #endif /* GETFEM_CONTINUATION_H__ */
getfem_model_solvers.h
Standard solvers for model bricks.
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::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:263
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
gmm::fill_random
void fill_random(L &l, double cfill)
*‍/
Definition: gmm_blas.h:154
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
dal::static_stored_object
base class for static stored objects
Definition: dal_static_stored_objects.h:206
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:977