branch: debug_assembly_breakage commit 70beae92f772491a27ee8fdb5b1d9cdb7dafa94f Author: Andriy.Andreykiv <a...@plaxis.com> Date: Thu Mar 14 13:29:31 2019 +0100
re-introducing the same changes as in the latest commits in the master, to be able to compile --- src/getfem/bgeot_rtree.h | 44 +++++++++++- src/getfem/getfem_model_solvers.h | 115 +++++++++++++++++++++++++++++++ src/getfem_generic_assembly_semantic.cc | 10 +-- src/getfem_model_solvers.cc | 116 -------------------------------- 4 files changed, 161 insertions(+), 124 deletions(-) diff --git a/src/getfem/bgeot_rtree.h b/src/getfem/bgeot_rtree.h index 397eec3..892d3b6 100644 --- a/src/getfem/bgeot_rtree.h +++ b/src/getfem/bgeot_rtree.h @@ -47,13 +47,51 @@ namespace bgeot { size_type id; base_node min, max; }; - + + /** Wraps "const box_index *" but overloads + * comparison operators based on id and not + * addresses. This ensures deterministic ordering in sets. + */ + struct box_index_ptr { + box_index_ptr(const box_index *p) + : p_{p} + {} + + box_index_ptr(const box_index_ptr&) = default; + + inline bool operator < (const box_index_ptr &bptr) const { + return p_->id < bptr.p_->id; + } + + inline bool operator == (const box_index_ptr &bptr) const { + return p_->id == bptr.p_->id; + } + + inline bool operator != (const box_index_ptr &bptr) const { + return !(*this == bptr); + } + + inline operator const box_index *() const { + return p_; + } + + inline const box_index * operator->() const { + return p_; + } + + inline const box_index& operator*() const { + return *p_; + } + + const box_index *p_; + }; + struct rtree_elt_base { enum { RECTS_PER_LEAF=8 }; bool isleaf_; bool isleaf() const { return isleaf_; } base_node rmin, rmax; - rtree_elt_base(bool leaf, const base_node& rmin_, const base_node& rmax_) + rtree_elt_base(bool leaf, const base_node& rmin_, const base_node& rmax_) : isleaf_(leaf), rmin(rmin_), rmax(rmax_) {} virtual ~rtree_elt_base() {} }; @@ -67,7 +105,7 @@ namespace bgeot { public: typedef std::deque<box_index> box_cont; typedef std::vector<const box_index*> pbox_cont; - typedef std::set<const box_index*> pbox_set; + typedef std::set<box_index_ptr> pbox_set; void add_box(base_node min, base_node max, size_type id=size_type(-1)) { box_index bi; bi.min = min; bi.max = max; diff --git a/src/getfem/getfem_model_solvers.h b/src/getfem/getfem_model_solvers.h index 6f1d8d4..6280bd2 100644 --- a/src/getfem/getfem_model_solvers.h +++ b/src/getfem/getfem_model_solvers.h @@ -719,6 +719,121 @@ namespace getfem { return select_linear_solver<model_complex_sparse_matrix, model_complex_plain_vector>(md, name); } + + /* ***************************************************************** */ + /* 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. diff --git a/src/getfem_generic_assembly_semantic.cc b/src/getfem_generic_assembly_semantic.cc index 905dc62..fbe8832 100644 --- a/src/getfem_generic_assembly_semantic.cc +++ b/src/getfem_generic_assembly_semantic.cc @@ -330,7 +330,7 @@ namespace getfem { break; case GA_NODE_INTERPOLATE_DERIVATIVE: c += 2.321*ga_hash_code(pnode->interpolate_name_der); - [[fallthrough]]; // The hash code is completed with the next item + //[[fallthrough]]; // The hash code is completed with the next item case GA_NODE_INTERPOLATE_VAL: case GA_NODE_INTERPOLATE_GRAD: case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG: case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST: @@ -525,7 +525,7 @@ namespace getfem { if (tree.secondary_domain.size() == 0) ga_throw_error(pnode->expr, pnode->pos, "Secondary domain used " "in a single domain term."); - [[fallthrough]]; + //[[fallthrough]]; case GA_NODE_INTERPOLATE: if (pnode->name.compare("X") == 0) { if (pnode->node_type == GA_NODE_INTERPOLATE) { @@ -549,7 +549,7 @@ namespace getfem { } break; } - [[fallthrough]]; + //[[fallthrough]]; case GA_NODE_ELEMENTARY: // and ... case GA_NODE_INTERPOLATE: case GA_NODE_XFEM_PLUS: case GA_NODE_XFEM_MINUS: @@ -2832,7 +2832,7 @@ namespace getfem { case GA_NODE_INTERPOLATE_FILTER: if (!child_0_is_constant) { is_constant = false; break; } - [[fallthrough]]; + //[[fallthrough]]; case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST: @@ -3807,7 +3807,7 @@ namespace getfem { if ((interpolate_node || interpolate_test_node) && (workspace.associated_mf(pnode->name) != 0)) marked = true; - + if (pnode->node_type == GA_NODE_INTERPOLATE_X || pnode->node_type == GA_NODE_INTERPOLATE_NORMAL) marked = true; 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. */ /* ***************************************************************** */