getfem-commits
[Top][All Lists]
Advanced

[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: Fri, 1 Sep 2017 05:17:29 -0400 (EDT)

branch: mb-transInversion
commit f237af9e65923bf720fdf4d6211e793d7adfddf0
Author: mb <address@hidden>
Date:   Fri Sep 1 11:17:25 2017 +0200

    Use dir_points for creation of linearised transformation.
---
 src/bgeot_geotrans_inv.cc       | 197 ++++++++--------------------------------
 src/getfem/bgeot_geotrans_inv.h |  36 +++++---
 2 files changed, 60 insertions(+), 173 deletions(-)

diff --git a/src/bgeot_geotrans_inv.cc b/src/bgeot_geotrans_inv.cc
index 22e66a1..b4422e0 100644
--- a/src/bgeot_geotrans_inv.cc
+++ b/src/bgeot_geotrans_inv.cc
@@ -108,167 +108,52 @@ bool point_in_convex(
     }
   };
 
-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)
+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,
+  dim_type P) : diff(real_nodes.begin()->size()), 
diff_ref(real_nodes.begin()->size())
 {
-  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);
+  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]);
   }
-  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);
+  auto N = direct_points_ref.begin()->size();
+  B_linear.base_resize(N, n_points - 1);
+  K_ref_linear.base_resize(N, n_points - 1);
+  P_linear = direct_points[0];
+  P_ref_linear = direct_points_ref[0];
 
-  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 + ")";
+  for (size_type i = 1; i < n_points; ++i) {
+    add(direct_points[i], scaled(P_linear, -1.0), mat_col(B_linear, i - 1));
+    add(direct_points_ref[i], scaled(P_ref_linear, -1.0), 
mat_col(K_ref_linear, i - 1));
   }
-  else GMM_THROW_DEFAULT("Could not determine type of " + original_pgt);
-
-  return linear_pgt;
-}
-
-pgeometric_trans create_linear_pgt(pgeometric_trans poriginal_gt) {
-  auto is_axisymmetric = is_torus_geom_trans(poriginal_gt);
-  auto name_original_pgt = name_of_geometric_trans(poriginal_gt);
-  auto plinear_pgt = 
geometric_trans_descriptor(create_linear_pgt_name(name_original_pgt));
-
-  return is_axisymmetric ? std::make_shared<const 
torus_geom_trans>(move(plinear_pgt))
-                         : plinear_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)));
+  if (N == P) lu_inverse(B_linear);
+  else {
+    base_matrix K_linear(n_points - 1, n_points - 1);
+    mult(transposed(B_linear), B_linear, K_linear);
+    lu_inverse(K_linear);
+    base_matrix temp(n_points - 1, n_points - 1);
+    mult(B_linear, K_linear, temp);
+    copy(temp, B_linear);
   }
-  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 &degrees, dim_type dim)
+void nonlinear_storage::linearised_structure::invert(
+  const base_node &x_real, base_node &x_ref, scalar_type IN_EPS) const
 {
-  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) / 2 + degrees[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, get_degrees_of_pgt(pgt_name1, dim1), 
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;
+  add(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) {
@@ -329,8 +214,8 @@ void find_initial_guess(
 
   res0 = std::numeric_limits<scalar_type>::max();
 
-  if (storage.plinear_inversion != nullptr) {
-    storage.plinear_inversion->invert(xreal, x, IN_EPS);
+  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);
   }
@@ -338,12 +223,6 @@ void find_initial_guess(
   if (res < res0) copy(storage.x_ref, x);
 }
 
-vector<size_type> get_linear_nodes_indices(pgeometric_trans pgt) {
-  auto pgt_name = name_of_geometric_trans(pgt);
-  auto original_dim = is_torus_geom_trans(pgt) ? 2 : pgt->dim();
-  return get_linear_nodes_indices(
-    pgt_name, get_degrees_of_pgt(pgt_name, pgt->dim()), original_dim);
-}
 
   /* inversion for non-linear geometric transformations 
      (Newton on Grad(pgt)(y - pgt(x)) = 0 )
diff --git a/src/getfem/bgeot_geotrans_inv.h b/src/getfem/bgeot_geotrans_inv.h
index 419da48..490d48e 100644
--- a/src/getfem/bgeot_geotrans_inv.h
+++ b/src/getfem/bgeot_geotrans_inv.h
@@ -64,7 +64,24 @@ namespace bgeot {
     base_node x_real;
     base_node x_ref;
     bool project_into_element;
-    std::shared_ptr<geotrans_inv_convex> plinear_inversion;
+
+    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,
+        dim_type P);
+      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;
+      base_node P_linear;
+      base_node P_ref_linear;
+      mutable base_node diff;
+      mutable base_node diff_ref;
+    };
+
+    std::shared_ptr<linearised_structure> plinearised_structure = nullptr;
   };
   /** 
       does the inversion of the geometric transformation for a given convex
@@ -153,10 +170,6 @@ namespace bgeot {
     friend class geotrans_inv_convex_bfgs;
   };
 
-  pgeometric_trans create_linear_pgt(pgeometric_trans poriginal_gt);
-  std::vector<size_type> get_linear_nodes_indices(pgeometric_trans pgt);
-
-
   template<class TAB>
   void geotrans_inv_convex::init(const TAB &nodes,  pgeometric_trans pgt_) {
     bool geotrans_changed = (pgt != pgt_); if (geotrans_changed) pgt = pgt_;
@@ -183,15 +196,10 @@ namespace bgeot {
       nonlinear_storage.x_ref.resize(P);
 
       if (pgt->complexity() > 1) {
-        auto plinear_pgt = create_linear_pgt(pgt);
-        std::vector<base_node> linear_nodes;
-
-        for (auto &&i : get_linear_nodes_indices(pgt)) {
-          linear_nodes.push_back(nodes[i]);
-        }
-
-        nonlinear_storage.plinear_inversion
-          = std::make_shared<geotrans_inv_convex>(linear_nodes, plinear_pgt);
+        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, P);
       }
     }
   }



reply via email to

[Prev in Thread] Current Thread [Next in Thread]