GetFEM  5.4.2
getfem_interpolation.cc
1 /*===========================================================================
2 
3  Copyright (C) 2001-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 
22 
24 
25 namespace getfem {
26 
27  size_type mesh_trans_inv::id_of_point(size_type ipt) const {
28 
29  if (!ids.empty()) {
30  map_iterator it=ids.find(ipt);
31  if (it != ids.end())
32  return it->second;
33  }
34  // otherwise assume that the point id is the point index
35  return ipt;
36  }
37 
38  void mesh_trans_inv::points_on_convex(size_type cv,
39  std::vector<size_type> &itab) const {
40  itab.resize(pts_cvx[cv].size()); size_type j = 0;
41  for (set_iterator it = pts_cvx[cv].begin(); it != pts_cvx[cv].end(); ++it)
42  itab[j++] = *it;
43  }
44 
45  size_type mesh_trans_inv::point_on_convex(size_type cv, size_type i) const {
46  set_iterator it = pts_cvx[cv].begin();
47  for (size_type j = 0; it != pts_cvx[cv].end() && j < i; ++it, ++j) {}
48  GMM_ASSERT1(it != pts_cvx[cv].end(), "internal error");
49  return *it;
50  }
51 
52  void mesh_trans_inv::distribute(int extrapolation, mesh_region rg_source) {
53 
54  rg_source.from_mesh(msh);
55  rg_source.error_if_not_convexes();
56  bool all_convexes = (rg_source.id() == mesh_region::all_convexes().id());
57 
58  size_type nbpts = nb_points();
59  size_type nbcvx = msh.nb_allocated_convex();
60  ref_coords.resize(nbpts);
61  std::vector<double> dist(nbpts);
62  std::vector<size_type> cvx_pts(nbpts);
63  pts_cvx.clear(); pts_cvx.resize(nbcvx);
64  base_node min, max, pt_ref; /* bound of the box enclosing the convex */
66  dal::bit_vector npt, cv_on_bound;
67  npt.add(0, nbpts);
68  scalar_type mult = scalar_type(1);
69 
70  bool projection_into_element(extrapolation == 0);
71 
72  do {
73  for (dal::bv_visitor j(rg_source.index()); !j.finished(); ++j) {
74  if (mult > scalar_type(1) && !(cv_on_bound.is_in(j))) continue;
75  bgeot::pgeometric_trans pgt = msh.trans_of_convex(j);
76  bounding_box(min, max, msh.points_of_convex(j), pgt);
77  for (size_type k=0; k < min.size(); ++k) { min[k]-=EPS; max[k]+=EPS; }
78  if (extrapolation == 2) {
79  if (mult == scalar_type(1))
80  for (short_type f = 0; f < msh.nb_faces_of_convex(j); ++f) {
81  size_type neighbor_cv = msh.neighbor_of_convex(j, f);
82  if (!all_convexes && neighbor_cv != size_type(-1)) {
83  // check if the neighbor is also contained in rg_source ...
84  if (!rg_source.is_in(neighbor_cv))
85  cv_on_bound.add(j); // ... if not, treat the element as a boundary one
86  }
87  else // boundary element of the overall mesh
88  cv_on_bound.add(j);
89  }
90  if (cv_on_bound.is_in(j)) {
91  scalar_type h = scalar_type(0);
92  for (size_type k=0; k < min.size(); ++k)
93  h = std::max(h, max[k] - min[k]);
94  for (size_type k=0; k < min.size(); ++k)
95  { min[k]-=mult*h; max[k]+=mult*h; }
96  }
97  }
98  points_in_box(boxpts, min, max);
99 
100  if (boxpts.size() > 0) gic.init(msh.points_of_convex(j), pgt);
101 
102  for (size_type l = 0; l < boxpts.size(); ++l) {
103  size_type ind = boxpts[l].i;
104  if (npt[ind] || dist[ind] > 0) {
105  bool converged;
106  bool gicisin = gic.invert(boxpts[l].n, pt_ref, converged, EPS,
107  projection_into_element);
108  bool toadd = extrapolation || gicisin;
109  double isin = pgt->convex_ref()->is_in(pt_ref);
110 
111  if (toadd && !(npt[ind])) {
112  if (isin < dist[ind]) pts_cvx[cvx_pts[ind]].erase(ind);
113  else toadd = false;
114  }
115  if (toadd) {
116 // if (mult > 1.5) {
117 // cout << "adding " << ind << "ref_coord = " << pt_ref
118 // << " cv = " << j << " mult = " << mult << endl;
119 // }
120  ref_coords[ind] = pt_ref;
121  dist[ind] = isin; cvx_pts[ind] = j;
122  pts_cvx[j].insert(ind);
123  npt.sup(ind);
124  }
125  }
126  }
127  }
128  mult *= scalar_type(2);
129  } while (npt.card() > 0 && extrapolation == 2);
130  }
131 } /* end of namespace getfem. */
132 
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
getfem::mesh_region::all_convexes
static mesh_region all_convexes()
provide a default value for the mesh_region parameters of assembly procedures etc.
Definition: getfem_mesh_region.h:142
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem_interpolation.h
Interpolation of fields from a mesh_fem onto another.
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186