GetFEM  5.4.1
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 
31  bool point_in_convex(const geometric_trans &geoTrans,
32  const base_node &x,
33  scalar_type res,
34  scalar_type IN_EPS) {
35  // Test un peu sev�re peut-�tre en ce qui concerne res.
36  return (geoTrans.convex_ref()->is_in(x) < IN_EPS) && (res < IN_EPS);
37  }
38 
39  void project_into_convex(base_node &x, const pgeometric_trans pgt) {
40 
41  for (auto &coord : x) {
42  if (coord < 0.0) coord = 0.0;
43  if (coord > 1.0) coord = 1.0;
44  }
45 
46 
47  auto pgt_torus = std::dynamic_pointer_cast<const torus_geom_trans>(pgt);
48  const pgeometric_trans
49  orig_pgt = pgt_torus ? pgt_torus->get_original_transformation()
50  : pgt;
51 
52  auto pbasic_convex_ref = basic_convex_ref(orig_pgt->convex_ref());
53  auto nb_simplices = pbasic_convex_ref->simplexified_convex()->nb_convex();
54 
55  if (nb_simplices == 1) { // simplex
56  auto sum_coordinates = 0.0;
57  for (const auto &coord : x) sum_coordinates += coord;
58 
59  if (sum_coordinates > 1.0) gmm::scale(x, 1.0 / sum_coordinates);
60  }
61  else if (pgt->dim() == 3 && nb_simplices != 4) { // prism
62  auto sum_coordinates = x[0] + x[1];
63  if (sum_coordinates > 1.0) {
64  x[0] /= sum_coordinates;
65  x[1] /= sum_coordinates;
66  }
67  }
68  }
69 
71  scalar_type IN_EPS,
72  bool project_into_element) {
73  bool converged = true;
74  return invert(n, n_ref, converged, IN_EPS, project_into_element);
75  }
76 
78  bool &converged,
79  scalar_type IN_EPS,
80  bool project_into_element) {
81  assert(pgt);
82  n_ref.resize(pgt->structure()->dim());
83  converged = true;
84  if (pgt->is_linear())
85  return invert_lin(n, n_ref, IN_EPS);
86  else
87  return invert_nonlin(n, n_ref, IN_EPS, converged, false,
88  project_into_element);
89  }
90 
91  /* inversion for linear geometric transformations */
92  bool geotrans_inv_convex::invert_lin(const base_node& n, base_node& n_ref,
93  scalar_type IN_EPS) {
94  base_node y(n); for (size_type i=0; i < N; ++i) y[i] -= G(i,0);
95  mult(transposed(B), y, n_ref);
96  y = pgt->transform(n_ref, G);
97  add(gmm::scaled(n, -1.0), y);
98 
99  return point_in_convex(*pgt, n_ref, gmm::vect_norm2(y), IN_EPS);
100  }
101 
102  void geotrans_inv_convex::update_B() {
103  if (P != N) {
104  pgt->compute_K_matrix(G, pc, K);
105  gmm::mult(gmm::transposed(K), K, CS);
106  bgeot::lu_inverse(&(*(CS.begin())), P);
107  gmm::mult(K, CS, B);
108  }
109  else {
110  // L'inversion peut �tre optimis�e par le non calcul global de B
111  // et la resolution d'un syst�me lin�aire.
112  base_matrix KT(K.nrows(), K.ncols());
113  pgt->compute_K_matrix(G, pc, KT);
114  gmm::copy(gmm::transposed(KT), K);
115  gmm::copy(K, B);
116  bgeot::lu_inverse(&(*(K.begin())), P); B.swap(K);
117  }
118  }
119 
120  class geotrans_inv_convex_bfgs {
121  geotrans_inv_convex &gic;
122  base_node xreal;
123  public:
124  geotrans_inv_convex_bfgs(geotrans_inv_convex &gic_,
125  const base_node &xr) : gic(gic_), xreal(xr) {}
126  scalar_type operator()(const base_node& x) const {
127  base_node r = gic.pgt->transform(x, gic.G) - xreal;
128  return gmm::vect_norm2_sqr(r)/2.;
129  }
130  void operator()(const base_node& x, base_small_vector& gr) const {
131  gic.pgt->poly_vector_grad(x, gic.pc);
132  gic.update_B();
133  base_node r = gic.pgt->transform(x, gic.G) - xreal;
134  gr.resize(x.size());
135  gmm::mult(gmm::transposed(gic.K), r, gr);
136  }
137  };
138 
139  void geotrans_inv_convex::update_linearization() {
140 
141  const convex_ind_ct &dir_pt_ind = pgt->structure()->ind_dir_points();
142  const stored_point_tab &ref_nodes = pgt->geometric_nodes();
143 
144  has_linearized_approx = true;
145 
146  auto n_points = dir_pt_ind.size();
147  auto N_ref = ref_nodes.begin()->size();
148 
149  std::vector<base_node> dir_pts, dir_pts_ref;
150  for (auto i : dir_pt_ind) {
151  dir_pts.push_back(base_node(N));
152  gmm::copy(mat_col(G, i), dir_pts.back());
153  dir_pts_ref.push_back(ref_nodes[i]);
154  }
155 
156  base_matrix K_lin(N, n_points - 1),
157  B_transp_lin(n_points - 1, N),
158  K_ref_lin(N_ref, n_points - 1);
159 
160  P_lin = dir_pts[0];
161  P_ref_lin = dir_pts_ref[0];
162 
163  for (size_type i = 1; i < n_points; ++i) {
164  add(dir_pts[i], gmm::scaled(P_lin, -1.0), mat_col(K_lin, i - 1));
165  add(dir_pts_ref[i], gmm::scaled(P_ref_lin, -1.0),
166  mat_col(K_ref_lin, i - 1));
167  }
168 
169  if (K_lin.nrows() == K_lin.ncols()) {
170  lu_inverse(K_lin);
171  gmm::copy(K_lin, B_transp_lin);
172  }
173  else {
174  base_matrix temp(n_points - 1, n_points - 1);
175  mult(transposed(K_lin), K_lin, temp);
176  lu_inverse(temp);
177  mult(temp, transposed(K_lin), B_transp_lin);
178  }
179 
180  K_ref_B_transp_lin.base_resize(N_ref, N);
181  mult(K_ref_lin, B_transp_lin, K_ref_B_transp_lin);
182  }
183 
184 
185  /* inversion for non-linear geometric transformations
186  (Newton on Grad(pgt)(y - pgt(x)) = 0 )
187  */
188  bool geotrans_inv_convex::invert_nonlin(const base_node& xreal,
189  base_node& x, scalar_type IN_EPS,
190  bool &converged,
191  bool /* throw_except */,
192  bool project_into_element) {
193  converged = true;
194 
195  base_node x0_ref(P), diff(N);
196 
197  { // find initial guess
198  x0_ref = pgt->geometric_nodes()[0];
199  scalar_type res = gmm::vect_dist2(mat_col(G, 0), xreal);
200  for (size_type j = 1; j < pgt->nb_points(); ++j) {
201  scalar_type res0 = gmm::vect_dist2(mat_col(G, j), xreal);
202  if (res0 < res) {
203  res = res0;
204  x0_ref = pgt->geometric_nodes()[j];
205  }
206  }
207 
208  scalar_type res0 = std::numeric_limits<scalar_type>::max();
209  if (has_linearized_approx) {
210 
211  add(xreal, gmm::scaled(P_lin, -1.0), diff);
212  mult(K_ref_B_transp_lin, diff, x);
213  gmm::add(P_ref_lin, x);
214 
215  if (project_into_element) project_into_convex(x, pgt);
216  res0 = gmm::vect_dist2(pgt->transform(x, G), xreal);
217  }
218 
219  if (res < res0) gmm::copy(x0_ref, x);
220  if (res < IN_EPS)
221  x *= 0.999888783; // For pyramid element to avoid the singularity
222  }
223 
224  add(pgt->transform(x, G), gmm::scaled(xreal, -1.0), diff);
225  auto res = gmm::vect_norm2(diff);
226  auto res0 = std::numeric_limits<scalar_type>::max();
227  double factor = 1.0;
228 
229  base_node x0_real(N);
230  while (res > IN_EPS) {
231  if ((gmm::abs(res - res0) < IN_EPS) || (factor < IN_EPS)) {
232  converged = false;
233  return point_in_convex(*pgt, x, res, IN_EPS);
234  }
235  if (res > res0) {
236  add(gmm::scaled(x0_ref, factor), x);
237  x0_real = pgt->transform(x, G);
238  add(x0_real, gmm::scaled(xreal, -1.0), diff);
239  factor *= 0.5;
240  }
241  else {
242  if (factor < 1.0-IN_EPS) factor = 2.0;
243  res0 = res;
244  }
245  pgt->poly_vector_grad(x, pc);
246  update_B();
247  mult(transposed(B), diff, x0_ref);
248  add(gmm::scaled(x0_ref, -1.0 * factor), x);
249  if (project_into_element) project_into_convex(x, pgt);
250  x0_real = pgt->transform(x, G);
251  add(x0_real, gmm::scaled(xreal, -1.0), diff);
252  res = gmm::vect_norm2(diff);
253  }
254 
255  return point_in_convex(*pgt, x, res, IN_EPS);
256  }
257 
258 } /* 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:70
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1268