[Getfem-commits] (no subject)

[Konstantinos Poulios]

branch: master
commit 65084cef321ce1d4824144042c77ac56815cf348
Author: Konstantinos Poulios <address@hidden>
Date:   Wed Jan 30 05:59:16 2019 +0100

    Move model_pb class from header to implementation file and remove 
scale_residual() function
---
 src/getfem/getfem_model_solvers.h | 125 +-------------------------------------
 src/getfem_model_solvers.cc       | 115 +++++++++++++++++++++++++++++++++++
 2 files changed, 117 insertions(+), 123 deletions(-)

diff --git a/src/getfem/getfem_model_solvers.h 
b/src/getfem/getfem_model_solvers.h
index 6aa3b5e..6f1d8d4 100644
--- a/src/getfem/getfem_model_solvers.h
+++ b/src/getfem/getfem_model_solvers.h
@@ -454,7 +454,7 @@ namespace getfem {
         while (is_singular) { // Linear system solve
           pb.compute_tangent_matrix();
           gmm::clear(dr);
-          gmm::copy(gmm::scaled(pb.residual(), pb.scale_residual()), b);
+          gmm::copy(pb.residual(), b);
           gmm::add(gmm::scaled(b0,alpha-R(1)), b);
           if (iter.get_noisy() > 1) cout << "starting linear solver" << endl;
           iter_linsolv.init();
@@ -575,7 +575,7 @@ namespace getfem {
       while (is_singular) {
         pb.compute_tangent_matrix();
         gmm::clear(dr);
-        gmm::copy(gmm::scaled(pb.residual(), pb.scale_residual()), b);
+        gmm::copy(pb.residual(), b);
         if (iter.get_noisy() > 1) cout << "starting linear solver" << endl;
         iter_linsolv.init();
         linear_solver(pb.tangent_matrix(), dr, b, iter_linsolv);
@@ -608,127 +608,6 @@ namespace getfem {
   }
 
 
-  /* ***************************************************************** */
-  /*  Intermediary structure for Newton algorithms with getfem::model. */
-  /* ***************************************************************** */
-
-  #define TRACE_SOL 0
-
-  template <typename MAT, typename VEC>
-  struct model_pb {
-
-    typedef MAT MATRIX;
-    typedef VEC VECTOR;
-    typedef typename gmm::linalg_traits<VECTOR>::value_type T;
-    typedef typename gmm::number_traits<T>::magnitude_type R;
-
-    model &md;
-    abstract_newton_line_search &ls;
-    VECTOR stateinit, &state;
-    const VECTOR &rhs;
-    const MATRIX &K;
-
-    void compute_tangent_matrix(void) {
-      md.to_variables(state);
-      md.assembly(model::BUILD_MATRIX);
-    }
-
-    const MATRIX &tangent_matrix(void) { return K; }
-
-    inline T scale_residual(void) const { return T(1); }
-
-    void compute_residual(void) {
-      md.to_variables(state);
-      md.assembly(model::BUILD_RHS);
-    }
-
-    void perturbation(void) {
-      R res = gmm::vect_norm2(state), ampl = std::max(res*R(1E-20), R(1E-50));
-      std::vector<R> V(gmm::vect_size(state));
-      gmm::fill_random(V);
-      gmm::add(gmm::scaled(V, ampl), state);
-    }
-
-    const VECTOR &residual(void) const { return rhs; }
-    const VECTOR &state_vector(void) const { return state; }
-    VECTOR &state_vector(void) { return state; }
-
-    R state_norm(void) const
-    { return gmm::vect_norm1(state); }
-
-    R approx_external_load_norm(void)
-    { return md.approx_external_load(); }
-
-    R residual_norm(void) { // A norm1 seems to be better than a norm2
-                            // at least for contact problems.
-      return gmm::vect_norm1(rhs);
-    }
-
-    R compute_res(bool comp = true) {
-      if (comp) compute_residual();
-      return residual_norm();
-    }
-
-
-    R line_search(VECTOR &dr, const gmm::iteration &iter) {
-      size_type nit = 0;
-      gmm::resize(stateinit, md.nb_dof());
-      gmm::copy(state, stateinit);
-      R alpha(1), res, /* res_init, */ R0;
-
-      /* res_init = */ res = compute_res(false);
-      // cout << "first residual = " << residual() << endl << endl;
-      R0 = gmm::real(gmm::vect_sp(dr, rhs));
-
-#if TRACE_SOL
-      static int trace_number = 0;
-      int trace_iter = 0;
-      {
-        std::stringstream trace_name;
-        trace_name << "line_search_state" << std::setfill('0')
-                   << std::setw(3) << trace_number << "_000_init";
-        gmm::vecsave(trace_name.str(),stateinit);
-      }
-      trace_number++;
-#endif
-
-      ls.init_search(res, iter.get_iteration(), R0);
-      do {
-        alpha = ls.next_try();
-        gmm::add(stateinit, gmm::scaled(dr, alpha), state);
-#if TRACE_SOL
-        {
-          trace_iter++;
-          std::stringstream trace_name;
-          trace_name  << "line_search_state" << std::setfill('0')
-                      << std::setw(3) << trace_number << "_"
-                      << std::setfill('0') << std::setw(3) << trace_iter;
-          gmm::vecsave(trace_name.str(), state);
-        }
-#endif
-        res = compute_res();
-        // cout << "residual = " << residual() << endl << endl;
-        R0 = gmm::real(gmm::vect_sp(dr, rhs));
-
-        ++ nit;
-      } while (!ls.is_converged(res, R0));
-
-      if (alpha != ls.converged_value()) {
-        alpha = ls.converged_value();
-        gmm::add(stateinit, gmm::scaled(dr, alpha), state);
-        res = ls.converged_residual();
-        compute_residual();
-      }
-
-      return alpha;
-    }
-
-    model_pb(model &m, abstract_newton_line_search &ls_, VECTOR &st,
-             const VECTOR &rhs_, const MATRIX &K_)
-      : md(m), ls(ls_), state(st), rhs(rhs_), K(K_) {}
-
-  };
-
   //---------------------------------------------------------------------
   // Default linear solver.
   //---------------------------------------------------------------------
diff --git a/src/getfem_model_solvers.cc b/src/getfem_model_solvers.cc
index 2a45398..9e70e14 100644
--- a/src/getfem_model_solvers.cc
+++ b/src/getfem_model_solvers.cc
@@ -119,6 +119,121 @@ namespace getfem {
 
 
   /* ***************************************************************** */
+  /*  Intermediary structure for Newton algorithms with getfem::model. */
+  /* ***************************************************************** */
+
+  #define TRACE_SOL 0
+
+  template <typename MAT, typename VEC>
+  struct model_pb {
+
+    typedef MAT MATRIX;
+    typedef VEC VECTOR;
+    typedef typename gmm::linalg_traits<VECTOR>::value_type T;
+    typedef typename gmm::number_traits<T>::magnitude_type R;
+
+    model &md;
+    abstract_newton_line_search &ls;
+    VECTOR stateinit, &state;
+    const VECTOR &rhs;
+    const MATRIX &K;
+
+    void compute_tangent_matrix() {
+      md.to_variables(state);
+      md.assembly(model::BUILD_MATRIX);
+    }
+
+    const MATRIX &tangent_matrix() { return K; }
+
+    void compute_residual() {
+      md.to_variables(state);
+      md.assembly(model::BUILD_RHS);
+    }
+
+    void perturbation() {
+      R res = gmm::vect_norm2(state), ampl = std::max(res*R(1E-20), R(1E-50));
+      std::vector<R> V(gmm::vect_size(state));
+      gmm::fill_random(V);
+      gmm::add(gmm::scaled(V, ampl), state);
+    }
+
+    const VECTOR &residual() const { return rhs; }
+    const VECTOR &state_vector() const { return state; }
+    VECTOR &state_vector() { return state; }
+
+    R state_norm() const { return gmm::vect_norm1(state); }
+
+    R approx_external_load_norm() { return md.approx_external_load(); }
+
+    R residual_norm() {
+      return gmm::vect_norm1(rhs); // A norm1 seems to be better than a norm2
+    }                              // at least for contact problems.
+
+    R compute_res(bool comp = true) {
+      if (comp) compute_residual();
+      return residual_norm();
+    }
+
+    R line_search(VECTOR &dr, const gmm::iteration &iter) {
+      size_type nit = 0;
+      gmm::resize(stateinit, md.nb_dof());
+      gmm::copy(state, stateinit);
+      R alpha(1), res, /* res_init, */ R0;
+
+      /* res_init = */ res = compute_res(false);
+      // cout << "first residual = " << residual() << endl << endl;
+      R0 = gmm::real(gmm::vect_sp(dr, rhs));
+
+#if TRACE_SOL
+      static int trace_number = 0;
+      int trace_iter = 0;
+      {
+        std::stringstream trace_name;
+        trace_name << "line_search_state" << std::setfill('0')
+                   << std::setw(3) << trace_number << "_000_init";
+        gmm::vecsave(trace_name.str(),stateinit);
+      }
+      trace_number++;
+#endif
+
+      ls.init_search(res, iter.get_iteration(), R0);
+      do {
+        alpha = ls.next_try();
+        gmm::add(stateinit, gmm::scaled(dr, alpha), state);
+#if TRACE_SOL
+        {
+          trace_iter++;
+          std::stringstream trace_name;
+          trace_name  << "line_search_state" << std::setfill('0')
+                      << std::setw(3) << trace_number << "_"
+                      << std::setfill('0') << std::setw(3) << trace_iter;
+          gmm::vecsave(trace_name.str(), state);
+        }
+#endif
+        res = compute_res();
+        // cout << "residual = " << residual() << endl << endl;
+        R0 = gmm::real(gmm::vect_sp(dr, rhs));
+
+        ++ nit;
+      } while (!ls.is_converged(res, R0));
+
+      if (alpha != ls.converged_value()) {
+        alpha = ls.converged_value();
+        gmm::add(stateinit, gmm::scaled(dr, alpha), state);
+        res = ls.converged_residual();
+        compute_residual();
+      }
+
+      return alpha;
+    }
+
+    model_pb(model &m, abstract_newton_line_search &ls_, VECTOR &st,
+             const VECTOR &rhs_, const MATRIX &K_)
+      : md(m), ls(ls_), state(st), rhs(rhs_), K(K_) {}
+
+  };
+
+  /* ***************************************************************** */
   /*     Standard solve.                                               */
   /* ***************************************************************** */