[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Markus Bürg |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Tue, 22 Aug 2017 06:19:04 -0400 (EDT) |
branch: mb-transInversion
commit 7e20e2c83bc977a3f4a31ff3972a0c8546ae30b7
Author: mb <address@hidden>
Date: Tue Aug 22 12:18:59 2017 +0200
Reimplement Newton method.
---
src/bgeot_geotrans_inv.cc | 91 ++++++++++++-----------------------------------
1 file changed, 22 insertions(+), 69 deletions(-)
diff --git a/src/bgeot_geotrans_inv.cc b/src/bgeot_geotrans_inv.cc
index 580bf43..eddcff9 100644
--- a/src/bgeot_geotrans_inv.cc
+++ b/src/bgeot_geotrans_inv.cc
@@ -298,80 +298,33 @@ vector<size_type> get_linear_nodes_indices(const string
&pgt_name, dim_type dim)
bool geotrans_inv_convex::invert_nonlin(const base_node& xreal,
base_node& x, scalar_type IN_EPS,
bool &converged, bool throw_except) {
- using namespace gmm;
-
converged = true;
- base_node xn(P), y, z,x0;
/* find an initial guess */
nonlinear_storage.plinear_inversion->invert(xreal, x);
project_into_convex(x, *nonlinear_storage.plinear_inversion->pgt);
-
- base_node vres(P);
- base_node rn(xreal); rn -= y;
- pgt->poly_vector_grad(x, pc);
- update_B();
-
- gmm::mult(gmm::transposed(K), rn, vres);
- scalar_type res = gmm::vect_norm2(vres);
-
- //cerr << "DEBUT: res0=" << res << ", X=" << xreal << "\nB=" << B << ",
K=" << K << "\n" << ", pc=" << pc << "\n";
- unsigned cnt = 50;
- while (res > EPS/10 && --cnt) {
- gmm::mult(gmm::transposed(B), rn, xn);
- scalar_type newres;
- for (unsigned i=1; i<=256; i*=2) {
- z = x + xn / scalar_type(i);
- y = pgt->transform(z, G);
-
- rn = xreal - y;
-
- pgt->poly_vector_grad(z, pc);
- update_B();
-
- // cout << "K = " << K << endl;
-
- if (P != N) {
- gmm::mult(gmm::transposed(K), rn, vres);
- newres = gmm::vect_norm2(vres);
- } else {
- newres = gmm::vect_norm2(rn); // "better" residu
- }
- if (newres < 1.5*res) break;
- }
- x = z; res = newres;
- // cerr << "cnt=" << cnt << ", x=" << x << ", res=" << res << "\n";
+ nonlinear_storage.x_real = pgt->transform(x, G);
+ add(nonlinear_storage.x_real, scaled(xreal, -1.0), nonlinear_storage.diff);
+
+ auto res = vect_norm2(nonlinear_storage.diff);
+ double res0 = 1e10;
+ auto cnt = 0;
+
+ while (res > IN_EPS) {
+ if (abs(res - res0) < IN_EPS) return false;
+
+ pgt->poly_vector_grad(x, pc);
+ update_B();
+ mult(transposed(B), nonlinear_storage.diff, nonlinear_storage.x_ref);
+ add(scaled(nonlinear_storage.x_ref, -1.0), x);
+ project_into_convex(x, *nonlinear_storage.plinear_inversion->pgt);
+ nonlinear_storage.x_real = pgt->transform(x, G);
+ add(nonlinear_storage.x_real, scaled(xreal, -1.0),
nonlinear_storage.diff);
+ res0 = res;
+ res = vect_norm2(nonlinear_storage.diff);
+ ++cnt;
}
- //cerr << " invert_nonlin done\n";
- //cerr << "cnt=" << cnt << ", P=" << P << ", N=" << N << ", G=" << G <<
"\nX=" << xreal << " Xref=" << x << "\nresidu=" << res << "\nB=" << B << ", K="
<< K << "\n" << ", pc=" << pc << "\n-------------------^^^^^^^^\n";
- if (cnt == 0) {
- //cout << "BFGS in geotrans_inv_convex!\n";
- geotrans_inv_convex_bfgs b(*this, xreal);
- gmm::iteration iter(EPS,0);
- x = x0;
- gmm::bfgs(b,b,x,10,iter);
- rn = pgt->transform(x, G) - xreal;
-
- if (pgt->convex_ref()->is_in(x) < IN_EPS &&
- (N==P && gmm::vect_norm2(rn) > IN_EPS)) {
- GMM_ASSERT1(!throw_except,
- "inversion of non-linear geometric transformation "
- "failed ! (too much iterations -- xreal=" << xreal
- << ", rn=" << rn
- << ", xref=" << x
- << ", is_in(x)=" << pgt->convex_ref()->is_in(x)
- << ", eps=" << IN_EPS << ")");
- converged = false;
- return false;
- }
- }
- // Test un peu sev�re peut-�tre en ce qui concerne rn.
- if (pgt->convex_ref()->is_in(x) < IN_EPS
- && (P == N || gmm::vect_norm2(rn) < IN_EPS)) {
- // cout << "point " << x << "is IN (" << pgt->convex_ref()->is_in(x)
- // << ")\n";
- return true;
- } // else cout << "point IS OUT\n";
- return false;
+
+ return true;
}
} /* end of namespace bgeot. */