[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Markus Bürg |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Mon, 21 Aug 2017 09:03:05 -0400 (EDT) |
branch: mb-transInversion
commit 8b6c956007208ca0e7aad4286911f63a59be544f
Author: mb <address@hidden>
Date: Mon Aug 21 15:03:02 2017 +0200
Use linear mapping to find initial guess.
---
src/bgeot_geotrans_inv.cc | 203 ++++++++++++++++++++++++++++++++++++++--
src/getfem/bgeot_geotrans_inv.h | 30 ++++++
2 files changed, 225 insertions(+), 8 deletions(-)
diff --git a/src/bgeot_geotrans_inv.cc b/src/bgeot_geotrans_inv.cc
index 9172f26..580bf43 100644
--- a/src/bgeot_geotrans_inv.cc
+++ b/src/bgeot_geotrans_inv.cc
@@ -20,7 +20,12 @@
===========================================================================*/
#include "getfem/bgeot_geotrans_inv.h"
+#include "getfem/bgeot_mesh_structure.h"
#include "gmm/gmm_solver_bfgs.h"
+
+using namespace gmm;
+using namespace std;
+
namespace bgeot
{
bool geotrans_inv_convex::invert(const base_node& n, base_node& n_ref,
@@ -99,6 +104,194 @@ namespace bgeot
}
};
+string element_type_of_pgt(const string &pgt_name) {
+ return pgt_name.substr(3, pgt_name.find("(") - 3);
+}
+
+pair<string, string> get_pgt_names(const string &pgt_name) {
+ auto pos_start = pgt_name.find("(") + 1;
+ auto pos_end = pgt_name.find(")") + 1;
+ auto pgt_name_1 = pgt_name.substr(pos_start, pos_end - pos_start);
+
+ pos_start = pos_end + 1;
+ pos_end = pgt_name.find(")", pos_start) + 1;
+ auto pgt_name_2 = pgt_name.substr(pos_start, pos_end - pos_start);
+
+ return {pgt_name_1, pgt_name_2};
+}
+
+dim_type get_dim_of_pgt(const string &pgt_name)
+{
+ auto element_type = element_type_of_pgt(pgt_name);
+ if ((element_type == "PK") || (element_type == "QK") || (element_type ==
"PRISM")
+ || (element_type == "Q2_INCOMPLETE")) {
+ auto pos_start = pgt_name.find("(") + 1;
+ auto end_symbol = (element_type == "Q2_INCOMPLETE") ? ")" : ",";
+ return stoi(pgt_name.substr(pos_start, pgt_name.find(end_symbol) -
pos_start));
+ }
+ else if (element_type == "PRODUCT") {
+ auto pgt_names = get_pgt_names(pgt_name);
+ return get_dim_of_pgt(pgt_names.first) + get_dim_of_pgt(pgt_names.second);
+ }
+ else
+ {
+ GMM_THROW_DEFAULT("Could not determine type of " + pgt_name);
+ return -1;
+ }
+}
+
+string create_linear_pgt_name(const string &original_pgt) {
+ string linear_pgt;
+ auto element_type = element_type_of_pgt(original_pgt);
+
+ if ((element_type == "PK") || (element_type == "QK") || (element_type ==
"PRISM"))
+ {
+ linear_pgt = original_pgt;
+
+ auto start_pos = linear_pgt.find(",") + 1;
+ linear_pgt.replace(start_pos, linear_pgt.find(")") - start_pos, "1");
+ }
+ else if (element_type == "Q2_INCOMPLETE") {
+ linear_pgt = "GT_QK(" + std::to_string(get_dim_of_pgt(original_pgt)) +
",1)";
+ }
+ else if (element_type == "PRODUCT") {
+ auto pgt_names = get_pgt_names(original_pgt);
+ auto linear_pgt1 = create_linear_pgt_name(pgt_names.first);
+ auto linear_pgt2 = create_linear_pgt_name(pgt_names.second);
+ linear_pgt = "GT_PRODUCT(" + linear_pgt1 + "," + linear_pgt2 + ")";
+ }
+ else GMM_THROW_DEFAULT("Could not determine type of " + original_pgt);
+
+ return linear_pgt;
+}
+
+pgeometric_trans create_linear_pgt(const string &original_pgt) {
+ return geometric_trans_descriptor(create_linear_pgt_name(original_pgt));
+}
+
+base_vector get_degrees_of_pgt(const string &pgt_name, dim_type dim) {
+ auto element_type = element_type_of_pgt(pgt_name);
+ base_vector degrees(dim);
+
+ if ((element_type == "PK") || (element_type == "QK") || (element_type ==
"PRISM")) {
+ auto pos = pgt_name.find(",") + 1;
+ fill(degrees, stoi(pgt_name.substr(pos, pgt_name.find(")") - pos)));
+ }
+ else if (element_type == "Q2_INCOMPLETE") fill(degrees, 2);
+ else if (element_type == "PRODUCT") {
+ auto pgt_names = get_pgt_names(pgt_name);
+ auto dim1 = get_dim_of_pgt(pgt_names.first);
+ auto degrees1 = sub_vector(degrees, sub_interval(0, dim1));
+ copy(get_degrees_of_pgt(pgt_names.first, dim1), degrees1);
+ auto dim2 = get_dim_of_pgt(pgt_names.second);
+ auto degrees2 = sub_vector(degrees, sub_interval(dim1, dim2));
+ copy(get_degrees_of_pgt(pgt_names.second, dim2), degrees2);
+ }
+ else GMM_THROW_DEFAULT("Could not determine type of " + pgt_name);
+
+ return degrees;
+}
+
+vector<size_type> get_linear_nodes_indices(
+ const string &pgt_name, const base_vector °rees)
+{
+ auto dim = degrees.size();
+ auto element_type = element_type_of_pgt(pgt_name);
+ vector<size_type> indices;
+
+ if (element_type == "PK") {
+ indices.push_back(0);
+ indices.push_back(degrees[0]);
+ if (dim > 1) indices.push_back((degrees[0] + 1) * (degrees[1] + 2) / 2 -
1);
+ if (dim > 2) {
+ indices.push_back(
+ (degrees[0] + 1) * (degrees[1] + 2) * ((2 * degrees[2] + 3) / 3 + 1) /
4 - 1);
+ }
+ }
+ else if (element_type == "QK") {
+ indices.push_back(0);
+ indices.push_back(degrees[0]);
+ if (dim > 1) {
+ indices.push_back(degrees[0] * (degrees[1] + 1));
+ indices.push_back((degrees[0] + 1) * (degrees[1] + 1) - 1);
+ }
+ if (dim > 2) {
+ auto nb_nodes_plane = (degrees[0] + 1) * (degrees[1] + 1);
+ indices.push_back(nb_nodes_plane * degrees[2]);
+ indices.push_back(nb_nodes_plane * degrees[2] + degrees[0]);
+ indices.push_back(nb_nodes_plane * degrees[2] + degrees[0] * (degrees[1]
+ 1));
+ indices.push_back(nb_nodes_plane * (degrees[2] + 1) - 1);
+ }
+ }
+ else if (element_type == "PRISM") {
+ indices.push_back(0);
+ indices.push_back(degrees[0]);
+ indices.push_back((degrees[0] + 1) * (degrees[1] + 2) / 2 - 1);
+ indices.push_back((degrees[0] + 1) * (degrees[1] + 2) / 2);
+ indices.push_back((degrees[0] + 1) * (degrees[1] + 2) * (degrees[2] + 1) /
2 - 1);
+ }
+ else if (element_type == "Q2_INCOMPLETE") {
+ indices.push_back(0);
+ indices.push_back(2);
+ if (dim > 1) {
+ indices.push_back(5);
+ indices.push_back(7);
+ }
+ if (dim > 2) {
+ indices.push_back(12);
+ indices.push_back(14);
+ indices.push_back(17);
+ indices.push_back(19);
+ }
+ }
+ else if (element_type == "PRODUCT") {
+ auto pgt_name1 = get_pgt_names(pgt_name).first;
+ auto dim1 = get_dim_of_pgt(pgt_name1);
+ auto indices_plane = get_linear_nodes_indices(pgt_name1, dim1);
+
+ for (auto i : indices_plane) indices.push_back(i);
+
+ for (auto d = dim1; d < dim; ++d) {
+ auto indices_old = indices;
+
+ for (auto i : indices_old) indices.push_back(degrees[d] *
indices_old.size() + i);
+ }
+ }
+
+ return indices;
+}
+
+void project_into_convex(base_node &x, const geometric_trans &geo_trans) {
+ auto dim = geo_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 nb_simplices =
geo_trans.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;
+ }
+ }
+}
+
+vector<size_type> get_linear_nodes_indices(const string &pgt_name, dim_type
dim) {
+ return get_linear_nodes_indices(pgt_name, get_degrees_of_pgt(pgt_name, dim));
+}
+
/* inversion for non-linear geometric transformations
(Newton on Grad(pgt)(y - pgt(x)) = 0 )
*/
@@ -110,14 +303,8 @@ namespace bgeot
converged = true;
base_node xn(P), y, z,x0;
/* find an initial guess */
- x0 = (pgt->geometric_nodes())[0]; copy(mat_col(G, 0), y);
- scalar_type d = gmm::vect_dist2_sqr(y, xreal);
- for (size_type j = 1; j < pgt->nb_points(); ++j) {
- scalar_type d2 = gmm::vect_dist2_sqr(mat_col(G, j), xreal);
- if (d2 < d)
- { d = d2; x0 = pgt->geometric_nodes()[j]; copy(mat_col(G, j), y); }
- }
- x = x0;
+ nonlinear_storage.plinear_inversion->invert(xreal, x);
+ project_into_convex(x, *nonlinear_storage.plinear_inversion->pgt);
base_node vres(P);
base_node rn(xreal); rn -= y;
diff --git a/src/getfem/bgeot_geotrans_inv.h b/src/getfem/bgeot_geotrans_inv.h
index 17cf465..198365c 100644
--- a/src/getfem/bgeot_geotrans_inv.h
+++ b/src/getfem/bgeot_geotrans_inv.h
@@ -57,6 +57,14 @@
#include "bgeot_kdtree.h"
namespace bgeot {
+ class geotrans_inv_convex;
+
+ struct nonlinear_storage {
+ base_node diff;
+ base_node x_real;
+ base_node x_ref;
+ std::shared_ptr<geotrans_inv_convex> plinear_inversion;
+ };
/**
does the inversion of the geometric transformation for a given convex
*/
@@ -120,9 +128,16 @@ namespace bgeot {
bool invert_nonlin(const base_node& n, base_node& n_ref,
scalar_type IN_EPS, bool &converged, bool throw_except);
void update_B();
+
+ nonlinear_storage nonlinear_storage;
+
friend class geotrans_inv_convex_bfgs;
};
+ pgeometric_trans create_linear_pgt(const std::string &original_pgt);
+ std::vector<size_type> get_linear_nodes_indices(const std::string &pgt_name,
dim_type dim);
+
+
template<class TAB>
void geotrans_inv_convex::init(const TAB &nodes, pgeometric_trans pgt_) {
bool geotrans_changed = (pgt != pgt_); if (geotrans_changed) pgt = pgt_;
@@ -143,6 +158,21 @@ namespace bgeot {
}
// computation of the pseudo inverse
update_B();
+ } else {
+ nonlinear_storage.diff.resize(P);
+ nonlinear_storage.x_real.resize(P);
+ nonlinear_storage.x_ref.resize(P);
+
+ auto name_pgt = name_of_geometric_trans(pgt);
+ auto linear_pgt = create_linear_pgt(name_pgt);
+ std::vector<base_node> linear_nodes;
+
+ for (auto &&i : get_linear_nodes_indices(name_pgt, P)) {
+ linear_nodes.push_back(nodes[i]);
+ }
+
+ nonlinear_storage.plinear_inversion
+ = std::make_shared<geotrans_inv_convex>(linear_nodes, linear_pgt);
}
}