[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Andriy Andreykiv |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Thu, 14 Mar 2019 08:29:45 -0400 (EDT) |
branch: debug_assembly_breakage
commit 70beae92f772491a27ee8fdb5b1d9cdb7dafa94f
Author: Andriy.Andreykiv <address@hidden>
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. */
/* ***************************************************************** */