[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: |
Wed, 13 Mar 2019 06:25:10 -0400 (EDT) |
branch: before_assembly_failure
commit 05030d5457b573a31257fb7f1fad89b72277f2da
Author: Andriy.Andreykiv <address@hidden>
Date: Wed Mar 13 11:24:48 2019 +0100
everting 2f7396a46dadad266ba882623d300d1a19e274dc and
ccb9a29bed98470d88ef896976b38e6eab147b07
that break assembly and re-introducing again the commits that follow
---
src/getfem/bgeot_rtree.h | 44 +++++++++++-
src/getfem/bgeot_tensor.h | 17 +++--
src/getfem/getfem_model_solvers.h | 115 +++++++++++++++++++++++++++++++
src/getfem_generic_assembly_semantic.cc | 8 +--
src/getfem_model_solvers.cc | 116 --------------------------------
5 files changed, 168 insertions(+), 132 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/bgeot_tensor.h b/src/getfem/bgeot_tensor.h
index 815f6eb..2629a79 100644
--- a/src/getfem/bgeot_tensor.h
+++ b/src/getfem/bgeot_tensor.h
@@ -222,14 +222,14 @@ namespace bgeot {
}
inline void init(size_type i, size_type j, size_type k) {
- sizes_.resize(3); sizes_[0] = i; sizes_[1] = j; sizes_[2] = k;
+ sizes_.resize(3); sizes_[0] = i; sizes_[1] = j; sizes_[2] = k;
coeff.resize(3); coeff[0] = 1; coeff[1] = i; coeff[2] = i*j;
this->resize(i*j*k);
}
inline void init(size_type i, size_type j, size_type k, size_type l) {
sizes_.resize(4);
- sizes_[0] = i; sizes_[1] = j; sizes_[2] = k; sizes_[3] = k;
+ sizes_[0] = i; sizes_[1] = j; sizes_[2] = k; sizes_[3] = k;
coeff.resize(4);
coeff[0] = 1; coeff[1] = i; coeff[2] = i*j; coeff[3] = i*j*k;
this->resize(i*j*k*l);
@@ -243,7 +243,7 @@ namespace bgeot {
{ init(i, j, k); }
inline void adjust_sizes(size_type i, size_type j, size_type k, size_type
l)
{ init(i, j, k, l); }
-
+
inline size_type adjust_sizes_changing_last(const tensor &t, size_type P) {
const multi_index &mi = t.sizes_; size_type d = mi.size();
sizes_.resize(d); coeff.resize(d);
@@ -323,12 +323,11 @@ namespace bgeot {
}
tensor(const multi_index &c) { init(c); }
- tensor(size_type i, size_type j)
- { init(multi_index(i, j)); }
- tensor(size_type i, size_type j, size_type k)
- { init(multi_index(i, j, k)); }
+ tensor(size_type i) = delete; // { init(i); }
+ tensor(size_type i, size_type j) { init(i, j); }
+ tensor(size_type i, size_type j, size_type k) { init(i, j, k); }
tensor(size_type i, size_type j, size_type k, size_type l)
- { init(multi_index(i, j, k, l)); }
+ { init(i, j, k, l); }
tensor() {}
};
@@ -406,7 +405,7 @@ namespace bgeot {
/* reduction du tenseur t par son indice ni et la matrice m. */
DEFINE_STATIC_THREAD_LOCAL(std::vector<T>, tmp);
DEFINE_STATIC_THREAD_LOCAL(multi_index, mi);
-
+
mi = t.sizes();
size_type dimt = mi[ni], dim = m.ncols();
GMM_ASSERT2(dimt, "Inconsistent dimension.");
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..09cc7e5 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:
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. */
/* ***************************************************************** */