[Getfem-commits] r4872 - in /trunk/getfem: interface/src/gf_cont_struct.cc interface/src/gf_cont_struct_get.cc src/getfem/getfem_continuation.h

[tomas . ligursky]

Author: ligut2am
Date: Fri Mar  6 10:10:11 2015
New Revision: 4872

URL: http://svn.gna.org/viewcvs/getfem?rev=4872&view=rev
Log:
minor changes in numerical branching

Modified:
    trunk/getfem/interface/src/gf_cont_struct.cc
    trunk/getfem/interface/src/gf_cont_struct_get.cc
    trunk/getfem/src/getfem/getfem_continuation.h

Modified: trunk/getfem/interface/src/gf_cont_struct.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_cont_struct.cc?rev=4872&r1=4871&r2=4872&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_cont_struct.cc        (original)
+++ trunk/getfem/interface/src/gf_cont_struct.cc        Fri Mar  6 10:10:11 2015
@@ -1,6 +1,6 @@
 /*===========================================================================
  
- Copyright (C) 2012-2014 Tomas Ligursky, Yves Renard.
+ Copyright (C) 2012-2015 Tomas Ligursky, Yves Renard.
  
  This file is a part of GETFEM++
  
@@ -83,7 +83,8 @@
        default value is 1e-8);
     - 'singularities', @int SING
        activates tools for detection and treatment of singular points (1 for
-       limit points, 2 for bifurcation points);
+       limit points, 2 for bifurcation points and points requiring special
+       branching techniques);
     - 'non-smooth'
        determines that some special methods for non-smooth problems can be
        used;

Modified: trunk/getfem/interface/src/gf_cont_struct_get.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_cont_struct_get.cc?rev=4872&r1=4871&r2=4872&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_cont_struct_get.cc    (original)
+++ trunk/getfem/interface/src/gf_cont_struct_get.cc    Fri Mar  6 10:10:11 2015
@@ -1,6 +1,6 @@
 /*===========================================================================
  
- Copyright (C) 2012-2012 Tomas Ligursky, Yves Renard.
+ Copyright (C) 2012-2015 Tomas Ligursky, Yves Renard.
  
  This file is a part of GETFEM++
  
@@ -138,6 +138,7 @@
       with the tangent given by `tangent_sol2` and address@hidden/
     sub_command
       ("non-smooth bifurcation test", 8, 8, 0, 1,
+
        size_type nbdof = ps->linked_model().nb_dof();
        darray x1 = in.pop().to_darray();
        std::vector<double> xx1(nbdof); gmm::copy(x1, xx1);
@@ -151,6 +152,7 @@
        darray t_x2 = in.pop().to_darray();
        std::vector<double> tt_x2(nbdof); gmm::copy(t_x2, tt_x2);
        scalar_type t_gamma2 = in.pop().to_scalar();
+
        ps->set_build(getfem::BUILD_ALL);
        ps->init_border();
        ps->clear_tau_bp_currentstep();
@@ -166,18 +168,47 @@
       of address@hidden/
     sub_command
       ("bifurcation test function", 0, 0, 0, 3,
+
        out.pop().from_scalar(ps->get_tau_bp_2());
        if (out.remaining()) out.pop().from_dcvector(ps->get_alpha_hist());
        if (out.remaining()) out.pop().from_dcvector(ps->get_tau_bp_hist());
        );
 
 
+    /* @GET ('non-smooth branching', @vec solution, @scalar parameter, @vec 
tangent_sol, @scalar tangent_par)
+       Approximate a non-smooth point close to the point given by `solution`
+       and `parameter` and locate one-sided smooth solution branches
+       emanating from there. Save the approximation of the non-smooth point
+       as a singular point into the @tcs object together with the array of
+       the tangents to the located solution branches that direct away from
+       the non-smooth point. It is supposed that the point given by
+       `solution` and `parameter` is a point on a smooth solution branch
+       within the distance equal to the minimum step size from the end point
+       of this branch, and the corresponding tangent that is directed towards
+       the end point is given by `tangent_sol` and address@hidden/
+    sub_command
+      ("non-smooth branching", 4, 4, 0, 0,
+
+       size_type nbdof = ps->linked_model().nb_dof();
+       darray x = in.pop().to_darray();
+       std::vector<double> xx(nbdof); gmm::copy(x, xx);
+       scalar_type gamma = in.pop().to_scalar();
+       darray t_x = in.pop().to_darray();
+       std::vector<double> tt_x(nbdof); gmm::copy(t_x, tt_x);
+       scalar_type t_gamma = in.pop().to_scalar();
+
+       ps->clear_sing_data();
+       getfem::treat_nonsmooth_point(*ps, xx, gamma, tt_x, t_gamma, 0);
+       );
+
+
     /address@hidden @CELL{X, gamma, T_X, T_gamma} = ('sing_data')
-      Return a singular point (`X`, `gamma`) encountered in the last
-      continuation step (if any) and a couple of arrays (`T_X`, `T_gamma`) of
-      tangents to all located solution branches, which emanate from 
address@hidden/
+      Return a singular point (`X`, `gamma`) stored in the @tcs object and a
+      couple of arrays (`T_X`, `T_gamma`) of tangents to all located solution
+      branches that emanate from address@hidden/
     sub_command
       ("sing_data", 0, 0, 0, 4,
+
        out.pop().from_dcvector(ps->get_x_sing());
        out.pop().from_scalar(ps->get_gamma_sing());
        out.pop().from_vector_container(ps->get_t_x_sing());

Modified: trunk/getfem/src/getfem/getfem_continuation.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_continuation.h?rev=4872&r1=4871&r2=4872&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_continuation.h       (original)
+++ trunk/getfem/src/getfem/getfem_continuation.h       Fri Mar  6 10:10:11 2015
@@ -1,7 +1,7 @@
 /* -*- c++ -*- (enables emacs c++ mode) */
 /*===========================================================================
  
- Copyright (C) 2011-2014 Tomas Ligursky, Yves Renard
+ Copyright (C) 2011-2015 Tomas Ligursky, Yves Renard
  
  This file is a part of GETFEM++
  
@@ -407,22 +407,24 @@
   }
   
   /* A tool for approximating a non-smooth point close to (x, gamma) and
-     locating (preferably) all smooth one-sided solution branches emanating
-     from there. It is supposed that (x, gamma) is a point on the previously
-     traversed smooth solution branch within the distance of S.h_min() from
-     the end point of this branch and (t_x, t_gamma) is the corresponding
-     tangent that is directed towards the end point. */
+     locating one-sided smooth solution branches emanating from there. It is
+     supposed that (x, gamma) is a point on a smooth solution branch within
+     the distance of S.h_min() from the end point of this branch and
+     (t_x, t_gamma) is the corresponding tangent that is directed towards the
+     end point. The value of version indicates whether the first new branch
+     found (if any) is to be chosen for further continuation (version = 1) or
+     not (version = 0). */
   template <typename CONT_S, typename VECT>
   void treat_nonsmooth_point(CONT_S &S, const VECT &x, double gamma,
                             const VECT &t_x, double t_gamma, int version) {
-    double gamma_end = gamma, Gamma, t_gamma0 = t_gamma, T_gamma = t_gamma,
+    double gamma_0 = gamma, Gamma, t_gamma0 = t_gamma, T_gamma = t_gamma,
       h = S.h_min(), cang, mcos = S.mincos();
-    VECT x_end(x), X(x), t_x0(t_x), T_x(t_x);
-
-    // approximate the non-smooth point by a bisection-like algorithm
+    VECT x_0(x), X(x), t_x0(t_x), T_x(t_x);
+
+    // approximate the end point more precisely by a bisection-like algorithm
     if (S.noisy() > 0)
       cout  << "starting locating a non-smooth point" << endl;
-    S.scaled_add(x, t_x, h, X); Gamma = gamma + h * t_gamma;
+    S.scaled_add(x_0, t_x0, h, X); Gamma = gamma_0 + h * t_gamma0;
     S.set_build(BUILD_ALL);
     if (newton_corr(S, X, Gamma, T_x, T_gamma, t_x0, t_gamma0)) {
       cang = S.cosang(T_x, t_x0, T_gamma, t_gamma0);
@@ -433,65 +435,66 @@
     h /= 2.;
     for (unsigned long i = 0; i < 15; i++) {
       if (S.noisy() > 0) cout << "prediction with h = " << h << endl;
-      S.scaled_add(x_end, t_x0, h, X); Gamma = gamma_end + h * t_gamma0;
+      S.scaled_add(x_0, t_x0, h, X); Gamma = gamma_0 + h * t_gamma0;
       S.set_build(BUILD_ALL);
       if (newton_corr(S, X, Gamma, T_x, T_gamma, t_x0, t_gamma0)
-         && (S.cosang(T_x, t_x, T_gamma, t_gamma) >= mcos)) {
-       S.copy(X, x_end); gamma_end = Gamma;
+         && (S.cosang(T_x, t_x0, T_gamma, t_gamma0) >= mcos)) {
+       S.copy(X, x_0); gamma_0 = Gamma;
        S.copy(T_x, t_x0); t_gamma0 = T_gamma;
       } else {
        S.copy(t_x0, T_x); T_gamma = t_gamma0;
       }
       h /= 2.;
     }
-    S.scaled_add(x_end, t_x0, h, x_end); gamma_end += h * t_gamma0;
-    S.set_sing_point(x_end, gamma_end);
-
-    // take two vectors to span a subspace of perturbations for the
-    // non-smooth point
+    S.set_sing_point(x_0, gamma_0);
+
+    // take two vectors to span a subspace of directions emanating from the
+    // end point
     if (S.noisy() > 0)
       cout  << "starting a thorough search for other branches" << endl;
     double t_gamma1 = t_gamma0, t_gamma2 = t_gamma0;
     VECT t_x1(t_x0), t_x2(t_x0);
     S.scale(t_x1, -1.); t_gamma1 *= -1.;
     S.insert_tangent_sing(t_x1, t_gamma1);
-
     h = S.h_min();
-    S.scaled_add(x_end, t_x0, h, X); Gamma = gamma_end + h * t_gamma0;
+    S.scaled_add(x_0, t_x0, h, X); Gamma = gamma_0 + h * t_gamma0;
     S.set_build(BUILD_ALL);
     compute_tangent(S, X, Gamma, t_x2, t_gamma2);
 
-    // perturb the non-smooth point systematically to find new tangent
-    // predictions
+    // try systematically various directions emanating from the end point for
+    // finding new possible tangent predictions
     unsigned long i1 = 0, i2 = 0, ncomb = 0;
     double a, a1, a2, no;
-    S.clear(t_x0); t_gamma0 = 0.;
 
     do {
       for (unsigned long i = 0; i < S.nbdir(); i++) {
        a = (2 * M_PI * double(i)) / double(S.nbdir());
-       a1 = h * sin(a); a2 = h * cos(a);
-       S.scaled_add(x_end, t_x1, a1, X); Gamma = gamma_end + a1 * t_gamma1;
-       S.scaled_add(X, t_x2, a2, X); Gamma += a2 * t_gamma2;
+       a1 = sin(a); a2 = cos(a);
+       S.scaled_add(t_x1, a1, t_x2, a2, T_x);
+       T_gamma = a1 * t_gamma1 + a2 * t_gamma2;
+       no = S.w_norm(T_x, T_gamma);
+       S.scaled_add(x_0, T_x, h / no, X);
+       Gamma = gamma_0 + h / no * T_gamma;
        S.set_build(BUILD_ALL);
        compute_tangent(S, X, Gamma, T_x, T_gamma);
 
-       if (S.abs(S.cosang(T_x, t_x0, T_gamma, t_gamma0)) < S.mincos()) {
+       if (S.abs(S.cosang(T_x, t_x0, T_gamma, t_gamma0)) < S.mincos()
+           || (i == 0 && ncomb == 0)) {
          S.copy(T_x, t_x0); t_gamma0 = T_gamma;
          if (S.insert_tangent_predict(T_x, T_gamma)) {
            if (S.noisy() > 0)
              cout << "new potential tangent vector found, "
                   << "trying one predictor-corrector step" << endl;
-           S.copy(x_end, X); Gamma = gamma_end;
-           
+           S.copy(x_0, X); Gamma = gamma_0;        
            if (test_predict_dir(S, X, Gamma, T_x, T_gamma)) {
              if (S.insert_tangent_sing(T_x, T_gamma)) {
-               if ((a == 0) && (ncomb == 0)
+               if ((i == 0) && (ncomb == 0) 
+                   // => (T_x, T_gamma) = (t_x2, t_gamma2)
                    && (S.abs(S.cosang(T_x, t_x0, T_gamma, t_gamma0))
-                       >= S.mincos())) { i2 = 1; ncomb = 1; }
+                       >= S.mincos())) { i2 = 1; }
                if (version) S.set_next_point(X, Gamma);
              }
-             S.copy(x_end, X); Gamma = gamma_end;
+             S.copy(x_0, X); Gamma = gamma_0;
              S.copy(t_x0, T_x); T_gamma = t_gamma0;
            }
            
@@ -504,12 +507,12 @@
       }
       
       // heuristics for varying the spanning vectors
-      bool index_changed;
-      if (i1 + 1 < i2) { ++i1; index_changed = true; }
+      bool perturb;
+      if (i1 + 1 < i2) { ++i1; perturb = false; }
       else if(i2 + 1 < S.nb_tangent_sing())
-       { ++i2; i1 = 0; index_changed = true; }
-      else index_changed = false;
-      if (index_changed) {
+       { ++i2; i1 = 0; perturb = false; }
+      else perturb = true;
+      if (!perturb) {
        S.copy(S.get_t_x_sing(i1), t_x1); t_gamma1 = S.get_t_gamma_sing(i1);
        S.copy(S.get_t_x_sing(i2), t_x2); t_gamma2 = S.get_t_gamma_sing(i2);
       } else {
@@ -517,7 +520,7 @@
        no = S.w_norm(T_x, T_gamma);
        S.scaled_add(t_x2, T_x, 0.1/no, t_x2);
        t_gamma2 += 0.1/no * T_gamma;
-       S.scaled_add(x_end, t_x2, h, X); Gamma = gamma_end + h * t_gamma2;
+       S.scaled_add(x_0, t_x2, h, X); Gamma = gamma_0 + h * t_gamma2;
        S.set_build(BUILD_ALL);
        compute_tangent(S, X, Gamma, t_x2, t_gamma2);
       }
@@ -628,40 +631,36 @@
     if (new_point) {
       S.copy(X, x); gamma = Gamma;
       S.copy(T_x, t_x); t_gamma = T_gamma;
-    } else if (S.non_smooth()) {
+    } else if (S.non_smooth() && S.singularities() > 1) {
       treat_nonsmooth_point(S, x, gamma, t_x0, t_gamma0, 1);
       if (S.next_point()) {
-       if (S.singularities() > 0) {
-         if (test_limit_point(S, T_gamma)) {
-           S.set_sing_label("limit point");
-           if (S.noisy() > 0) cout << "Limit point detected!" << endl;
-         }
-         if (S.singularities() > 1) {
-           if (S.noisy() > 0)
-             cout << "starting computing a test function for bifurcations"
-                  << endl;
-           S.set_build(BUILD_ALL);
-           bool bifurcation_detected = (S.nb_tangent_sing() > 2);
-           if (bifurcation_detected) {
-             // update the stored values of the test function only
-             S.set_tau_bp_1(tau_bp_init);
-             S.set_tau_bp_2(test_function_bp(S, S.get_x_next(),
-                                             S.get_gamma_next(),
-                                             S.get_t_x_sing(1),
-                                             S.get_t_gamma_sing(1)));
-           } else
-             bifurcation_detected
-               = test_nonsmooth_bifurcation(S, x, gamma, t_x, t_gamma,
-                                            S.get_x_next(),
-                                            S.get_gamma_next(),
-                                            S.get_t_x_sing(1),
-                                            S.get_t_gamma_sing(1));
-           if (bifurcation_detected) {
-             S.set_sing_label("non-smooth bifurcation point");
-             if (S.noisy() > 0)
-               cout << "Non-smooth bifurcation point detected!" << endl;
-           }
-         }
+       if (test_limit_point(S, T_gamma)) {
+         S.set_sing_label("limit point");
+         if (S.noisy() > 0) cout << "Limit point detected!" << endl;
+       }
+       if (S.noisy() > 0)
+         cout << "starting computing a test function for bifurcations"
+              << endl;
+       S.set_build(BUILD_ALL);
+       bool bifurcation_detected = (S.nb_tangent_sing() > 2);
+       if (bifurcation_detected) {
+         // update the stored values of the test function only
+         S.set_tau_bp_1(tau_bp_init);
+         S.set_tau_bp_2(test_function_bp(S, S.get_x_next(),
+                                         S.get_gamma_next(),
+                                         S.get_t_x_sing(1),
+                                         S.get_t_gamma_sing(1)));
+       } else
+         bifurcation_detected
+           = test_nonsmooth_bifurcation(S, x, gamma, t_x, t_gamma,
+                                        S.get_x_next(),
+                                        S.get_gamma_next(),
+                                        S.get_t_x_sing(1),
+                                        S.get_t_gamma_sing(1));
+       if (bifurcation_detected) {
+         S.set_sing_label("non-smooth bifurcation point");
+         if (S.noisy() > 0)
+           cout << "Non-smooth bifurcation point detected!" << endl;
        }
        
        S.copy(S.get_x_next(), x); gamma = S.get_gamma_next();