[Getfem-commits] (no subject)

[Andriy Andreykiv]

branch: consistent_rtree_box_order
commit 00ddea3d79633e1a519535d11eaf84bc75d5ec42
Author: Andriy.Andreykiv <address@hidden>
Date:   Fri Mar 8 15:03:32 2019 +0100

    returning model_pb to the header as it's being used elsewhere
---
 src/getfem/getfem_model_solvers.h | 115 +++++++++++++++++++++++++++++++++++++
 src/getfem_model_solvers.cc       | 116 --------------------------------------
 2 files changed, 115 insertions(+), 116 deletions(-)

diff --git a/src/getfem/getfem_model_solvers.h 
b/src/getfem/getfem_model_solvers.h
index 6f1d8d4..3ac5381 100644
--- a/src/getfem/getfem_model_solvers.h
+++ b/src/getfem/getfem_model_solvers.h
@@ -608,6 +608,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_) {}
+
+  };
+
   //---------------------------------------------------------------------
   // Default linear solver.
   //---------------------------------------------------------------------
diff --git a/src/getfem_model_solvers.cc b/src/getfem_model_solvers.cc
index 9e70e14..b43c0ec 100644
--- a/src/getfem_model_solvers.cc
+++ b/src/getfem_model_solvers.cc
@@ -117,122 +117,6 @@ namespace getfem {
     md.set_time_integration(1);
   }
 
-
-  /* ***************************************************************** */
-  /*  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.                                               */
   /* ***************************************************************** */