GetFEM  5.4.2
bgeot_geotrans_inv.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2000-2020 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file bgeot_geotrans_inv.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>
34  @date December 20, 2000.
35  @brief Inversion of geometric transformations.
36 
37  Inversion means: given a set of convexes and a point, find:
38 
39  - a subset of candidate convexes, which are likely to contain the
40  point (using bgeot::kdtree).
41 
42  - on these candidate convexes, invert the geometric
43  transformation, i.e. find the corresponding coordinates on the
44  reference element.
45 
46  Inversion of a geometric transformation is not a trivial task,
47  especially with non-linear geometric transformations. This is the
48  central part of interpolation routines from a getfem::mesh_fem onto
49  another.
50 */
51 
52 #ifndef BGEOT_GEOTRANS_INV_H__
53 #define BGEOT_GEOTRANS_INV_H__
54 
55 #include "bgeot_geometric_trans.h"
56 #include "bgeot_small_vector.h"
57 #include "bgeot_kdtree.h"
58 
59 namespace bgeot {
60 
61  /**
62  does the inversion of the geometric transformation for a given convex
63  */
65  size_type N, P;
66  base_matrix G, pc, K, B, CS;
67  pgeometric_trans pgt = nullptr;
68  scalar_type EPS;
69 
70  bool has_linearized_approx = false;
71  base_matrix K_ref_B_transp_lin;
72  base_node P_lin, P_ref_lin;
73 
74  public:
75  const base_matrix &get_G() const { return G; }
76 
77  geotrans_inv_convex(scalar_type e=10e-12)
78  : N(0), P(0), pgt(0), EPS(e) {}
79 
80  template<class TAB> geotrans_inv_convex(const convex<base_node, TAB> &cv,
81  pgeometric_trans pgt_,
82  scalar_type e=10e-12)
83  : N(0), P(0), pgt(0), EPS(e)
84  {
85  init(cv.points(), pgt_);
86  }
87 
88  geotrans_inv_convex(const std::vector<base_node> &nodes,
89  pgeometric_trans pgt_,
90  scalar_type e=10e-12)
91  : N(0), P(0), pgt(0), EPS(e)
92  {
93  init(nodes, pgt_);
94  }
95 
96  template<class TAB> void init(const TAB &nodes, pgeometric_trans pgt_);
97 
98  /**
99  given the node on the real element, returns the node on the
100  reference element (even if it is outside of the ref. convex).
101 
102  If the geometric transformation is not invertible at point n,
103  an exception is thrown.
104 
105  @return true if the n is inside the convex
106 
107  @param n node on the real element
108 
109  @param n_ref computed node on the reference convex
110 
111  @param IN_EPS a threshold.
112  */
113  bool invert(const base_node& n, base_node& n_ref,
114  scalar_type IN_EPS=1e-12, bool project_into_element=false);
115 
116  /**
117  given the node on the real element, returns the node
118  on the reference element (even if it is outside of the ref. convex).
119 
120  This version will not throw an exception if the geometric
121  transformation is not invertible at point n.
122 
123  @return true if the n is inside the convex
124 
125  @param n node on the real element
126 
127  @param n_ref computed node on the reference convex
128 
129  @param converged on output, will be set to true if the
130  geometric transformation could be inverted.
131 
132  @param IN_EPS a threshold.
133  */
134  bool invert(const base_node& n, base_node& n_ref, bool &converged,
135  scalar_type IN_EPS=1e-12, bool project_into_element=false);
136  private:
137  bool invert_lin(const base_node& n, base_node& n_ref, scalar_type IN_EPS);
138  bool invert_nonlin(const base_node& n, base_node& n_ref,
139  scalar_type IN_EPS, bool &converged, bool throw_except,
140  bool project_into_element);
141  void update_B();
142  void update_linearization();
143 
144  friend class geotrans_inv_convex_bfgs;
145  };
146 
147  template<class TAB>
148  void geotrans_inv_convex::init(const TAB &nodes, pgeometric_trans pgt_) {
149  bool geotrans_changed = (pgt != pgt_); if (geotrans_changed) pgt = pgt_;
150  GMM_ASSERT3(!nodes.empty(), "empty points!");
151  if (N != nodes[0].size())
152  { N = nodes[0].size(); geotrans_changed = true; }
153  if (geotrans_changed) {
154  P = pgt->structure()->dim();
155  pc.resize(pgt->nb_points() , P);
156  K.resize(N,P); B.resize(N,P); CS.resize(P,P);
157  G.resize(N, pgt->nb_points());
158  }
159  vectors_to_base_matrix(G, nodes);
160  if (pgt->is_linear()) {
161  if (geotrans_changed) {
162  base_node Dummy(P);
163  pgt->poly_vector_grad(Dummy, pc);
164  }
165  // computation of the pseudo inverse
166  update_B();
167  } else {
168  if (pgt->complexity() > 1)
169  update_linearization();
170  }
171  }
172 
173 
174  /**
175  handles the geometric inversion for a given (supposedly quite large)
176  set of points
177  */
179  {
180  protected :
181  mutable kdtree tree;
182  scalar_type EPS;
184  public :
185  void clear(void) { tree.clear(); }
186  /// Add the points contained in c to the list of points.
187  template<class CONT> void add_points(const CONT &c) {
188  tree.reserve(std::distance(c.begin(),c.end()));
189  typename CONT::const_iterator it = c.begin(), ite = c.end();
190  for (; it != ite; ++it) tree.add_point(*it);
191  }
192 
193  /// Number of points.
194  size_type nb_points(void) const { return tree.nb_points(); }
195  /// Add point p to the list of points.
196  size_type add_point(base_node p) { return tree.add_point(p); }
197  void add_point_with_id(base_node p,size_type id)
198  { tree.add_point_with_id(p,id); }
199 
200  /// Find all the points present in the box between min and max.
202  const base_node &min,
203  const base_node &max) const {
204  tree.points_in_box(ipts, min, max);
205  return ipts.size();
206  }
207 
208  /** Search all the points in the convex cv, which is the transformation
209  * of the convex cref via the geometric transformation pgt.
210  *
211  * IMPORTANT : It is assumed that the whole convex is include in the
212  * minmax box of its nodes times a factor 1.2. If the transformation is
213  * linear, the factor is 1.0.
214  *
215  * @param cv the convex points (as given by getfem_mesh::convex(ic)).
216  *
217  * @param pgt the geometric trans (as given by
218  * getfem_mesh::trans_of_convex(ic)).
219  *
220  * @param pftab container for the coordinates of points in the reference
221  * convex (should be of size nb_points())
222  *
223  * @param itab the indices of points found in the convex.
224  *
225  * @param bruteforce use a brute force search(only for debugging purposes).
226  *
227  * @return the number of points in the convex (i.e. the size of itab,
228  * and pftab)
229  */
230  template<class TAB, class CONT1, class CONT2>
232  pgeometric_trans pgt,
233  CONT1 &pftab, CONT2 &itab,
234  bool bruteforce=false);
235 
236  geotrans_inv(scalar_type EPS_ = 10E-12) : EPS(EPS_) {}
237  };
238 
239 
240 
241  template<class TAB, class CONT1, class CONT2>
243  pgeometric_trans pgt,
244  CONT1 &pftab, CONT2 &itab,
245  bool bruteforce) {
246  base_node min, max; /* bound of the box enclosing the convex */
247  size_type nbpt = 0; /* nb of points in the convex */
248  kdtree_tab_type boxpts;
249  bounding_box(min, max, cv.points(), pgt);
250  for (size_type k=0; k < min.size(); ++k) { min[k] -= EPS; max[k] += EPS; }
251  gic.init(cv.points(),pgt);
252  /* get the points in a box enclosing the convex */
253  if (!bruteforce) points_in_box(boxpts, min, max);
254  else boxpts = tree.points();
255  /* and invert the geotrans, and check if the obtained point is
256  inside the reference convex */
257  for (size_type l = 0; l < boxpts.size(); ++l) {
258  // base_node pt_ref;
259  if (gic.invert(boxpts[l].n, pftab[nbpt], EPS)) {
260  itab[nbpt++] = boxpts[l].i;
261  }
262  }
263  return nbpt;
264  }
265 
266 
267 
268 
269 
270 } /* end of namespace bgeot. */
271 
272 
273 #endif /* BGEOT_GEOMETRIC_TRANSFORMATION_H__ */
bgeot::geotrans_inv::points_in_convex
size_type points_in_convex(const convex< base_node, TAB > &cv, pgeometric_trans pgt, CONT1 &pftab, CONT2 &itab, bool bruteforce=false)
Search all the points in the convex cv, which is the transformation of the convex cref via the geomet...
Definition: bgeot_geotrans_inv.h:242
bgeot::geotrans_inv
handles the geometric inversion for a given (supposedly quite large) set of points
Definition: bgeot_geotrans_inv.h:178
bgeot_kdtree.h
Simple implementation of a KD-tree.
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:64
bgeot::kdtree_tab_type
std::vector< index_node_pair > kdtree_tab_type
store a set of points with associated indexes.
Definition: bgeot_kdtree.h:68
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
bgeot::kdtree::add_point
size_type add_point(const base_node &n)
insert a new point
Definition: bgeot_kdtree.h:116
bgeot::convex
generic definition of a convex ( bgeot::convex_structure + vertices coordinates )
Definition: bgeot_convex.h:50
bgeot::geotrans_inv::points_in_box
size_type points_in_box(kdtree_tab_type &ipts, const base_node &min, const base_node &max) const
Find all the points present in the box between min and max.
Definition: bgeot_geotrans_inv.h:201
bgeot::kdtree
Balanced tree over a set of points.
Definition: bgeot_kdtree.h:102
bgeot::kdtree::clear
void clear()
reset the tree, remove all points
Definition: bgeot_kdtree.h:113
bgeot::geotrans_inv::nb_points
size_type nb_points(void) const
Number of points.
Definition: bgeot_geotrans_inv.h:194
bgeot::small_vector< scalar_type >
bgeot_geometric_trans.h
Geometric transformations on convexes.
bgeot::geotrans_inv::add_points
void add_points(const CONT &c)
Add the points contained in c to the list of points.
Definition: bgeot_geotrans_inv.h:187
bgeot_small_vector.h
Small (dim < 8) vectors.
bgeot::geotrans_inv::add_point
size_type add_point(base_node p)
Add point p to the list of points.
Definition: bgeot_geotrans_inv.h:196
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
bgeot::kdtree::add_point_with_id
void add_point_with_id(const base_node &n, size_type i)
insert a new point, with an associated number.
Definition: bgeot_kdtree.h:120