GetFEM  5.4.2
bgeot_geotrans_inv.cc
1 /*===========================================================================
2 
3  Copyright (C) 2000-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
24 #include "getfem/bgeot_torus.h"
25 #include "gmm/gmm_solver_bfgs.h"
26 
27 
28 namespace bgeot
29 {
30  void project_into_convex(base_node &x, const pgeometric_trans pgt) {
31 
32  for (auto &coord : x) {
33  if (coord < 0.0) coord = 0.0;
34  if (coord > 1.0) coord = 1.0;
35  }
36 
37 
38  auto pgt_torus = std::dynamic_pointer_cast<const torus_geom_trans>(pgt);
39  const pgeometric_trans
40  orig_pgt = pgt_torus ? pgt_torus->get_original_transformation()
41  : pgt;
42 
43  auto pbasic_convex_ref = basic_convex_ref(orig_pgt->convex_ref());
44  auto nb_simplices = pbasic_convex_ref->simplexified_convex()->nb_convex();
45 
46  if (nb_simplices == 1) { // simplex
47  auto sum_coordinates = 0.0;
48  for (const auto &coord : x) sum_coordinates += coord;
49 
50  if (sum_coordinates > 1.0) gmm::scale(x, 1.0 / sum_coordinates);
51  }
52  else if (pgt->dim() == 3 && nb_simplices != 4) { // prism
53  auto sum_coordinates = x[0] + x[1];
54  if (sum_coordinates > 1.0) {
55  x[0] /= sum_coordinates;
56  x[1] /= sum_coordinates;
57  }
58  }
59  }
60 
62  scalar_type IN_EPS,
63  bool project_into_element) {
64  bool converged = true;
65  return invert(n, n_ref, converged, IN_EPS, project_into_element);
66  }
67 
69  bool &converged,
70  scalar_type IN_EPS,
71  bool project_into_element) {
72  assert(pgt);
73  n_ref.resize(pgt->structure()->dim());
74  converged = true;
75  if (pgt->is_linear())
76  return invert_lin(n, n_ref, IN_EPS);
77  else
78  return invert_nonlin(n, n_ref, IN_EPS, converged, false,
79  project_into_element);
80  }
81 
82  /* inversion for linear geometric transformations */
83  bool geotrans_inv_convex::invert_lin(const base_node& n, base_node& n_ref,
84  scalar_type IN_EPS) {
85  base_node y(n); for (size_type i=0; i < N; ++i) y[i] -= G(i,0);
86  mult(transposed(B), y, n_ref);
87  y = pgt->transform(n_ref, G);
88  add(gmm::scaled(n, -1.0), y);
89 
90  return (pgt->convex_ref()->is_in(n_ref) < IN_EPS) &&
91  (gmm::vect_norm2(y) < IN_EPS);
92  }
93 
94  void geotrans_inv_convex::update_B() {
95  if (P != N) {
96  pgt->compute_K_matrix(G, pc, K);
97  gmm::mult(gmm::transposed(K), K, CS);
98  bgeot::lu_inverse(&(*(CS.begin())), P);
99  gmm::mult(K, CS, B);
100  }
101  else {
102  // L'inversion peut être optimisée par le non calcul global de B
103  // et la resolution d'un système linéaire.
104  base_matrix KT(K.nrows(), K.ncols());
105  pgt->compute_K_matrix(G, pc, KT);
106  gmm::copy(gmm::transposed(KT), K);
107  gmm::copy(K, B);
108  bgeot::lu_inverse(&(*(K.begin())), P); B.swap(K);
109  }
110  }
111 
112  class geotrans_inv_convex_bfgs {
113  geotrans_inv_convex &gic;
114  base_node xreal;
115  public:
116  geotrans_inv_convex_bfgs(geotrans_inv_convex &gic_,
117  const base_node &xr) : gic(gic_), xreal(xr) {}
118  scalar_type operator()(const base_node& x) const {
119  base_node r = gic.pgt->transform(x, gic.G) - xreal;
120  return gmm::vect_norm2_sqr(r)/2.;
121  }
122  void operator()(const base_node& x, base_small_vector& gr) const {
123  gic.pgt->poly_vector_grad(x, gic.pc);
124  gic.update_B();
125  base_node r = gic.pgt->transform(x, gic.G) - xreal;
126  gr.resize(x.size());
127  gmm::mult(gmm::transposed(gic.K), r, gr);
128  }
129  };
130 
131  void geotrans_inv_convex::update_linearization() {
132 
133  const convex_ind_ct &dir_pt_ind = pgt->structure()->ind_dir_points();
134  const stored_point_tab &ref_nodes = pgt->geometric_nodes();
135 
136  has_linearized_approx = true;
137 
138  auto n_points = dir_pt_ind.size();
139  auto N_ref = ref_nodes.begin()->size();
140 
141  std::vector<base_node> dir_pts, dir_pts_ref;
142  for (auto i : dir_pt_ind) {
143  dir_pts.push_back(base_node(N));
144  gmm::copy(mat_col(G, i), dir_pts.back());
145  dir_pts_ref.push_back(ref_nodes[i]);
146  }
147 
148  base_matrix K_lin(N, n_points - 1),
149  B_transp_lin(n_points - 1, N),
150  K_ref_lin(N_ref, n_points - 1);
151 
152  P_lin = dir_pts[0];
153  P_ref_lin = dir_pts_ref[0];
154 
155  for (size_type i = 1; i < n_points; ++i) {
156  add(dir_pts[i], gmm::scaled(P_lin, -1.0), mat_col(K_lin, i - 1));
157  add(dir_pts_ref[i], gmm::scaled(P_ref_lin, -1.0),
158  mat_col(K_ref_lin, i - 1));
159  }
160 
161  if (K_lin.nrows() == K_lin.ncols()) {
162  lu_inverse(K_lin);
163  gmm::copy(K_lin, B_transp_lin);
164  }
165  else {
166  base_matrix temp(n_points - 1, n_points - 1);
167  mult(transposed(K_lin), K_lin, temp);
168  lu_inverse(temp);
169  mult(temp, transposed(K_lin), B_transp_lin);
170  }
171 
172  K_ref_B_transp_lin.base_resize(N_ref, N);
173  mult(K_ref_lin, B_transp_lin, K_ref_B_transp_lin);
174  }
175 
176 
177  /* inversion for non-linear geometric transformations
178  (Newton on Grad(pgt)(y - pgt(x)) = 0 )
179  */
180  bool geotrans_inv_convex::invert_nonlin(const base_node& xreal,
181  base_node& x, scalar_type IN_EPS,
182  bool &converged,
183  bool /* throw_except */,
184  bool project_into_element) {
185  converged = true;
186  base_node x0_ref(P), diff(N);
187 
188  { // find initial guess
189  x0_ref = pgt->geometric_nodes()[0];
190  scalar_type res = gmm::vect_dist2(mat_col(G, 0), xreal);
191  for (size_type j = 1; j < pgt->nb_points(); ++j) {
192  scalar_type res0 = gmm::vect_dist2(mat_col(G, j), xreal);
193  if (res0 < res) {
194  res = res0;
195  x0_ref = pgt->geometric_nodes()[j];
196  }
197  }
198 
199  scalar_type res0 = std::numeric_limits<scalar_type>::max();
200  if (has_linearized_approx) {
201 
202  add(xreal, gmm::scaled(P_lin, -1.0), diff);
203  mult(K_ref_B_transp_lin, diff, x);
204  gmm::add(P_ref_lin, x);
205 
206  if (project_into_element) project_into_convex(x, pgt);
207  res0 = gmm::vect_dist2(pgt->transform(x, G), xreal);
208  }
209 
210  if (res < res0) gmm::copy(x0_ref, x);
211  if (res < IN_EPS)
212  x *= 0.999888783; // For pyramid element to avoid the singularity
213  }
214 
215  add(pgt->transform(x, G), gmm::scaled(xreal, -1.0), diff);
216  scalar_type res = gmm::vect_norm2(diff);
217  scalar_type res0 = std::numeric_limits<scalar_type>::max();
218  scalar_type factor = 1.0;
219 
220  base_node x0_real(N);
221  while (res > IN_EPS/100.) {
222  if ((gmm::abs(res - res0) < IN_EPS/100.) || (factor < IN_EPS)) {
223  converged = false;
224  return (pgt->convex_ref()->is_in(x) < IN_EPS) && (res < IN_EPS);
225  }
226  if (res > res0) {
227  add(gmm::scaled(x0_ref, factor), x);
228  x0_real = pgt->transform(x, G);
229  add(x0_real, gmm::scaled(xreal, -1.0), diff);
230  factor *= 0.5;
231  }
232  else {
233  if (factor < 1.0-IN_EPS) factor *= 2.0;
234  res0 = res;
235  }
236  pgt->poly_vector_grad(x, pc);
237  update_B();
238  mult(transposed(B), diff, x0_ref);
239  add(gmm::scaled(x0_ref, -factor), x);
240  if (project_into_element) project_into_convex(x, pgt);
241  x0_real = pgt->transform(x, G);
242  add(x0_real, gmm::scaled(xreal, -1.0), diff);
243  res = gmm::vect_norm2(diff);
244  }
245  return (pgt->convex_ref()->is_in(x) < IN_EPS) && (res < IN_EPS);
246  }
247 
248 } /* end of namespace bgeot. */
bgeot_torus.h
Provides mesh of torus.
bgeot::basic_convex_ref
pconvex_ref basic_convex_ref(pconvex_ref cvr)
return the associated order 1 reference convex.
Definition: bgeot_convex_ref.h:138
gmm_solver_bfgs.h
Implements BFGS (Broyden, Fletcher, Goldfarb, Shanno) algorithm.
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::vect_dist2
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist2(const V1 &v1, const V2 &v2)
Euclidean distance between two vectors.
Definition: gmm_blas.h:597
gmm::vect_norm2_sqr
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2_sqr(const V &v)
squared Euclidean norm of a vector.
Definition: gmm_blas.h:544
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
bgeot::small_vector< scalar_type >
bgeot_geotrans_inv.h
Inversion of geometric transformations.
bgeot_mesh_structure.h
Mesh structure definition.
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:27
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
bgeot::geotrans_inv_convex::invert
bool invert(const base_node &n, base_node &n_ref, scalar_type IN_EPS=1e-12, bool project_into_element=false)
given the node on the real element, returns the node on the reference element (even if it is outside ...
Definition: bgeot_geotrans_inv.cc:61
gmm::add
void add(const L1 &l1, L2 &l2)
*â€/
Definition: gmm_blas.h:1268