[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r4872 - in /trunk/getfem: interface/src/gf_cont_struct.
From: |
tomas . ligursky |
Subject: |
[Getfem-commits] r4872 - in /trunk/getfem: interface/src/gf_cont_struct.cc interface/src/gf_cont_struct_get.cc src/getfem/getfem_continuation.h |
Date: |
Fri, 06 Mar 2015 09:10:12 -0000 |
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();
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [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 <=