GetFEM  5.4.4
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 <[email protected]>
34  @author Yves Renard <[email protected]>
35  @author Konstantinos Poulios <[email protected]>
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);
236  double v_gamma;
237  return test_function_bp(A, g, tx, tgamma, v_x, v_gamma);
238  }
239 
240  /* Test function for bifurcation points for the gradient computed at
241  (x, gamma). */
242  double test_function_bp(const VECT &x, double gamma,
243  const VECT &tx, double tgamma,
244  VECT &v_x, double &v_gamma) {
245  MAT A; VECT g(x);
246  F_x(x, gamma, A);
247  F_gamma(x, gamma, g);
248  return test_function_bp(A, g, tx, tgamma, v_x, v_gamma);
249  }
250 
251  double test_function_bp(const VECT &x, double gamma,
252  const VECT &tx, double tgamma) {
253  VECT v_x(x);
254  double v_gamma;
255  return test_function_bp(x, gamma, tx, tgamma, v_x, v_gamma);
256  }
257 
258  public:
259  /* Test for smooth bifurcation points. */
260  bool test_smooth_bifurcation(const VECT &x, double gamma,
261  const VECT &tx, double tgamma) {
262  double tau_bp_1_old = get_tau_bp_1();
263  set_tau_bp_1(get_tau_bp_2());
264  set_tau_bp_2(test_function_bp(x, gamma, tx, tgamma));
265  return (get_tau_bp_2() * get_tau_bp_1() < 0) &&
266  (gmm::abs(get_tau_bp_1()) < gmm::abs(tau_bp_1_old));
267  }
268 
269  /* Test for non-smooth bifurcation points. */
270  bool test_nonsmooth_bifurcation(const VECT &x1, double gamma1,
271  const VECT &tx1, double tgamma1,
272  const VECT &x2, double gamma2,
273  const VECT &tx2, double tgamma2) {
274  VECT g1(x1), g2(x1), g(x1), tx(x1);
275 
276  // compute gradients at the two given points
277  MAT A1, A2;
278  F_x(x2, gamma2, A2);
279  F_gamma(x2, gamma2, g2);
280  F_x(x1, gamma1, A1);
281  F_gamma(x1, gamma1, g1);
282  double tau1 = test_function_bp(A1, g1, tx1, tgamma1);
283  double tau2 = test_function_bp(A2, g2, tx2, tgamma2);
284  double tau_var_ref = std::max(gmm::abs(tau2 - tau1), 1.e-8);
285  set_tau_bp_2(tau1);
286  init_tau_bp_graph();
287  MAT A(A2);
288 
289  // monitor the sign changes of the test function on the convex
290  // combination
291  size_type nb_changes = 0;
292  double delta = delta_min(), tau0 = tau_bp_init, tgamma;
293  for (double alpha=0.; alpha < 1.; ) {
294  alpha = std::min(alpha + delta, 1.);
295  scaled_add(A1, 1.-alpha, A2, alpha, A); // A = (1-alpha)*A1 + alpha*A2
296  scaled_add(g1, 1.-alpha, g2, alpha, g); // g = (1-alpha)*g1 + alpha*g2
297  scaled_add(tx1, tgamma1, 1.-alpha, tx2, tgamma2, alpha, tx, tgamma);
298  //[tx,tgamma] = (1-alpha)*[tx1,tgamma1] + alpha*[tx2,tgamma2]
299 
300  tau2 = test_function_bp(A, g, tx, tgamma);
301  if ((tau2 * tau1 < 0) && (gmm::abs(tau1) < gmm::abs(tau0)))
302  ++nb_changes;
303  insert_tau_bp_graph(alpha, tau2);
304 
305  if (gmm::abs(tau2 - tau1) < 0.5 * thrvar() * tau_var_ref)
306  delta = std::min(2 * delta, delta_max());
307  else if (gmm::abs(tau2 - tau1) > thrvar() * tau_var_ref)
308  delta = std::max(0.1 * delta, delta_min());
309  tau0 = tau1;
310  tau1 = tau2;
311  }
312 
313  set_tau_bp_1(tau_bp_init);
314  set_tau_bp_2(tau2);
315  return nb_changes % 2;
316  }
317 
318  private:
319  /* Newton-type corrections for the couple ((X, Gamma), (tX, tGamma)).
320  The current direction of (tX, tGamma) is informatively compared with
321  (tx, tgamma). */
322  bool newton_corr(VECT &X, double &Gamma, VECT &tX,
323  double &tGamma, const VECT &tx, double tgamma,
324  size_type &it) {
325  bool converged = false;
326  double Delta_Gamma, res(0), diff;
327  VECT f(X), g(X), Delta_X(X), y(X);
328 
329  if (noisy() > 1) cout << "Starting correction" << endl;
330  F(X, Gamma, f); // f = F(X, Gamma) = -rhs(X, Gamma)
331 //CHANGE 1: line search
332 //double res0 = norm(f);
333 
334  for (it=0; it < maxit() && res < 1.e8; ++it) {
335  F_gamma(X, Gamma, f, g); // g = F_gamma(X, Gamma)
336  solve_grad(X, Gamma, Delta_X, y, f, g); // y = F_x(X, Gamma)^-1 * g
337  // Delta_X = F_x(X, Gamma)^-1 * f
338  Delta_Gamma = sp(tX, Delta_X) / (sp(tX, y) - tGamma); // Delta_Gamma = tX.Delta_X / (tX.y - tGamma)
339  if (isnan(Delta_Gamma)) {
340  if (noisy() > 1) cout << "Newton correction failed with NaN" << endl;
341  return false;
342  }
343  gmm::add(gmm::scaled(y, -Delta_Gamma), Delta_X); // Delta_X -= Delta_Gamma * y
344  scaled_add(X, Gamma, Delta_X, Delta_Gamma, -1,
345  X, Gamma); // [X,Gamma] -= [Delta_X,Delta_Gamma]
346  F(X, Gamma, f); // f = F(X, gamma) = -rhs(X, gamma)
347  res = norm(f);
348 
349 //CHANGE 1: line search
350 //for (size_type ii=0; ii < 4 && (isnan(res) || res > res0); ++ii) { // some basic linesearch
351 // scale(Delta_X, Delta_Gamma, 0.5);
352 // scaled_add(X, Gamma, Delta_X, Delta_Gamma, 1, X, Gamma); // [X,Gamma] += [Delta_X,Delta_Gamma]
353 // F(X, Gamma, f); // f = F(X, gamma) = -rhs(X, gamma)
354 // res = norm(f);
355 //}
356 
357  tGamma = 1. / (tGamma - w_sp(tX, y)); // tGamma = 1 / (tGamma - k*tX.y)
358  gmm::copy(gmm::scaled(y, -tGamma), tX); // tX = -tGamma * y
359  scale(tX, tGamma, 1./w_norm(tX, tGamma)); // [tX,tGamma] /= w_norm(tX,tGamma)
360 
361  diff = w_norm(Delta_X, Delta_Gamma);
362  if (noisy() > 1)
363  cout << " Correction " << std::setw(3) << it << ":"
364  << " Gamma = " << std::fixed << std::setprecision(6) << Gamma
365  << " residual = " << std::scientific << std::setprecision(3) << res
366  << " difference = " << std::scientific << std::setprecision(3) << diff
367  << " cosang = " << std::fixed << std::setprecision(6)
368  << cosang(tX, tx, tGamma, tgamma) << endl;
369 
370  if (res <= maxres() && diff <= maxdiff()) {
371  converged = true;
372  // recalculate the final tangent, for sure
373  compute_tangent(X, Gamma, tX, tGamma);
374  break;
375  }
376  }
377  if (noisy() > 1) cout << "Correction finished with Gamma = "
378  << Gamma << endl;
379  return converged;
380  }
381 
382  bool newton_corr(VECT &X, double &Gamma, VECT &tX,
383  double &tGamma, const VECT &tx, double tgamma) {
384  size_type it;
385  return newton_corr(X, Gamma, tX, tGamma, tx, tgamma, it);
386  }
387 
388  /* Try to perform one predictor-corrector step starting from the couple
389  ((x, gamma), (tx, tgamma)). Return the resulting couple in the case of
390  convergence. */
391  bool test_predict_dir(VECT &x, double &gamma,
392  VECT &tx, double &tgamma) {
393  bool converged = false;
394  double h = h_init(), Gamma, tGamma;
395  VECT X(x), tX(x);
396  while (!converged) { //step control
397  // prediction
398  scaled_add(x, gamma, tx, tgamma, h, X, Gamma); // [X,Gamma] = [x,gamma] + h * [tx,tgamma]
399  if (noisy() > 1)
400  cout << "(TPD) Prediction : Gamma = " << Gamma
401  << " (for h = " << h << ", tgamma = " << tgamma << ")" << endl;
402  copy(tx, tgamma, tX, tGamma);
403  //correction
404  converged = newton_corr(X, Gamma, tX, tGamma, tx, tgamma);
405 
406  if (h > h_min())
407  h = std::max(0.199 * h_dec() * h, h_min());
408  else
409  break;
410  }
411  if (converged) {
412  // check the direction of the tangent found
413  scaled_add(X, Gamma, x, gamma, -1., tx, tgamma); // [tx,tgamma] = [X,Gamma] - [x,gamma]
414  if (sp(tX, tx, tGamma, tgamma) < 0)
415  scale(tX, tGamma, -1.); // [tX,tGamma] *= -1
416  copy(X, Gamma, x, gamma);
417  copy(tX, tGamma, tx, tgamma);
418  }
419  return converged;
420  }
421 
422  /* A tool for approximating a smooth bifurcation point close to (x, gamma)
423  and locating the two branches emanating from there. */
424  void treat_smooth_bif_point(const VECT &x, double gamma,
425  const VECT &tx, double tgamma, double h) {
426  double tau0(get_tau_bp_1()), tau1(get_tau_bp_2());
427  double gamma0(gamma), Gamma,
428  tgamma0(tgamma), tGamma(tgamma), v_gamma;
429  VECT x0(x), X(x), tx0(tx), tX(tx), v_x(tx);
430 
431  if (noisy() > 0)
432  cout << "Starting locating a bifurcation point" << endl;
433 
434  // predictor-corrector steps with a secant-type step-length adaptation
435  h *= tau1 / (tau0 - tau1);
436  for (size_type i=0; i < 10 && (gmm::abs(h) >= h_min()); ++i) {
437  scaled_add(x0, gamma0, tx0, tgamma0, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h * [tx0,tgamma0]
438  if (noisy() > 1)
439  cout << "(TSBP) Prediction : Gamma = " << Gamma
440  << " (for h = " << h << ", tgamma = " << tgamma << ")" << endl;
441  if (newton_corr(X, Gamma, tX, tGamma, tx0, tgamma0)) {
442  copy(X, Gamma, x0, gamma0);
443  if (cosang(tX, tx0, tGamma, tgamma0) >= mincos())
444  copy(tX, tGamma, tx0, tgamma0);
445  tau0 = tau1;
446  tau1 = test_function_bp(X, Gamma, tx0, tgamma0, v_x, v_gamma);
447  h *= tau1 / (tau0 - tau1);
448  } else {
449  scaled_add(x0, gamma0, tx0, tgamma0, h, x0, gamma0); // [x0,gamma0] += h*[tx0,tgamma0]
450  test_function_bp(x0, gamma0, tx0, tgamma0, v_x, v_gamma);
451  break;
452  }
453  }
454  if (noisy() > 0)
455  cout << "Bifurcation point located" << endl;
456  set_sing_point(x0, gamma0);
457  insert_tangent_sing(tx0, tgamma0);
458 
459  if (noisy() > 0)
460  cout << "Starting searching for the second branch" << endl;
461  double no = w_norm(v_x, v_gamma);
462  scale(v_x, v_gamma, 1./no); // [v_x,v_gamma] /= no
463  if (test_predict_dir(x0, gamma0, v_x, v_gamma)
464  && insert_tangent_sing(v_x, v_gamma))
465  { if (noisy() > 0) cout << "Second branch found" << endl; }
466  else if (noisy() > 0) cout << "Second branch not found!" << endl;
467  }
468 
469  public:
470 
471  /* A tool for approximating a non-smooth point close to (x, gamma) and
472  consequent locating one-sided smooth solution branches emanating from
473  there. It is supposed that (x, gamma) is a point on a smooth solution
474  branch within the distance of h_min() from the end point of this
475  branch and (tx, tgamma) is the corresponding tangent directed towards
476  the end point. The boolean set_next determines whether the first new
477  branch found (if any) is to be chosen for further continuation. */
478  void treat_nonsmooth_point(const VECT &x, double gamma,
479  const VECT &tx, double tgamma, bool set_next) {
480  double gamma0(gamma), Gamma(gamma);
481  double tgamma0(tgamma), tGamma(tgamma);
482  double h = h_min(), mcos = mincos();
483  VECT x0(x), X(x), tx0(tx), tX(tx);
484 
485  // approximate the end point more precisely by a bisection-like algorithm
486  if (noisy() > 0)
487  cout << "Starting locating a non-smooth point" << endl;
488 
489  scaled_add(x0, gamma0, tx0, tgamma0, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h*[tx0,tgamma0]
490  if (newton_corr(X, Gamma, tX, tGamma, tx0, tgamma0)) { // --> X, Gamma, tX, tGamma
491  double cang = cosang(tX, tx0, tGamma, tgamma0);
492  if (cang >= mcos) mcos = (cang + 1.) / 2.;
493  }
494 
495  copy(tx0, tgamma0, tX, tGamma);
496  h /= 2.;
497  for (size_type i = 0; i < 15; i++) {
498  scaled_add(x0, gamma0, tx0, tgamma0, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h*[tx0,tgamma0]
499  if (noisy() > 1)
500  cout << "(TNSBP) Prediction : Gamma = " << Gamma
501  << " (for h = " << h << ", tgamma = " << tgamma << ")" << endl;
502  if (newton_corr(X, Gamma, tX, tGamma, tx0, tgamma0)
503  && (cosang(tX, tx0, tGamma, tgamma0) >= mcos)) {
504  copy(X, Gamma, x0, gamma0);
505  copy(tX, tGamma, tx0, tgamma0);
506  } else {
507  copy(tx0, tgamma0, tX, tGamma);
508  }
509  h /= 2.;
510  }
511  if (noisy() > 0)
512  cout << "Non-smooth point located" << endl;
513  set_sing_point(x0, gamma0);
514 
515  // take two reference vectors to span a subspace of directions emanating
516  // from the end point
517  if (noisy() > 0)
518  cout << "Starting a thorough search for different branches" << endl;
519  double tgamma1 = tgamma0, tgamma2 = tgamma0;
520  VECT tx1(tx0), tx2(tx0);
521  scale(tx1, tgamma1, -1.); // [tx1,tgamma1] *= -1
522  insert_tangent_sing(tx1, tgamma1);
523  h = h_min();
524  scaled_add(x0, gamma0, tx0, tgamma0, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h*[tx0,tgamma0]
525  compute_tangent(X, Gamma, tx2, tgamma2);
526 
527  // try systematically the directions of linear combinations of the couple
528  // of the reference vectors for finding new possible tangent predictions
529  // emanating from the end point
530  size_type i1 = 0, i2 = 0, nspan = 0;
531  double a, a1, a2, no;
532 
533  do {
534  for (size_type i = 0; i < nbdir(); i++) {
535  a = (2 * M_PI * double(i)) / double(nbdir());
536  a1 = sin(a);
537  a2 = cos(a);
538  scaled_add(tx1, tgamma1, a1, tx2, tgamma2, a2, tX, tGamma); // [tX,tGamma] = a1*[tx1,tgamma1] + a2*[tx2,tgamma2]
539  no = w_norm(tX, tGamma);
540  scaled_add(x0, gamma0, tX, tGamma, h/no, X, Gamma); // [X,Gamma] = [x0,gamma0] + h/no * [tX,tGamma]
541  compute_tangent(X, Gamma, tX, tGamma);
542 
543  if (gmm::abs(cosang(tX, tx0, tGamma, tgamma0)) < mincos()
544  || (i == 0 && nspan == 0)) {
545  copy(tX, tGamma, tx0, tgamma0);
546  if (insert_tangent_predict(tX, tGamma)) {
547  if (noisy() > 1)
548  cout << "New potential tangent vector found, "
549  << "trying one predictor-corrector step" << endl;
550  copy(x0, gamma0, X, Gamma);
551 
552  if (test_predict_dir(X, Gamma, tX, tGamma)) {
553  if (insert_tangent_sing(tX, tGamma)) {
554  if ((i == 0) && (nspan == 0)
555  // => (tX, tGamma) = (tx2, tgamma2)
556  && (gmm::abs(cosang(tX, tx0, tGamma, tgamma0))
557  >= mincos())) { i2 = 1; }
558  if (noisy() > 0) cout << "A new branch located (for nspan = "
559  << nspan << ")" << endl;
560  if (set_next) set_next_point(X, Gamma);
561 
562  }
563  copy(x0, gamma0, X, Gamma);
564  copy(tx0, tgamma0, tX, tGamma);
565  }
566 
567  scale(tX, tGamma, -1.); // [tX,tGamma] *= -1
568  if (test_predict_dir(X, Gamma, tX, tGamma)
569  && insert_tangent_sing(tX, tGamma)) {
570  if (noisy() > 0) cout << "A new branch located (for nspan = "
571  << nspan << ")" << endl;
572  if (set_next) set_next_point(X, Gamma);
573  }
574  }
575  }
576  }
577 
578  // heuristics for varying the reference vectors
579  bool perturb = true;
580  if (i1 + 1 < i2) { ++i1; perturb = false; }
581  else if(i2 + 1 < nb_tangent_sing())
582  { ++i2; i1 = 0; perturb = false; }
583  if (!perturb) {
584  copy(get_tx_sing(i1), get_tgamma_sing(i1), tx1, tgamma1);
585  copy(get_tx_sing(i2), get_tgamma_sing(i2), tx2, tgamma2);
586  } else {
587  gmm::fill_random(tX);
588  tGamma = gmm::random(1.);
589  no = w_norm(tX, tGamma);
590  scaled_add(tx2, tgamma2, tX, tGamma, 0.1/no, tx2, tgamma2);
591  // [tx2,tgamma2] += 0.1/no * [tX,tGamma]
592  scaled_add(x0, gamma0, tx2, tgamma2, h, X, Gamma); // [X,Gamma] = [x0,gamma0] + h*[tx2,tgamma2]
593  compute_tangent(X, Gamma, tx2, tgamma2);
594  }
595  } while (++nspan < nbspan());
596 
597  if (noisy() > 0)
598  cout << "Number of branches emanating from the non-smooth point "
599  << nb_tangent_sing() << endl;
600  }
601 
602 
603  void init_Moore_Penrose_continuation(const VECT &x,
604  double gamma, VECT &tx,
605  double &tgamma, double &h) {
606  gmm::clear(tx);
607  tgamma = (tgamma >= 0) ? 1. : -1.;
608  if (noisy() > 1)
609  cout << "Computing an initial tangent" << endl;
610  compute_tangent(x, gamma, tx, tgamma);
611  h = h_init();
612  if (this->singularities > 0)
613  init_test_functions(x, gamma, tx, tgamma);
614  }
615 
616 
617  /* Perform one step of the (non-smooth) Moore-Penrose continuation.
618  NOTE: The new point need not to be saved in the model in the end! */
619  void Moore_Penrose_continuation(VECT &x, double &gamma,
620  VECT &tx, double &tgamma,
621  double &h, double &h0) {
622  bool converged, new_point = false, tangent_switched = false;
623  size_type it, step_dec = 0;
624  double tgamma0 = tgamma, Gamma, tGamma;
625  VECT tx0(tx), X(x), tX(x);
626 
627  clear_tau_bp_currentstep();
628  clear_sing_data();
629 
630  do {
631  h0 = h;
632  // prediction
633  scaled_add(x, gamma, tx, tgamma, h, X, Gamma); // [X,Gamma] = [x,gamma] + h*[tx,tgamma]
634  if (noisy() > 1)
635  cout << " Prediction : Gamma = " << Gamma
636  << " (for h = " << std::scientific << std::setprecision(3) << h
637  << ", tgamma = " << tgamma << ")" << endl;
638  copy(tx, tgamma, tX, tGamma);
639 
640  // correction
641  converged = newton_corr(X, Gamma, tX, tGamma, tx, tgamma, it);
642  double cang(converged ? cosang(tX, tx, tGamma, tgamma) : 0);
643 
644  if (converged && cang >= mincos()) {
645  new_point = true;
646  if (this->singularities > 0) {
647  if (test_limit_point(tGamma)) {
648  set_sing_label("limit point");
649  if (noisy() > 0) cout << "Limit point detected!" << endl;
650  }
651  if (this->singularities > 1) { // Treat bifurcations
652  if (noisy() > 1)
653  cout << "New point found, computing a test function "
654  << "for bifurcations" << endl;
655  if (!tangent_switched) {
656  if (test_smooth_bifurcation(X, Gamma, tX, tGamma)) {
657  set_sing_label("smooth bifurcation point");
658  if (noisy() > 0)
659  cout << "Smooth bifurcation point detected!" << endl;
660  treat_smooth_bif_point(X, Gamma, tX, tGamma, h);
661  }
662  } else if (test_nonsmooth_bifurcation(x, gamma, tx0,
663  tgamma0, X, Gamma, tX,
664  tGamma)) {
665  set_sing_label("non-smooth bifurcation point");
666  if (noisy() > 0)
667  cout << "Non-smooth bifurcation point detected!" << endl;
668  treat_nonsmooth_point(x, gamma, tx0, tgamma0, false);
669  }
670  }
671  }
672 
673 //CHANGE 2: avoid false step increases
674 //if (step_dec == 0 && it < thrit() && h_inc()*(1-cang) < (1-mincos()))
675  if (step_dec == 0 && it < thrit())
676  h = std::min(h_inc() * h, h_max());
677  } else if (h > h_min()) {
678  h = std::max(h_dec() * h, h_min());
679  step_dec++;
680  } else if (this->non_smooth && !tangent_switched) {
681  if (noisy() > 0)
682  cout << "Classical continuation has failed!" << endl;
683  if (switch_tangent(x, gamma, tx, tgamma, h)) {
684  tangent_switched = true;
685  step_dec = (h >= h_init()) ? 0 : 1;
686  if (noisy() > 0)
687  cout << "Restarting the classical continuation" << endl;
688  } else break;
689  } else break;
690  } while (!new_point);
691 
692  if (new_point) {
693  copy(X, Gamma, x, gamma);
694  copy(tX, tGamma, tx, tgamma);
695  } else if (this->non_smooth && this->singularities > 1) {
696  h0 = h_min();
697  treat_nonsmooth_point(x, gamma, tx0, tgamma0, true);
698  if (gmm::vect_size(get_x_next()) > 0) {
699  if (test_limit_point(tGamma)) {
700  set_sing_label("limit point");
701  if (noisy() > 0) cout << "Limit point detected!" << endl;
702  }
703  if (noisy() > 1)
704  cout << "Computing a test function for bifurcations"
705  << endl;
706  bool bifurcation_detected = (nb_tangent_sing() > 2);
707  if (bifurcation_detected) {
708  // update the stored values of the test function only
709  set_tau_bp_1(tau_bp_init);
710  set_tau_bp_2(test_function_bp(get_x_next(),
711  get_gamma_next(),
712  get_tx_sing(1),
713  get_tgamma_sing(1)));
714  } else
715  bifurcation_detected
716  = test_nonsmooth_bifurcation(x, gamma, tx, tgamma,
717  get_x_next(),
718  get_gamma_next(),
719  get_tx_sing(1),
720  get_tgamma_sing(1));
721  if (bifurcation_detected) {
722  set_sing_label("non-smooth bifurcation point");
723  if (noisy() > 0)
724  cout << "Non-smooth bifurcation point detected!" << endl;
725  }
726 
727  copy(get_x_next(), get_gamma_next(), x, gamma);
728  copy(get_tx_sing(1), get_tgamma_sing(1), tx, tgamma);
729  h = h_init();
730  new_point = true;
731  }
732  }
733 
734  if (!new_point) {
735  cout << "Continuation has failed!" << endl;
736  h0 = h = 0;
737  }
738  }
739 
740  void Moore_Penrose_continuation(VECT &x, double &gamma,
741  VECT &tx, double &tgamma, double &h) {
742  double h0;
743  Moore_Penrose_continuation(x, gamma, tx, tgamma, h, h0);
744  }
745 
746  protected:
747  // Linear algebra functions
748  void copy(const VECT &v1, const double &a1, VECT &v, double &a) const
749  { gmm::copy(v1, v); a = a1; }
750  void scale(VECT &v, double &a, double c) const { gmm::scale(v, c); a *= c; }
751  void scaled_add(const VECT &v1, const VECT &v2, double c2, VECT &v) const
752  { gmm::add(v1, gmm::scaled(v2, c2), v); }
753  void scaled_add(const VECT &v1, double c1,
754  const VECT &v2, double c2, VECT &v) const
755  { gmm::add(gmm::scaled(v1, c1), gmm::scaled(v2, c2), v); }
756  void scaled_add(const VECT &v1, const double &a1,
757  const VECT &v2, const double &a2, double c2,
758  VECT &v, double &a) const
759  { gmm::add(v1, gmm::scaled(v2, c2), v); a = a1 + c2*a2; }
760  void scaled_add(const VECT &v1, const double &a1, double c1,
761  const VECT &v2, const double &a2, double c2,
762  VECT &v, double &a) const {
763  gmm::add(gmm::scaled(v1, c1), gmm::scaled(v2, c2), v);
764  a = c1*a1 + c2*a2;
765  }
766  void scaled_add(const MAT &m1, double c1,
767  const MAT &m2, double c2, MAT &m) const
768  { gmm::add(gmm::scaled(m1, c1), gmm::scaled(m2, c2), m); }
769  void mult(const MAT &A, const VECT &v1, VECT &v) const
770  { gmm::mult(A, v1, v); }
771 
772  double norm(const VECT &v) const
773  { return gmm::vect_norm2(v); }
774 
775  double sp(const VECT &v1, const VECT &v2) const
776  { return gmm::vect_sp(v1, v2); }
777  double sp(const VECT &v1, const VECT &v2, double w1, double w2) const
778  { return sp(v1, v2) + w1 * w2; }
779 
780  virtual double intrv_sp(const VECT &v1, const VECT &v2) const = 0;
781 
782  double w_sp(const VECT &v1, const VECT &v2) const
783  { return scfac * intrv_sp(v1, v2); }
784  double w_sp(const VECT &v1, const VECT &v2, double w1, double w2) const
785  { return w_sp(v1, v2) + w1 * w2; }
786  double w_norm(const VECT &v, double w) const
787  { return sqrt(w_sp(v, v) + w * w); }
788 
789  double cosang(const VECT &v1, const VECT &v2) const {
790  double no = sqrt(intrv_sp(v1, v1) * intrv_sp(v2, v2));
791  return (no == 0) ? 0: intrv_sp(v1, v2) / no;
792  }
793  double cosang(const VECT &v1, const VECT &v2, double w1, double w2) const {
794 //CHANGE 3: new definition of cosang
795 //double wgamma(0.1);
796 //double no = w_norm(v1, wgamma*w1) * w_norm(v2, wgamma*w2);
797 //return (no == 0) ? 0 : w_sp(v1, v2, wgamma*w1, wgamma*w2) / no;
798  double no = sqrt((intrv_sp(v1, v1) + w1*w1)*
799  (intrv_sp(v2, v2) + w2*w2));
800  return (no == 0) ? 0 : (intrv_sp(v1, v2) + w1*w2) / no;
801  }
802 
803  public:
804 
805  // Misc. for accessing private data
806  int noisy() const { return noisy_; }
807  double h_init() const { return h_init_; }
808  double h_min() const { return h_min_; }
809  double h_max() const { return h_max_; }
810  double h_dec() const { return h_dec_; }
811  double h_inc() const { return h_inc_; }
812  size_type maxit() const { return maxit_; }
813  size_type thrit() const { return thrit_; }
814  double maxres() const { return maxres_; }
815  double maxdiff() const { return maxdiff_; }
816  double mincos() const { return mincos_; }
817  double delta_max() const { return delta_max_; }
818  double delta_min() const { return delta_min_; }
819  double thrvar() const { return thrvar_; }
820  size_type nbdir() const { return nbdir_; }
821  size_type nbspan() const { return nbspan_; }
822 
823  void set_tau_lp(double tau) { tau_lp = tau; }
824  double get_tau_lp() const { return tau_lp; }
825  void set_tau_bp_1(double tau) { tau_bp_1 = tau; }
826  double get_tau_bp_1() const { return tau_bp_1; }
827  void set_tau_bp_2(double tau) { tau_bp_2 = tau; }
828  double get_tau_bp_2() const { return tau_bp_2; }
829  void clear_tau_bp_currentstep() {
830  tau_bp_graph.clear();
831  }
832  void init_tau_bp_graph() { tau_bp_graph[0.] = tau_bp_2; }
833  void insert_tau_bp_graph(double alpha, double tau) {
834  tau_bp_graph[alpha] = tau;
835  gmm::resize(alpha_hist, 0);
836  gmm::resize(tau_bp_hist, 0);
837  }
838  const VECT &get_alpha_hist() {
839  size_type i = 0;
840  gmm::resize(alpha_hist, tau_bp_graph.size());
841  for (std::map<double, double>::iterator it = tau_bp_graph.begin();
842  it != tau_bp_graph.end(); it++) {
843  alpha_hist[i] = (*it).first; i++;
844  }
845  return alpha_hist;
846  }
847  const VECT &get_tau_bp_hist() {
848  size_type i = 0;
849  gmm::resize(tau_bp_hist, tau_bp_graph.size());
850  for (std::map<double, double>::iterator it = tau_bp_graph.begin();
851  it != tau_bp_graph.end(); it++) {
852  tau_bp_hist[i] = (*it).second; i++;
853  }
854  return tau_bp_hist;
855  }
856 
857  void clear_sing_data() {
858  sing_label = "";
859  gmm::resize(x_sing, 0);
860  gmm::resize(x_next, 0);
861  tx_sing.clear();
862  tgamma_sing.clear();
863  tx_predict.clear();
864  tgamma_predict.clear();
865  }
866  void set_sing_label(std::string label) { sing_label = label; }
867  const std::string get_sing_label() const { return sing_label; }
868  void set_sing_point(const VECT &x, double gamma) {
869  gmm::resize(x_sing, gmm::vect_size(x));
870  copy(x, gamma, x_sing, gamma_sing);
871  }
872  const VECT &get_x_sing() const { return x_sing; }
873  double get_gamma_sing() const { return gamma_sing; }
874  size_type nb_tangent_sing() const { return tx_sing.size(); }
875  bool insert_tangent_sing(const VECT &tx, double tgamma){
876  bool is_included = false;
877  for (size_type i = 0; (i < tx_sing.size()) && (!is_included); ++i) {
878  double cang = cosang(tx_sing[i], tx, tgamma_sing[i], tgamma);
879  is_included = (cang >= mincos_);
880  }
881  if (!is_included) {
882  tx_sing.push_back(tx);
883  tgamma_sing.push_back(tgamma);
884  }
885  return !is_included;
886  }
887  const VECT &get_tx_sing(size_type i) const { return tx_sing[i]; }
888  double get_tgamma_sing(size_type i) const { return tgamma_sing[i]; }
889  const std::vector<VECT> &get_tx_sing() const { return tx_sing; }
890  const std::vector<double> &get_tgamma_sing() const { return tgamma_sing; }
891 
892  void set_next_point(const VECT &x, double gamma) {
893  if (gmm::vect_size(x_next) == 0) {
894  gmm::resize(x_next, gmm::vect_size(x));
895  copy(x, gamma, x_next, gamma_next);
896  }
897  }
898  const VECT &get_x_next() const { return x_next; }
899  double get_gamma_next() const { return gamma_next; }
900 
901  bool insert_tangent_predict(const VECT &tx, double tgamma) {
902  bool is_included = false;
903  for (size_type i = 0; (i < tx_predict.size()) && (!is_included); ++i) {
904  double cang = gmm::abs(cosang(tx_predict[i], tx, tgamma_predict[i], tgamma));
905  is_included = (cang >= mincos_);
906  }
907  if (!is_included) {
908  tx_predict.push_back(tx);
909  tgamma_predict.push_back(tgamma);
910  }
911  return !is_included;
912  }
913 
914  void init_border(size_type nbdof) {
915  srand(unsigned(time(NULL)));
916  gmm::resize(bb_x_, nbdof); gmm::fill_random(bb_x_);
917  gmm::resize(cc_x_, nbdof); gmm::fill_random(cc_x_);
918  bb_gamma = gmm::random(1.)/scalar_type(nbdof);
919  cc_gamma = gmm::random(1.)/scalar_type(nbdof);
920  dd = gmm::random(1.)/scalar_type(nbdof);
921  gmm::scale(bb_x_, scalar_type(1)/scalar_type(nbdof));
922  gmm::scale(cc_x_, scalar_type(1)/scalar_type(nbdof));
923  }
924 
925  protected:
926 
927  const VECT &bb_x(size_type nbdof)
928  { if (gmm::vect_size(bb_x_) != nbdof) init_border(nbdof); return bb_x_; }
929  const VECT &cc_x(size_type nbdof)
930  { if (gmm::vect_size(cc_x_) != nbdof) init_border(nbdof); return cc_x_; }
931 
932  size_type estimated_memsize() {
933  size_type szd = sizeof(double);
934  return (this->singularities == 0) ? 0
935  : (2 * gmm::vect_size(bb_x_) * szd
936  + 4 * gmm::vect_size(get_tau_bp_hist()) * szd
937  + (1 + nb_tangent_sing()) * gmm::vect_size(get_x_sing()) * szd);
938  }
939 
940  // virtual methods
941 
942  // solve A * g = L
943  virtual void solve(const MAT &A, VECT &g, const VECT &L) const = 0;
944  // solve A * (g1|g2) = (L1|L2)
945  virtual void solve(const MAT &A, VECT &g1, VECT &g2,
946  const VECT &L1, const VECT &L2) const = 0;
947  // F(x, gamma) --> f
948  virtual void F(const VECT &x, double gamma, VECT &f) const = 0;
949  // (F(x, gamma + eps) - f0) / eps --> g
950  virtual void F_gamma(const VECT &x, double gamma, const VECT &f0,
951  VECT &g) const = 0;
952  // (F(x, gamma + eps) - F(x, gamma)) / eps --> g
953  virtual void F_gamma(const VECT &x, double gamma, VECT &g) const = 0;
954  // F_x(x, gamma) --> A
955  virtual void F_x(const VECT &x, double gamma, MAT &A) const = 0;
956  // solve F_x(x, gamma) * g = L
957  virtual void solve_grad(const VECT &x, double gamma,
958  VECT &g, const VECT &L) const = 0;
959  // solve F_x(x, gamma) * (g1|g2) = (L1|L2)
960  virtual void solve_grad(const VECT &x, double gamma, VECT &g1,
961  VECT &g2, const VECT &L1, const VECT &L2) const = 0;
962  // F_x(x, gamma) * w --> y
963  virtual void mult_grad(const VECT &x, double gamma,
964  const VECT &w, VECT &y) const = 0;
965 
966  public:
967 
968  virtual_cont_struct
969  (int sing = 0, bool nonsm = false, double sfac=0.,
970  double hin = 1.e-2, double hmax = 1.e-1, double hmin = 1.e-5,
971  double hinc = 1.3, double hdec = 0.5,
972  size_type mit = 10, size_type tit = 4,
973  double mres = 1.e-6, double mdiff = 1.e-6, double mcos = 0.9,
974  double dmax = 0.005, double dmin = 0.00012, double tvar = 0.02,
975  size_type ndir = 40, size_type nspan = 1, int noi = 0)
976  : singularities(sing), non_smooth(nonsm), scfac(sfac),
977  h_init_(hin), h_max_(hmax), h_min_(hmin), h_inc_(hinc), h_dec_(hdec),
978  maxit_(mit), thrit_(tit), maxres_(mres), maxdiff_(mdiff), mincos_(mcos),
979  delta_max_(dmax), delta_min_(dmin), thrvar_(tvar),
980  nbdir_(ndir), nbspan_(nspan), noisy_(noi),
981  tau_lp(0.), tau_bp_1(tau_bp_init), tau_bp_2(tau_bp_init),
982  gamma_sing(0.), gamma_next(0.)
983  {}
984  virtual ~virtual_cont_struct() {}
985 
986  };
987 
988 
989 
990 
991 
992 
993  //=========================================================================
994  // Moore-Penrose continuation method for Getfem models
995  //=========================================================================
996 
997 
998 #ifdef GETFEM_MODELS_H__
999 
1000  class cont_struct_getfem_model
1001  : public virtual_cont_struct<base_vector, model_real_sparse_matrix>,
1002  virtual public dal::static_stored_object {
1003 
1004  private:
1005  mutable model *md;
1006  std::string parameter_name;
1007  std::string initdata_name, finaldata_name, currentdata_name;
1008  gmm::sub_interval I; // for continuation based on a subset of model variables
1009  rmodel_plsolver_type lsolver;
1010  double maxres_solve;
1011 
1012  void set_variables(const base_vector &x, double gamma) const;
1013  void update_matrix(const base_vector &x, double gamma) const;
1014 
1015  // implemented virtual methods
1016 
1017  double intrv_sp(const base_vector &v1, const base_vector &v2) const {
1018  return (I.size() > 0) ? gmm::vect_sp(gmm::sub_vector(v1,I),
1019  gmm::sub_vector(v2,I))
1020  : gmm::vect_sp(v1, v2);
1021  }
1022 
1023  // solve A * g = L
1024  void solve(const model_real_sparse_matrix &A, base_vector &g, const base_vector &L) const;
1025  // solve A * (g1|g2) = (L1|L2)
1026  void solve(const model_real_sparse_matrix &A, base_vector &g1, base_vector &g2,
1027  const base_vector &L1, const base_vector &L2) const;
1028  // F(x, gamma) --> f
1029  void F(const base_vector &x, double gamma, base_vector &f) const;
1030  // (F(x, gamma + eps) - f0) / eps --> g
1031  void F_gamma(const base_vector &x, double gamma, const base_vector &f0,
1032  base_vector &g) const;
1033  // (F(x, gamma + eps) - F(x, gamma)) / eps --> g
1034  void F_gamma(const base_vector &x, double gamma, base_vector &g) const;
1035 
1036  // F_x(x, gamma) --> A
1037  void F_x(const base_vector &x, double gamma, model_real_sparse_matrix &A) const;
1038  // solve F_x(x, gamma) * g = L
1039  void solve_grad(const base_vector &x, double gamma,
1040  base_vector &g, const base_vector &L) const;
1041  // solve F_x(x, gamma) * (g1|g2) = (L1|L2)
1042  void solve_grad(const base_vector &x, double gamma, base_vector &g1,
1043  base_vector &g2,
1044  const base_vector &L1, const base_vector &L2) const;
1045  // F_x(x, gamma) * w --> y
1046  void mult_grad(const base_vector &x, double gamma,
1047  const base_vector &w, base_vector &y) const;
1048 
1049  public:
1050  size_type estimated_memsize();
1051  const model &linked_model() { return *md; }
1052 
1053  void set_parametrised_data_names
1054  (const std::string &in, const std::string &fn, const std::string &cn) {
1055  initdata_name = in;
1056  finaldata_name = fn;
1057  currentdata_name = cn;
1058  }
1059 
1060  void set_interval_from_variable_name(const std::string &varname) {
1061  if (varname == "") I = gmm::sub_interval(0,0);
1062  else I = md->interval_of_variable(varname);
1063  }
1064 
1065  cont_struct_getfem_model
1066  (model &md_, const std::string &pn, double sfac, rmodel_plsolver_type ls,
1067  double hin = 1.e-2, double hmax = 1.e-1, double hmin = 1.e-5,
1068  double hinc = 1.3, double hdec = 0.5, size_type mit = 10,
1069  size_type tit = 4, double mres = 1.e-6, double mdiff = 1.e-6,
1070  double mcos = 0.9, double mress = 1.e-8, int noi = 0, int sing = 0,
1071  bool nonsm = false, double dmax = 0.005, double dmin = 0.00012,
1072  double tvar = 0.02, size_type ndir = 40, size_type nspan = 1)
1073  : virtual_cont_struct(sing, nonsm, sfac, hin, hmax, hmin, hinc, hdec,
1074  mit, tit, mres, mdiff, mcos, dmax, dmin, tvar,
1075  ndir, nspan, noi),
1076  md(&md_), parameter_name(pn),
1077  initdata_name(""), finaldata_name(""), currentdata_name(""),
1078  I(0,0), lsolver(ls), maxres_solve(mress)
1079  {
1080  GMM_ASSERT1(!md->is_complex(),
1081  "Continuation has only a real version, sorry.");
1082  }
1083 
1084  };
1085 
1086 #endif
1087 
1088 
1089 } /* end of namespace getfem. */
1090 
1091 
1092 #endif /* GETFEM_CONTINUATION_H__ */
dal::static_stored_object
base class for static stored objects
Definition: dal_static_stored_objects.h:206
getfem_model_solvers.h
Standard solvers for model bricks.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:977
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)
fill a vector or matrix with random value (uniform [-1,1]).
Definition: gmm_blas.h:130
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:210
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1664
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:264
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1276
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
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
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:47