GetFEM  5.4.2
getfem_interpolation_on_deformed_domains.cc
1 /*===========================================================================
2 
3  Copyright (C) 2013-2020 Yves Renard, Konstantinos Poulios and Andriy Andreykiv.
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 
23 #include "getfem/getfem_models.h"
24 
25 namespace getfem {
26 
27 // Structure describing a contact boundary (or contact body)
28 struct contact_boundary {
29  size_type region; // boundary region for the slave (source)
30  // and volume region for the master (target)
31  const getfem::mesh_fem *mfu; // F.e.m. for the displacement.
32  std::string dispname; // Variable name for the displacement
33  mutable const model_real_plain_vector *U;// Displacement
34  mutable model_real_plain_vector U_unred; // Unreduced displacement
35 
36  contact_boundary(size_type r, const mesh_fem *mf, const std::string &dn)
37  : region(r), mfu(mf), dispname(dn)
38  {}
39 };
40 
41 //extract element displacements from a contact boundary object
42 base_small_vector element_U(const contact_boundary &cb, size_type cv)
43 {
44  auto U_elm = base_small_vector{};
45  slice_vector_on_basic_dof_of_element(*(cb.mfu), *cb.U, cv, U_elm);
46  return U_elm;
47 }
48 
49 //Returns an iterator of a box which centre is closest to the given point
50 auto most_central_box(const bgeot::rtree::pbox_set &bset,
51  const bgeot::base_node &pt) -> decltype(begin(bset))
52 {
53  using namespace std;
54 
55  auto itmax = begin(bset);
56 
57  auto it = itmax;
58  if (bset.size() > 1) {
59  auto rate_max = scalar_type{-1};
60  for (; it != end(bset); ++it) {
61  auto rate_box = scalar_type{1};
62  for (size_type i = 0; i < pt.size(); ++i) {
63  auto h = (*it)->max->at(i) - (*it)->min->at(i);
64  if (h > 0.) {
65  auto rate = min((*it)->max->at(i) - pt[i], pt[i] - (*it)->min->at(i)) / h;
66  rate_box = min(rate, rate_box);
67  }
68  }
69  if (rate_box > rate_max) {
70  itmax = it;
71  rate_max = rate_box;
72  }
73  }
74  }
75 
76  return itmax;
77 }
78 
79 //Transformation that creates identity mapping between two contact boundaries,
80 //deformed with provided displacement fields
81 class interpolate_transformation_on_deformed_domains
82  : public virtual_interpolate_transformation {
83 
84  contact_boundary master;//also marked with a target or Y prefix/suffix
85  contact_boundary slave; //also marked with a source or X prefix/suffix
86 
87  mutable bgeot::rtree element_boxes;
88  mutable std::map<size_type, std::vector<size_type>> box_to_convex; //index to obtain a convex
89  //number from a box number
90  mutable bgeot::geotrans_inv_convex gic;
91  mutable fem_precomp_pool fppool;
92 
93  //Create a box tree based on the deformed elements of the master (target)
94  void compute_element_boxes() const { // called by init
95  base_matrix G;
96  model_real_plain_vector Uelm; //element displacement
97  element_boxes.clear();
98 
99  auto bnum = master.region;
100  auto &mfu = *(master.mfu);
101  auto &U = *(master.U);
102  auto &m = mfu.linked_mesh();
103  auto N = m.dim();
104 
105  base_node Xdeformed(N), bmin(N), bmax(N);
106  auto region = m.region(bnum);
107 
108  //the box tree creation and subsequent transformation inversion
109  //should be done for all elements of the master, while integration
110  //will be performed only on a thread partition of the slave
111  region.prohibit_partitioning();
112 
113  GMM_ASSERT1(mfu.get_qdim() == N, "Wrong mesh_fem qdim");
114 
115  dal::bit_vector points_already_interpolated;
116  std::vector<base_node> transformed_points(m.nb_max_points());
117  box_to_convex.clear();
118 
119  for (getfem::mr_visitor v(region, m); !v.finished(); ++v) {
120  auto cv = v.cv();
121  auto pgt = m.trans_of_convex(cv);
122  auto pf_s = mfu.fem_of_element(cv);
123  auto pfp = fppool(pf_s, pgt->pgeometric_nodes());
124 
125  slice_vector_on_basic_dof_of_element(mfu, U, cv, Uelm);
126  mfu.linked_mesh().points_of_convex(cv, G);
127 
128  auto ctx = fem_interpolation_context{pgt, pfp, size_type(-1), G, cv};
129  auto nb_pt = pgt->structure()->nb_points();
130 
131  for (size_type k = 0; k < nb_pt; ++k) {
132  auto ind = m.ind_points_of_convex(cv)[k];
133 
134  // computation of a transformed vertex
135  ctx.set_ii(k);
136  if (points_already_interpolated.is_in(ind)) {
137  Xdeformed = transformed_points[ind];
138  } else {
139  pf_s->interpolation(ctx, Uelm, Xdeformed, dim_type{N});
140  Xdeformed += ctx.xreal(); //Xdeformed = U + Xo
141  transformed_points[ind] = Xdeformed;
142  points_already_interpolated.add(ind);
143  }
144 
145  if (k == 0) // computation of bounding box
146  bmin = bmax = Xdeformed;
147  else {
148  for (size_type l = 0; l < N; ++l) {
149  bmin[l] = std::min(bmin[l], Xdeformed[l]);
150  bmax[l] = std::max(bmax[l], Xdeformed[l]);
151  }
152  }
153  }
154 
155  // Store the bounding box and additional information.
156  box_to_convex[element_boxes.add_box(bmin, bmax)].push_back(cv);
157  }
158  element_boxes.build_tree();
159  }
160 
161  fem_interpolation_context deformed_master_context(size_type cv) const
162  {
163  auto &mfu = *(master.mfu);
164  auto G = base_matrix{};
165  auto pfu = mfu.fem_of_element(cv);
166  auto pgt = master.mfu->linked_mesh().trans_of_convex(cv);
167  auto pfp = fppool(pfu, pgt->pgeometric_nodes());
168  master.mfu->linked_mesh().points_of_convex(cv, G);
169  return {pgt, pfp, size_type(-1), G, cv};
170  }
171 
172  std::vector<bgeot::base_node> deformed_master_nodes(size_type cv) const {
173  using namespace bgeot;
174  using namespace std;
175 
176  auto nodes = vector<base_node>{};
177 
178  auto U_elm = element_U(master, cv);
179  auto &mfu = *(master.mfu);
180  auto G = base_matrix{};
181  auto pfu = mfu.fem_of_element(cv);
182  auto pgt = master.mfu->linked_mesh().trans_of_convex(cv);
183  auto pfp = fppool(pfu, pgt->pgeometric_nodes());
184  auto N = mfu.linked_mesh().dim();
185  auto pt = base_node(N);
186  auto U = base_small_vector(N);
187  master.mfu->linked_mesh().points_of_convex(cv, G);
188  auto ctx = fem_interpolation_context{pgt, pfp, size_type(-1), G, cv};
189  auto nb_pt = pgt->structure()->nb_points();
190  nodes.reserve(nb_pt);
191  for (size_type k = 0; k < nb_pt; ++k) {
192  ctx.set_ii(k);
193  pfu->interpolation(ctx, U_elm, U, dim_type{N});
194  gmm::add(ctx.xreal(), U, pt);
195  nodes.push_back(pt);
196  }
197 
198  return nodes;
199  }
200 
201 public:
202 
203  interpolate_transformation_on_deformed_domains(
204  size_type source_region,
205  const getfem::mesh_fem &mf_source,
206  const std::string &source_displacements,
207  size_type target_region,
208  const getfem::mesh_fem &mf_target,
209  const std::string &target_displacements)
210  :
211  slave{source_region, &mf_source, source_displacements},
212  master{target_region, &mf_target, target_displacements}
213 {}
214 
215 
216  void extract_variables(const ga_workspace &workspace,
217  std::set<var_trans_pair> &vars,
218  bool ignore_data,
219  const mesh &m_x,
220  const std::string &interpolate_name) const override {
221  if (!ignore_data || !(workspace.is_constant(master.dispname))){
222  vars.emplace(master.dispname, interpolate_name);
223  vars.emplace(slave.dispname, "");
224  }
225  }
226 
227  void init(const ga_workspace &workspace) const override {
228 
229  for (auto pcb : std::list<const contact_boundary*>{&master, &slave}) {
230  auto &mfu = *(pcb->mfu);
231  if (mfu.is_reduced()) {
232  gmm::resize(pcb->U_unred, mfu.nb_basic_dof());
233  mfu.extend_vector(workspace.value(pcb->dispname), pcb->U_unred);
234  pcb->U = &(pcb->U_unred);
235  } else {
236  pcb->U = &(workspace.value(pcb->dispname));
237  }
238  }
239  compute_element_boxes();
240  };
241 
242  void finalize() const override {
243  element_boxes.clear();
244  box_to_convex.clear();
245  master.U_unred.clear();
246  slave.U_unred.clear();
247  fppool.clear();
248  }
249 
250  int transform(const ga_workspace &workspace,
251  const mesh &m_x,
252  fem_interpolation_context &ctx_x,
253  const base_small_vector &/*Normal*/,
254  const mesh **m_t,
255  size_type &cv,
256  short_type &face_num,
257  base_node &P_ref,
258  base_small_vector &N_y,
259  std::map<var_trans_pair, base_tensor> &derivatives,
260  bool compute_derivatives) const override {
261 
262  auto &target_mesh = master.mfu->linked_mesh();
263  *m_t = &target_mesh;
264  auto transformation_success = false;
265 
266  using namespace gmm;
267  using namespace bgeot;
268  using namespace std;
269 
270  //compute a deformed point of the slave
271  auto cv_x = ctx_x.convex_num();
272  auto U_elm_x = element_U(slave, cv_x);
273  auto &mfu_x = *(slave.mfu);
274  auto pfu_x = mfu_x.fem_of_element(cv_x);
275  auto N = mfu_x.linked_mesh().dim();
276  auto U_x = base_small_vector(N);
277  auto G_x = base_matrix{}; //coordinates of the source element nodes
278  m_x.points_of_convex(cv_x, G_x);
279  ctx_x.set_pf(pfu_x);
280  pfu_x->interpolation(ctx_x, U_elm_x, U_x, dim_type{N});
281  auto pt_x = base_small_vector(N); //deformed point of the slave
282  add(ctx_x.xreal(), U_x, pt_x);
283 
284  //Find the best box from the master (target) that
285  //corresponds to this point (The box which centre is the closest to the point).
286  //Obtain the corresponding element number using the box id and box_to_convex
287  //indices. Compute deformed nodes of the target element. Invert the geometric
288  //transformation of the target element with deformed nodes, obtaining this way
289  //reference coordinates of the target element
290  auto bset = rtree::pbox_set{};
291  element_boxes.find_boxes_at_point(pt_x, bset);
292  while (!bset.empty())
293  {
294  auto itmax = most_central_box(bset, pt_x);
295 
296  for (auto i : box_to_convex.at((*itmax)->id))
297  {
298  auto deformed_nodes_y = deformed_master_nodes(i);
299  gic.init(deformed_nodes_y, target_mesh.trans_of_convex(i));
300  auto converged = true;
301  auto is_in = gic.invert(pt_x, P_ref, converged);
302  if (is_in && converged) {
303  cv = i;
304  face_num = static_cast<short_type>(-1);
305  transformation_success = true;
306  break;
307  }
308  }
309  if (transformation_success || (bset.size() == 1)) break;
310  bset.erase(itmax);
311  }
312 
313  //Since this transformation can be seen as Xsource + Usource - Utarget,
314  //the corresponding stiffnesses are identity matrix for Usource and
315  //minus identity for Utarget. The required answer in this function is
316  //stiffness X shape function. Hence, returning shape function for Usource
317  //and min shape function for Utarget
318  if (compute_derivatives && transformation_success) {
319  GMM_ASSERT2(derivatives.size() == 2,
320  "Expecting to return derivatives only for Umaster and Uslave");
321 
322  for (auto &pair : derivatives)
323  {
324  if (pair.first.varname == slave.dispname)
325  {
326  auto base_ux = base_tensor{};
327  auto vbase_ux = base_matrix{} ;
328  ctx_x.base_value(base_ux);
329  auto qdim_ux = pfu_x->target_dim();
330  auto ndof_ux = pfu_x->nb_dof(cv_x) * N / qdim_ux;
331  vectorize_base_tensor(base_ux, vbase_ux, ndof_ux, qdim_ux, N);
332  pair.second.adjust_sizes(ndof_ux, N);
333  copy(vbase_ux.as_vector(), pair.second.as_vector());
334  }
335  else
336  if (pair.first.varname == master.dispname)
337  {
338  auto ctx_y = deformed_master_context(cv);
339  ctx_y.set_xref(P_ref);
340  auto base_uy = base_tensor{};
341  auto vbase_uy = base_matrix{} ;
342  ctx_y.base_value(base_uy);
343  auto pfu_y = master.mfu->fem_of_element(cv);
344  auto dim_y = master.mfu->linked_mesh().dim();
345  auto qdim_uy = pfu_y->target_dim();
346  auto ndof_uy = pfu_y->nb_dof(cv) * dim_y / qdim_uy;
347  vectorize_base_tensor(base_uy, vbase_uy, ndof_uy, qdim_uy, dim_y);
348  pair.second.adjust_sizes(ndof_uy, dim_y);
349  copy(vbase_uy.as_vector(), pair.second.as_vector());
350  scale(pair.second.as_vector(), -1.);
351  }
352  else GMM_ASSERT2(false, "unexpected derivative variable");
353  }
354  }
355 
356  return transformation_success ? 1 : 0;
357  }
358 
359 };
360 
362  (ga_workspace &workspace, const std::string &transname,
363  const mesh &source_mesh, const std::string &source_displacements,
364  const mesh_region &source_region, const mesh &target_mesh,
365  const std::string &target_displacements, const mesh_region &target_region)
366  {
367  auto pmf_source = workspace.associated_mf(source_displacements);
368  auto pmf_target = workspace.associated_mf(target_displacements);
369  auto p_transformation
370  = std::make_shared<interpolate_transformation_on_deformed_domains>(source_region.id(),
371  *pmf_source,
372  source_displacements,
373  target_region.id(),
374  *pmf_target,
375  target_displacements);
376  workspace.add_interpolate_transformation(transname, p_transformation);
377  }
378 
380  (model &md, const std::string &transname,
381  const mesh &source_mesh, const std::string &source_displacements,
382  const mesh_region &source_region, const mesh &target_mesh,
383  const std::string &target_displacements, const mesh_region &target_region)
384  {
385  auto &mf_source = md.mesh_fem_of_variable(source_displacements);
386  auto mf_target = md.mesh_fem_of_variable(target_displacements);
387  auto p_transformation
388  = std::make_shared<interpolate_transformation_on_deformed_domains>(source_region.id(),
389  mf_source,
390  source_displacements,
391  target_region.id(),
392  mf_target,
393  target_displacements);
394  md.add_interpolate_transformation(transname, p_transformation);
395  }
396 
397 } /* end of namespace getfem. */
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:64
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::add_interpolate_transformation_on_deformed_domains
void add_interpolate_transformation_on_deformed_domains(ga_workspace &workspace, const std::string &transname, const mesh &source_mesh, const std::string &source_displacements, const mesh_region &source_region, const mesh &target_mesh, const std::string &target_displacements, const mesh_region &target_region)
Add a transformation to the workspace that creates an identity mapping between two meshes in deformed...
Definition: getfem_interpolation_on_deformed_domains.cc:362
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::rtree
Balanced tree of n-dimensional rectangles.
Definition: bgeot_rtree.h:98
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem_models.h
Model representation in Getfem.
getfem_generic_assembly.h
A language for generic assembly of pde boundary value problems.
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
getfem::model::add_interpolate_transformation
void add_interpolate_transformation(const std::string &name, pinterpolate_transformation ptrans)
Add an interpolate transformation to the model to be used with the generic assembly.
Definition: getfem_models.h:1106
getfem::model::mesh_fem_of_variable
const mesh_fem & mesh_fem_of_variable(const std::string &name) const
Gives the access to the mesh_fem of a variable if any.
Definition: getfem_models.cc:2891
bgeot::small_vector< scalar_type >
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:27
getfem::slice_vector_on_basic_dof_of_element
void slice_vector_on_basic_dof_of_element(const mesh_fem &mf, const VEC1 &vec, size_type cv, VEC2 &coeff, size_type qmult1=size_type(-1), size_type qmult2=size_type(-1))
Given a mesh_fem.
Definition: getfem_mesh_fem.h:659
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:977
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
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