[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Wed, 6 Sep 2017 14:07:17 -0400 (EDT) |
branch: mb-transInversion
commit 06b29a6bac25b5ca7b52a73fcfaf491c591c636a
Author: Yves Renard <address@hidden>
Date: Wed Sep 6 20:06:59 2017 +0200
modifications for compilation problems on g++ and cosmetic changes
---
src/bgeot_geotrans_inv.cc | 244 ++++++++++++++++++++--------------------
src/getfem/bgeot_geotrans_inv.h | 58 +++++-----
2 files changed, 150 insertions(+), 152 deletions(-)
diff --git a/src/bgeot_geotrans_inv.cc b/src/bgeot_geotrans_inv.cc
index cdbe996..5d2744a 100644
--- a/src/bgeot_geotrans_inv.cc
+++ b/src/bgeot_geotrans_inv.cc
@@ -52,17 +52,17 @@ namespace bgeot
} else return invert_nonlin(n, n_ref,IN_EPS,converged, false);
}
-bool point_in_convex(
- const geometric_trans &geoTrans,
- const base_node &x,
- scalar_type res,
- scalar_type IN_EPS) {
- // Test un peu sev�re peut-�tre en ce qui concerne res.
- return (geoTrans.convex_ref()->is_in(x) < IN_EPS) && (res < IN_EPS);
-}
+ bool point_in_convex(const geometric_trans &geoTrans,
+ const base_node &x,
+ scalar_type res,
+ scalar_type IN_EPS) {
+ // Test un peu sev�re peut-�tre en ce qui concerne res.
+ return (geoTrans.convex_ref()->is_in(x) < IN_EPS) && (res < IN_EPS);
+ }
/* inversion for linear geometric transformations */
- bool geotrans_inv_convex::invert_lin(const base_node& n, base_node& n_ref,
scalar_type IN_EPS) {
+ bool geotrans_inv_convex::invert_lin(const base_node& n, base_node& n_ref,
+ scalar_type IN_EPS) {
base_node y(n); for (size_type i=0; i < N; ++i) y[i] -= G(i,0);
mult(transposed(B), y, n_ref);
y = pgt->transform(n_ref, G);
@@ -108,132 +108,130 @@ bool point_in_convex(
}
};
-nonlinear_storage::linearised_structure::linearised_structure(
- const convex_ind_ct &direct_points_indices,
- const stored_point_tab &reference_nodes,
- const std::vector<base_node> &real_nodes) :
- diff(real_nodes.begin()->size()), diff_ref(direct_points_indices.size() -
1)
-{
- auto n_points = direct_points_indices.size();
- std::vector<base_node> direct_points;
- std::vector<base_node> direct_points_ref;
- direct_points.reserve(n_points);
- direct_points_ref.reserve(n_points);
-
- for (auto i : direct_points_indices) {
- direct_points.push_back(real_nodes[i]);
- direct_points_ref.push_back(reference_nodes[i]);
- }
-
- auto N = direct_points.begin()->size();
- auto N_ref = direct_points_ref.begin()->size();
- base_matrix K_linear(N, n_points - 1);
- B_linear.base_resize(N, n_points - 1);
- K_ref_linear.base_resize(N_ref, n_points - 1);
- P_linear = direct_points[0];
- P_ref_linear = direct_points_ref[0];
-
- for (size_type i = 1; i < n_points; ++i) {
- add(direct_points[i], scaled(P_linear, -1.0), mat_col(K_linear, i - 1));
- add(direct_points_ref[i], scaled(P_ref_linear, -1.0),
mat_col(K_ref_linear, i - 1));
- }
-
- if (K_linear.nrows() == K_linear.ncols()) {
- lu_inverse(K_linear);
- copy(transposed(K_linear), B_linear);
- }
- else {
- base_matrix temp(n_points - 1, n_points - 1);
- mult(transposed(K_linear), K_linear, temp);
- lu_inverse(temp);
- mult(K_linear, temp, B_linear);
- }
-}
-
-void nonlinear_storage::linearised_structure::invert(
- const base_node &x_real, base_node &x_ref, scalar_type IN_EPS) const
-{
- add(x_real, scaled(P_linear, -1.0), diff);
- mult(transposed(B_linear), diff, diff_ref);
- mult(K_ref_linear, diff_ref, x_ref);
- add(P_ref_linear, x_ref);
-}
-
-void project_into_convex(base_node &x, const geometric_trans *pgeo_trans, bool
project) {
- if (project == false) return;
-
- auto dim = pgeo_trans->dim();
-
- for (auto d = 0; d < dim; ++d) {
- if (x[d] < 0.0) x[d] = 0.0;
- if (x[d] > 1.0) x[d] = 1.0;
- }
-
- auto poriginal_trans = pgeo_trans;
-
- if (auto ptorus_trans = dynamic_cast<const torus_geom_trans*>(pgeo_trans)) {
- poriginal_trans = ptorus_trans->get_original_transformation().get();
+ nonlinear_storage_struct::linearised_structure::linearised_structure
+ (const convex_ind_ct &direct_points_indices,
+ const stored_point_tab &reference_nodes,
+ const std::vector<base_node> &real_nodes) :
+ diff(real_nodes.begin()->size()), diff_ref(direct_points_indices.size()-1)
{
+ auto n_points = direct_points_indices.size();
+ std::vector<base_node> direct_points;
+ std::vector<base_node> direct_points_ref;
+ direct_points.reserve(n_points);
+ direct_points_ref.reserve(n_points);
+
+ for (auto i : direct_points_indices) {
+ direct_points.push_back(real_nodes[i]);
+ direct_points_ref.push_back(reference_nodes[i]);
+ }
+
+ auto N = direct_points.begin()->size();
+ auto N_ref = direct_points_ref.begin()->size();
+ base_matrix K_linear(N, n_points - 1);
+ B_linear.base_resize(N, n_points - 1);
+ K_ref_linear.base_resize(N_ref, n_points - 1);
+ P_linear = direct_points[0];
+ P_ref_linear = direct_points_ref[0];
+
+ for (size_type i = 1; i < n_points; ++i) {
+ add(direct_points[i], scaled(P_linear, -1.0), mat_col(K_linear, i - 1));
+ add(direct_points_ref[i], scaled(P_ref_linear, -1.0),
+ mat_col(K_ref_linear, i - 1));
+ }
+
+ if (K_linear.nrows() == K_linear.ncols()) {
+ lu_inverse(K_linear);
+ copy(transposed(K_linear), B_linear);
+ }
+ else {
+ base_matrix temp(n_points - 1, n_points - 1);
+ mult(transposed(K_linear), K_linear, temp);
+ lu_inverse(temp);
+ mult(K_linear, temp, B_linear);
+ }
}
- auto pbasic_convex_ref = basic_convex_ref(poriginal_trans->convex_ref());
- auto nb_simplices = pbasic_convex_ref->simplexified_convex()->nb_convex();
-
- if (nb_simplices == 1) { //simplex
- auto sum_coordinates = 0.0;
-
- for (auto d = 0; d < dim; ++d) sum_coordinates += x[d];
-
- if (sum_coordinates > 1.0) scale(x, 1.0 / sum_coordinates);
+ void nonlinear_storage_struct::linearised_structure::invert
+ (const base_node &x_real, base_node &x_ref, scalar_type /* IN_EPS */) const {
+ add(x_real, scaled(P_linear, -1.0), diff);
+ mult(transposed(B_linear), diff, diff_ref);
+ mult(K_ref_linear, diff_ref, x_ref);
+ add(P_ref_linear, x_ref);
}
- else if ((dim == 3) && (nb_simplices != 4)) { //prism
- auto sum_coordinates = x[0] + x[1];
- if (sum_coordinates > 1.0) {
- x[0] /= sum_coordinates;
- x[1] /= sum_coordinates;
+ void project_into_convex(base_node &x, const geometric_trans *pgeo_trans,
+ bool project) {
+ if (project == false) return;
+
+ auto dim = pgeo_trans->dim();
+
+ for (auto d = 0; d < dim; ++d) {
+ if (x[d] < 0.0) x[d] = 0.0;
+ if (x[d] > 1.0) x[d] = 1.0;
}
- }
-}
-void find_initial_guess(
- base_node &x,
- nonlinear_storage_struct &storage,
- const base_node &xreal,
- const base_matrix &G,
- const geometric_trans *pgt,
- scalar_type IN_EPS)
-{
- storage.x_ref = pgt->geometric_nodes()[0];
-
- auto res = vect_dist2(mat_col(G, 0), xreal);
- double res0;
-
- for (size_type j = 1; j < pgt->nb_points(); ++j) {
- res0 = vect_dist2(mat_col(G, j), xreal);
- if (res > res0) {
- res = res0;
- storage.x_ref = pgt->geometric_nodes()[j];
+ auto poriginal_trans = pgeo_trans;
+
+ if (auto ptorus_trans = dynamic_cast<const torus_geom_trans*>(pgeo_trans))
{
+ poriginal_trans = ptorus_trans->get_original_transformation().get();
+ }
+
+ auto pbasic_convex_ref = basic_convex_ref(poriginal_trans->convex_ref());
+ auto nb_simplices = pbasic_convex_ref->simplexified_convex()->nb_convex();
+
+ if (nb_simplices == 1) { // simplex
+ auto sum_coordinates = 0.0;
+
+ for (auto d = 0; d < dim; ++d) sum_coordinates += x[d];
+
+ if (sum_coordinates > 1.0) scale(x, 1.0 / sum_coordinates);
+ }
+ else if ((dim == 3) && (nb_simplices != 4)) { // prism
+ auto sum_coordinates = x[0] + x[1];
+
+ if (sum_coordinates > 1.0) {
+ x[0] /= sum_coordinates;
+ x[1] /= sum_coordinates;
+ }
}
}
-
- res0 = std::numeric_limits<scalar_type>::max();
-
- if (storage.plinearised_structure != nullptr) {
- storage.plinearised_structure->invert(xreal, x, IN_EPS);
- project_into_convex(x, pgt, storage.project_into_element);
- res0 = vect_dist2(pgt->transform(x, G), xreal);
+
+ void find_initial_guess(base_node &x,
+ nonlinear_storage_struct &storage,
+ const base_node &xreal,
+ const base_matrix &G,
+ const geometric_trans *pgt,
+ scalar_type IN_EPS) {
+ storage.x_ref = pgt->geometric_nodes()[0];
+
+ auto res = vect_dist2(mat_col(G, 0), xreal);
+ double res0;
+
+ for (size_type j = 1; j < pgt->nb_points(); ++j) {
+ res0 = vect_dist2(mat_col(G, j), xreal);
+ if (res > res0) {
+ res = res0;
+ storage.x_ref = pgt->geometric_nodes()[j];
+ }
+ }
+
+ res0 = std::numeric_limits<scalar_type>::max();
+
+ if (storage.plinearised_structure != nullptr) {
+ storage.plinearised_structure->invert(xreal, x, IN_EPS);
+ project_into_convex(x, pgt, storage.project_into_element);
+ res0 = vect_dist2(pgt->transform(x, G), xreal);
+ }
+
+ if (res < res0) copy(storage.x_ref, x);
}
-
- if (res < res0) copy(storage.x_ref, x);
-}
-
+
/* inversion for non-linear geometric transformations
(Newton on Grad(pgt)(y - pgt(x)) = 0 )
*/
bool geotrans_inv_convex::invert_nonlin(const base_node& xreal,
base_node& x, scalar_type IN_EPS,
- bool &converged, bool throw_except) {
+ bool &converged, bool /* throw_except */) {
converged = true;
find_initial_guess(x, nonlinear_storage, xreal, G, pgt.get(), IN_EPS);
add(pgt->transform(x, G), scaled(xreal, -1.0), nonlinear_storage.diff);
@@ -251,7 +249,8 @@ void find_initial_guess(
if (res > res0) {
add(scaled(nonlinear_storage.x_ref, factor), x);
nonlinear_storage.x_real = pgt->transform(x, G);
- add(nonlinear_storage.x_real, scaled(xreal, -1.0),
nonlinear_storage.diff);
+ add(nonlinear_storage.x_real, scaled(xreal, -1.0),
+ nonlinear_storage.diff);
factor *= 0.5;
}
else {
@@ -265,7 +264,8 @@ void find_initial_guess(
add(scaled(nonlinear_storage.x_ref, -1.0 * factor), x);
project_into_convex(x, pgt.get(),
nonlinear_storage.project_into_element);
nonlinear_storage.x_real = pgt->transform(x, G);
- add(nonlinear_storage.x_real, scaled(xreal, -1.0),
nonlinear_storage.diff);
+ add(nonlinear_storage.x_real, scaled(xreal, -1.0),
+ nonlinear_storage.diff);
res = vect_norm2(nonlinear_storage.diff);
++cnt;
}
diff --git a/src/getfem/bgeot_geotrans_inv.h b/src/getfem/bgeot_geotrans_inv.h
index bd8f4ba..04d0123 100644
--- a/src/getfem/bgeot_geotrans_inv.h
+++ b/src/getfem/bgeot_geotrans_inv.h
@@ -70,7 +70,8 @@ namespace bgeot {
const convex_ind_ct &direct_points_indices,
const stored_point_tab &reference_nodes,
const std::vector<base_node> &real_nodes);
- void invert(const base_node &x_real, base_node &x_ref, scalar_type
IN_EPS) const;
+ void invert(const base_node &x_real, base_node &x_ref,
+ scalar_type IN_EPS) const;
base_matrix K_ref_linear;
base_matrix B_linear;
@@ -90,31 +91,29 @@ namespace bgeot {
base_matrix G, pc, K, B, CS;
pgeometric_trans pgt = nullptr;
scalar_type EPS;
+ nonlinear_storage_struct nonlinear_storage;
+
public:
const base_matrix &get_G() const { return G; }
- geotrans_inv_convex(scalar_type e=10e-12, bool project_into_element =
false) :
+ geotrans_inv_convex(scalar_type e=10e-12, bool project_into_element=false)
:
N(0), P(0), pgt(0), EPS(e)
- {
- nonlinear_storage.project_into_element = project_into_element;
- };
-
- template<class TAB> geotrans_inv_convex(
- const convex<base_node, TAB> &cv,
- pgeometric_trans pgt_,
- scalar_type e=10e-12,
- bool project_into_element = false) : N(0), P(0), pgt(0), EPS(e)
- {
- nonlinear_storage.project_into_element = project_into_element;
- init(cv.points(),pgt_);
+ { this->nonlinear_storage.project_into_element = project_into_element; }
+
+ template<class TAB> geotrans_inv_convex(const convex<base_node, TAB> &cv,
+ pgeometric_trans pgt_,
+ scalar_type e=10e-12,
+ bool project_into_element = false)
+ : N(0), P(0), pgt(0), EPS(e) {
+ this->nonlinear_storage.project_into_element = project_into_element;
+ init(cv.points(),pgt_);
}
- geotrans_inv_convex(
- const std::vector<base_node> &nodes,
- pgeometric_trans pgt_,
- scalar_type e=10e-12,
- bool project_into_element = false) : N(0), P(0), pgt(0), EPS(e)
- {
- nonlinear_storage.project_into_element = project_into_element;
+ geotrans_inv_convex(const std::vector<base_node> &nodes,
+ pgeometric_trans pgt_,
+ scalar_type e=10e-12,
+ bool project_into_element = false)
+ : N(0), P(0), pgt(0), EPS(e) {
+ this->nonlinear_storage.project_into_element = project_into_element;
init(nodes,pgt_);
}
@@ -164,8 +163,6 @@ namespace bgeot {
scalar_type IN_EPS, bool &converged, bool throw_except);
void update_B();
- nonlinear_storage_struct nonlinear_storage;
-
friend class geotrans_inv_convex_bfgs;
};
@@ -190,15 +187,16 @@ namespace bgeot {
// computation of the pseudo inverse
update_B();
} else {
- nonlinear_storage.diff.resize(N);
- nonlinear_storage.x_real.resize(N);
- nonlinear_storage.x_ref.resize(P);
+ this->nonlinear_storage.diff.resize(N);
+ this->nonlinear_storage.x_real.resize(N);
+ this->nonlinear_storage.x_ref.resize(P);
if (pgt->complexity() > 1) {
std::vector<base_node> real_nodes(nodes.begin(), nodes.end());
- nonlinear_storage.plinearised_structure
- = std::make_shared<nonlinear_storage::linearised_structure>(
- pgt->structure()->ind_dir_points(), pgt->geometric_nodes(),
real_nodes);
+ this->nonlinear_storage.plinearised_structure
+ = std::make_shared<nonlinear_storage_struct::linearised_structure>
+ (pgt->structure()->ind_dir_points(), pgt->geometric_nodes(),
+ real_nodes);
}
}
}
@@ -255,7 +253,7 @@ namespace bgeot {
*
* @param itab the indices of points found in the convex.
*
- * @param bruteforce use a brute force search (only for debugging
purposes).
+ * @param bruteforce use a brute force search(only for debugging
purposes).
*
* @return the number of points in the convex (i.e. the size of itab,
* and pftab)