branch: consistent_rtree_box_order commit 00ddea3d79633e1a519535d11eaf84bc75d5ec42 Author: Andriy.Andreykiv <a...@plaxis.com> 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. */ /* ***************************************************************** */