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