GetFEM  5.4.3
getfem_contact_and_friction_common.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2011-2020 Yves Renard, Konstantinos Poulios.
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 getfem_contact_and_friction_common.h
33  @author Yves Renard <[email protected]>
34  @author Konstantinos Poulios <[email protected]>
35  @date November, 2011.
36  @brief Comomon tools for unilateral contact and Coulomb friction bricks.
37  */
38 #ifndef GETFEM_CONTACT_AND_FRICTION_COMMON_H__
39 #define GETFEM_CONTACT_AND_FRICTION_COMMON_H__
40 
41 #include "getfem_models.h"
42 #include "getfem_assembling_tensors.h"
43 #include "getfem/bgeot_rtree.h"
44 #include <getfem/getfem_mesher.h>
45 
46 
47 namespace getfem {
48 
49  //=========================================================================
50  //
51  // Projection on a ball and gradient of the projection.
52  //
53  //=========================================================================
54 
55  template<typename VEC> void ball_projection(const VEC &x,
56  scalar_type radius) {
57  if (radius <= scalar_type(0))
58  gmm::clear(const_cast<VEC&>(x));
59  else {
60  scalar_type a = gmm::vect_norm2(x);
61  if (a > radius) gmm::scale(const_cast<VEC&>(x), radius/a);
62  }
63  }
64 
65  template<typename VEC, typename VECR>
66  void ball_projection_grad_r(const VEC &x, scalar_type radius,
67  VECR &g) {
68  if (radius > scalar_type(0)) {
69  scalar_type a = gmm::vect_norm2(x);
70  if (a >= radius) {
71  gmm::copy(x, g); gmm::scale(g, scalar_type(1)/a);
72  return;
73  }
74  }
75  gmm::clear(g);
76  }
77 
78  template <typename VEC, typename MAT>
79  void ball_projection_grad(const VEC &x, scalar_type radius, MAT &g) {
80  if (radius <= scalar_type(0)) { gmm::clear(g); return; }
81  gmm::copy(gmm::identity_matrix(), g);
82  scalar_type a = gmm::vect_norm2(x);
83  if (a >= radius) {
84  gmm::scale(g, radius/a);
85  // gmm::rank_one_update(g, gmm::scaled(x, -radius/(a*a*a)), x);
86  for (size_type i = 0; i < x.size(); ++i)
87  for (size_type j = 0; j < x.size(); ++j)
88  g(i,j) -= radius*x[i]*x[j] / (a*a*a);
89  }
90  }
91 
92  template <typename VEC, typename VECR>
93  void coupled_projection(const VEC &x, const VEC &n,
94  scalar_type f, VECR &g) {
95  scalar_type xn = gmm::vect_sp(x, n);
96  scalar_type xnm = gmm::neg(xn);
97  scalar_type th = f * xnm;
98  scalar_type xtn = gmm::sqrt(gmm::vect_norm2_sqr(x) - xn*xn);
99 
100  gmm::copy(gmm::scaled(n, -xnm), g);
101  if (th > scalar_type(0)) {
102  if (xtn <= th) {
103  gmm::add(x, g);
104  gmm::add(gmm::scaled(n, -xn), g);
105  } else {
106  gmm::add(gmm::scaled(x, f*xnm/xtn), g);
107  gmm::add(gmm::scaled(n, -f*xnm*xn/xtn), g);
108  }
109  }
110  }
111 
112 
113  template <typename VEC, typename MAT>
114  void coupled_projection_grad(const VEC &x, const VEC &n,
115  scalar_type f, MAT &g) {
116  scalar_type xn = gmm::vect_sp(x, n);
117  scalar_type xnm = gmm::neg(xn);
118  scalar_type th = f * xnm;
119  scalar_type xtn = gmm::sqrt(gmm::vect_norm2_sqr(x) - xn*xn);
120  size_type N = gmm::vect_size(x);
121  gmm::clear(g);
122 
123  if (th > scalar_type(0)) {
124  if (xtn <= th) {
125  gmm::copy(gmm::identity_matrix(), g);
126  gmm::rank_one_update(g, gmm::scaled(n, -scalar_type(1)), n);
127  } else if (xn < scalar_type(0)) {
128  static base_small_vector t; gmm::resize(t, N);
129  gmm::add(x, gmm::scaled(n, -xn), t);
130  gmm::scale(t, scalar_type(1)/xtn);
131  if (N > 2) {
132  gmm::copy(gmm::identity_matrix(), g);
133  gmm::rank_one_update(g, gmm::scaled(t, -scalar_type(1)), t);
134  gmm::rank_one_update(g, gmm::scaled(n, -scalar_type(1)), n);
135  gmm::scale(g, -xn*th/xtn);
136  }
137  gmm::rank_one_update(g, gmm::scaled(t, -f), n);
138  }
139  }
140 
141  if (xn < scalar_type(0)) gmm::rank_one_update(g, n, n);
142  }
143 
144  //=========================================================================
145  //
146  // De Saxce projection and its gradients.
147  //
148  //=========================================================================
149 
150 
151  template<typename VEC>
152  void De_Saxce_projection(const VEC &x, const VEC &n_, scalar_type f) {
153  static base_small_vector n; // For more robustness, n_ is not supposed unitary
154  size_type N = gmm::vect_size(x);
155  gmm::resize(n, N);
156  gmm::copy(gmm::scaled(n_, scalar_type(1)/gmm::vect_norm2(n_)), n);
157  scalar_type xn = gmm::vect_sp(x, n);
158  scalar_type nxt = sqrt(gmm::abs(gmm::vect_norm2_sqr(x) - xn*xn));
159  if (xn >= scalar_type(0) && f * nxt <= xn) {
160  gmm::clear(const_cast<VEC&>(x));
161  } else if (xn > scalar_type(0) || nxt > -f*xn) {
162  gmm::add(gmm::scaled(n, -xn), const_cast<VEC&>(x));
163  gmm::scale(const_cast<VEC&>(x), -f / nxt);
164  gmm::add(n, const_cast<VEC&>(x));
165  gmm::scale(const_cast<VEC&>(x), (xn - f * nxt) / (f*f+scalar_type(1)));
166  }
167  }
168 
169  template<typename VEC, typename MAT>
170  void De_Saxce_projection_grad(const VEC &x, const VEC &n_,
171  scalar_type f, MAT &g) {
172  static base_small_vector n;
173  size_type N = gmm::vect_size(x);
174  gmm::resize(n, N);
175  gmm::copy(gmm::scaled(n_, scalar_type(1)/gmm::vect_norm2(n_)), n);
176  scalar_type xn = gmm::vect_sp(x, n);
177  scalar_type nxt = sqrt(gmm::abs(gmm::vect_norm2_sqr(x) - xn*xn));
178 
179 
180  if (xn > scalar_type(0) && f * nxt <= xn) {
181  gmm::clear(g);
182  } else if (xn > scalar_type(0) || nxt > -f*xn) {
183  static base_small_vector xt;
184  gmm::resize(xt, N);
185  gmm::add(x, gmm::scaled(n, -xn), xt);
186  gmm::scale(xt, scalar_type(1)/nxt);
187 
188  if (N > 2) {
189  gmm::copy(gmm::identity_matrix(), g);
190  gmm::rank_one_update(g, gmm::scaled(n, -scalar_type(1)), n);
191  gmm::rank_one_update(g, gmm::scaled(xt, -scalar_type(1)), xt);
192  gmm::scale(g, f*(f - xn/nxt));
193  } else {
194  gmm::clear(g);
195  }
196 
197  gmm::scale(xt, -f); gmm::add(n, xt);
198  gmm::rank_one_update(g, xt, xt);
199  gmm::scale(g, scalar_type(1) / (f*f+scalar_type(1)));
200  } else {
201  gmm::copy(gmm::identity_matrix(), g);
202  }
203  }
204 
205 
206  template<typename VEC, typename MAT>
207  static void De_Saxce_projection_gradn(const VEC &x, const VEC &n_,
208  scalar_type f, MAT &g) {
209  static base_small_vector n;
210  size_type N = gmm::vect_size(x);
211  scalar_type nn = gmm::vect_norm2(n_);
212  gmm::resize(n, N);
213  gmm::copy(gmm::scaled(n_, scalar_type(1)/nn), n);
214  scalar_type xn = gmm::vect_sp(x, n);
215  scalar_type nxt = sqrt(gmm::abs(gmm::vect_norm2_sqr(x) - xn*xn));
216  gmm::clear(g);
217 
218  if (!(xn > scalar_type(0) && f * nxt <= xn)
219  && (xn > scalar_type(0) || nxt > -f*xn)) {
220  static base_small_vector xt, aux;
221  gmm::resize(xt, N); gmm::resize(aux, N);
222  gmm::add(x, gmm::scaled(n, -xn), xt);
223  gmm::scale(xt, scalar_type(1)/nxt);
224 
225  scalar_type c = (scalar_type(1) + f*xn/nxt)/nn;
226  for (size_type i = 0; i < N; ++i) g(i,i) = c;
227  gmm::rank_one_update(g, gmm::scaled(n, -c), n);
228  gmm::rank_one_update(g, gmm::scaled(n, f/nn), xt);
229  gmm::rank_one_update(g, gmm::scaled(xt, -f*xn/(nn*nxt)), xt);
230  gmm::scale(g, xn - f*nxt);
231 
232  gmm::add(gmm::scaled(xt, -f), n, aux);
233  gmm::rank_one_update(g, aux, gmm::scaled(xt, (nxt+f*xn)/nn));
234 
235  gmm::scale(g, scalar_type(1) / (f*f+scalar_type(1)));
236  }
237  }
238 
239  //=========================================================================
240  //
241  // Some basic assembly functions.
242  //
243  //=========================================================================
244 
245  template <typename MAT1, typename MAT2>
246  void mat_elem_assembly(const MAT1 &M_, const MAT2 &Melem,
247  const mesh_fem &mf1, size_type cv1,
248  const mesh_fem &mf2, size_type cv2) {
249  MAT1 &M = const_cast<MAT1 &>(M_);
250  typedef typename gmm::linalg_traits<MAT1>::value_type T;
251  T val;
252  mesh_fem::ind_dof_ct cvdof1 = mf1.ind_basic_dof_of_element(cv1);
253  mesh_fem::ind_dof_ct cvdof2 = mf2.ind_basic_dof_of_element(cv2);
254 
255  GMM_ASSERT1(cvdof1.size() == gmm::mat_nrows(Melem)
256  && cvdof2.size() == gmm::mat_ncols(Melem),
257  "Dimensions mismatch");
258 
259  if (mf1.is_reduced()) {
260  if (mf2.is_reduced()) {
261  for (size_type i = 0; i < cvdof1.size(); ++i)
262  for (size_type j = 0; j < cvdof2.size(); ++j)
263  if ((val = Melem(i,j)) != T(0))
264  asmrankoneupdate
265  (M, gmm::mat_row(mf1.extension_matrix(), cvdof1[i]),
266  gmm::mat_row(mf2.extension_matrix(), cvdof2[j]), val);
267  } else {
268  for (size_type i = 0; i < cvdof1.size(); ++i)
269  for (size_type j = 0; j < cvdof2.size(); ++j)
270  if ((val = Melem(i,j)) != T(0))
271  asmrankoneupdate
272  (M, gmm::mat_row(mf1.extension_matrix(), cvdof1[i]),
273  cvdof2[j], val);
274  }
275  } else {
276  if (mf2.is_reduced()) {
277  for (size_type i = 0; i < cvdof1.size(); ++i)
278  for (size_type j = 0; j < cvdof2.size(); ++j)
279  if ((val = Melem(i,j)) != T(0))
280  asmrankoneupdate
281  (M, cvdof1[i],
282  gmm::mat_row(mf2.extension_matrix(), cvdof2[j]), val);
283  } else {
284  for (size_type i = 0; i < cvdof1.size(); ++i)
285  for (size_type j = 0; j < cvdof2.size(); ++j)
286  if ((val = Melem(i,j)) != T(0))
287  M(cvdof1[i], cvdof2[j]) += val;
288  }
289  }
290  }
291 
292 
293  template <typename VEC1, typename VEC2>
294  void vec_elem_assembly(const VEC1 &V_, const VEC2 &Velem,
295  const mesh_fem &mf, size_type cv) {
296  VEC1 &V = const_cast<VEC1 &>(V_);
297  typedef typename gmm::linalg_traits<VEC1>::value_type T;
298  std::vector<size_type> cvdof(mf.ind_basic_dof_of_element(cv).begin(),
299  mf.ind_basic_dof_of_element(cv).end());
300 
301  GMM_ASSERT1(cvdof.size() == gmm::vect_size(Velem), "Dimensions mismatch");
302 
303  if (mf.is_reduced()) {
304  T val;
305  for (size_type i = 0; i < cvdof.size(); ++i)
306  if ((val = Velem[i]) != T(0))
307  gmm::add(gmm::scaled(gmm::mat_row(mf.extension_matrix(), cvdof[i]),
308  val), V);
309  } else {
310  for (size_type i = 0; i < cvdof.size(); ++i) V[cvdof[i]] += Velem[i];
311  }
312  }
313 
314  template <typename MAT1, typename MAT2>
315  void mat_elem_assembly(const MAT1 &M_, const gmm::sub_interval &I1,
316  const gmm::sub_interval &I2,
317  const MAT2 &Melem,
318  const mesh_fem &mf1, size_type cv1,
319  const mesh_fem &mf2, size_type cv2) {
320  MAT1 &M = const_cast<MAT1 &>(M_);
321  typedef typename gmm::linalg_traits<MAT1>::value_type T;
322  T val;
323 
324  mesh_fem::ind_dof_ct cvdof1 = mf1.ind_basic_dof_of_element(cv1);
325  mesh_fem::ind_dof_ct cvdof2 = mf2.ind_basic_dof_of_element(cv2);
326 
327  GMM_ASSERT1(cvdof1.size() == gmm::mat_nrows(Melem)
328  && cvdof2.size() == gmm::mat_ncols(Melem),
329  "Dimensions mismatch");
330 
331  if (mf1.is_reduced()) {
332  if (mf2.is_reduced()) {
333  for (size_type i = 0; i < cvdof1.size(); ++i)
334  for (size_type j = 0; j < cvdof2.size(); ++j)
335  if ((val = Melem(i,j)) != T(0))
336  asmrankoneupdate
337  (gmm::sub_matrix(M, I1, I2),
338  gmm::mat_row(mf1.extension_matrix(), cvdof1[i]),
339  gmm::mat_row(mf2.extension_matrix(), cvdof2[j]), val);
340  } else {
341  for (size_type i = 0; i < cvdof1.size(); ++i)
342  for (size_type j = 0; j < cvdof2.size(); ++j)
343  if ((val = Melem(i,j)) != T(0))
344  asmrankoneupdate
345  (gmm::sub_matrix(M, I1, I2),
346  gmm::mat_row(mf1.extension_matrix(), cvdof1[i]),
347  cvdof2[j], val);
348  }
349  } else {
350  if (mf2.is_reduced()) {
351  for (size_type i = 0; i < cvdof1.size(); ++i)
352  for (size_type j = 0; j < cvdof2.size(); ++j)
353  if ((val = Melem(i,j)) != T(0))
354  asmrankoneupdate
355  (gmm::sub_matrix(M, I1, I2), cvdof1[i],
356  gmm::mat_row(mf2.extension_matrix(), cvdof2[j]), val);
357  } else {
358  for (size_type i = 0; i < cvdof1.size(); ++i)
359  for (size_type j = 0; j < cvdof2.size(); ++j)
360  if ((val = Melem(i,j)) != T(0))
361  M(cvdof1[i]+I1.first(), cvdof2[j]+I2.first()) += val;
362  }
363  }
364  }
365 
366  template <typename VEC1, typename VEC2>
367  void vec_elem_assembly(const VEC1 &V_, const gmm::sub_interval &I,
368  const VEC2 &Velem, const mesh_fem &mf, size_type cv) {
369  VEC1 &V = const_cast<VEC1 &>(V_);
370  typedef typename gmm::linalg_traits<VEC1>::value_type T;
371  std::vector<size_type> cvdof(mf.ind_basic_dof_of_element(cv).begin(),
372  mf.ind_basic_dof_of_element(cv).end());
373 
374  GMM_ASSERT1(cvdof.size() == gmm::vect_size(Velem), "Dimensions mismatch");
375 
376  if (mf.is_reduced()) {
377  T val;
378  for (size_type i = 0; i < cvdof.size(); ++i)
379  if ((val = Velem[i]) != T(0))
380  gmm::add(gmm::scaled(gmm::mat_row(mf.extension_matrix(), cvdof[i]),
381  val), gmm::sub_vector(V, I));
382  } else {
383  for (size_type i = 0; i < cvdof.size(); ++i)
384  V[I.first()+cvdof[i]] += Velem[i];
385  }
386  }
387 
388 
389  void compute_normal(const fem_interpolation_context &ctx, size_type face,
390  bool in_reference_conf, const model_real_plain_vector &coeff,
391  base_node &n0, base_node &n, base_matrix &grad);
392 
393  //=========================================================================
394  //
395  // Structure which stores the contact boundaries, rigid obstacles and
396  // computes the contact pairs in large sliding/large deformation
397  //
398  // Deprecated
399  //
400  //=========================================================================
401 
402 
403  class multi_contact_frame {
404 
405  // Structure describing a contact boundary
406  struct contact_boundary {
407  size_type region; // Boundary number
408  const getfem::mesh_fem *mfu; // F.e.m. for the displacement.
409  const getfem::mesh_fem *mflambda; // F.e.m. for the displacement.
410  const getfem::mesh_im *mim; // Integration method for the boundary.
411  std::string multname; // Name of the optional contact stress
412  // multiplier when linked to a model.
413  size_type ind_U; // Index of displacement.
414  size_type ind_lambda; // Index of multiplier (if any).
415  bool slave;
416  contact_boundary(void) {}
417  contact_boundary(size_type r, const mesh_fem *mf,
418  const mesh_im &mi, size_type i_U, const mesh_fem *mfl,
419  size_type i_l = size_type(-1))
420  : region(r), mfu(mf), mflambda(mfl), mim(&mi),
421  ind_U(i_U), ind_lambda(i_l), slave(false) {}
422  };
423 
424 
425  size_type N; // Meshes dimensions
426  bool self_contact; // Self-contact is searched or not.
427  bool ref_conf; // Contact in reference configuration
428  // for linear elasticity small sliding contact.
429  bool use_delaunay; // Use delaunay to detect the contact pairs instead
430  // of influence boxes.
431  int nodes_mode; // 0 = Use Gauss points for both slave and master
432  // 1 = Use finite element nodes for slave and
433  // Gauss points for master.
434  // 2 = Use finite element nodes for both slave
435  // and master
436  bool raytrace; // Use raytrace instead of projection.
437 
438  scalar_type release_distance; // Limit distance beyond which the contact
439  // will not be considered. CAUTION: should be comparable to the element
440  // size (if it is too large, a too large set of influence boxes will be
441  // detected and the computation will be slow, except for delaunay option)
442 
443  scalar_type cut_angle; // Cut angle (in radian) for normal cones
444  scalar_type EPS; // Should be typically hmin/1000 (for computing
445  // gradients with finite differences
446  const model *md; // The model if the structure is linked to a model.
447 
448  typedef model_real_plain_vector VECTOR;
449  std::vector<const VECTOR *> Us; // Displacement vectors
450  std::vector<const VECTOR *> Ws; // "Velocity" vectors
451  std::vector<std::string> Unames; // Displacement vectors names.
452  std::vector<std::string> Wnames; // "Velocity" vectors names.
453  std::vector<VECTOR> ext_Us; // Unreduced displacement vectors
454  std::vector<VECTOR> ext_Ws; // Unreduced "velocity" vectors
455  std::vector<const VECTOR *> lambdas; // Displacement vectors
456  std::vector<std::string> lambdanames; // Displacement vectors names.
457  std::vector<VECTOR> ext_lambdas; // Unreduced displacement vectors
458 
459  std::vector<contact_boundary> contact_boundaries;
460 
461  std::vector<std::string> coordinates;
462  model_real_plain_vector pt, ptx, pty, ptz, ptw;
463  std::list<ga_workspace> obstacles_gw;
464  std::vector<ga_function> obstacles_f;
465  std::vector<std::string> obstacles;
466  std::vector<std::string> obstacles_velocities;
467 
468 
469  struct normal_cone : public std::vector<base_small_vector> {
470 
471  void add_normal(const base_small_vector &n)
472  { std::vector<base_small_vector>::push_back(n);}
473  normal_cone(void) {}
474  normal_cone(const base_small_vector &n)
475  : std::vector<base_small_vector>(1, n) { }
476  };
477 
478  //
479  // Influence boxes
480  //
481  struct influence_box { // Additional information for an influence box
482  size_type ind_boundary; // Boundary number
483  size_type ind_element; // Element number
484  short_type ind_face; // Face number in element
485  base_small_vector mean_normal; // Mean outward normal unit vector
486  influence_box(void) {}
487  influence_box(size_type ib, size_type ie,
488  short_type iff, const base_small_vector &n)
489  : ind_boundary(ib), ind_element(ie), ind_face(iff), mean_normal(n) {}
490  };
491 
492  bgeot::rtree element_boxes; // influence boxes
493  std::vector<influence_box> element_boxes_info;
494 
495  //
496  // Stored points (for Delaunay and slave nodal boundaries)
497  //
498 
499  struct boundary_point { // Additional information for a boundary point
500  base_node ref_point; // Point coordinate in reference configuration
501  size_type ind_boundary; // Boundary number
502  size_type ind_element; // Element number
503  short_type ind_face; // Face number in element
504  size_type ind_pt; // Dof number for fem nodes or point number
505  // of integration method (depending on nodes_mode)
506  normal_cone normals; // Set of outward unit normal vectors
507  boundary_point(void) {}
508  boundary_point(const base_node &rp, size_type ib, size_type ie,
509  short_type iff, size_type id, const base_small_vector &n)
510  : ref_point(rp), ind_boundary(ib), ind_element(ie), ind_face(iff),
511  ind_pt(id), normals(n) {}
512  };
513 
514  std::vector<base_node> boundary_points;
515  std::vector<boundary_point> boundary_points_info;
516 
517 
518  size_type add_U(const model_real_plain_vector *U, const std::string &name,
519  const model_real_plain_vector *w, const std::string &wname);
520  size_type add_lambda(const model_real_plain_vector *lambda,
521  const std::string &name);
522 
523  void extend_vectors(void);
524 
525  void normal_cone_simplification(void);
526 
527  bool test_normal_cones_compatibility(const normal_cone &nc1,
528  const normal_cone &nc2);
529 
530  bool test_normal_cones_compatibility(const base_small_vector &n,
531  const normal_cone &nc2);
532 
533  dal::bit_vector aux_dof_cv; // An auxiliary variable for are_dof_linked
534  // function (in order to be of constant complexity).
535 
536  bool are_dof_linked(size_type ib1, size_type idof1,
537  size_type ib2, size_type idof2);
538 
539  bool is_dof_linked(size_type ib1, size_type idof1,
540  size_type ib2, size_type cv);
541  public:
542 
543  struct face_info {
544  size_type ind_boundary;
545  size_type ind_element;
546  short_type ind_face;
547  face_info(void) {}
548  face_info(size_type ib, size_type ie, short_type iff)
549  : ind_boundary(ib), ind_element(ie), ind_face(iff) {}
550  };
551 
552  protected:
553 
554  std::vector<std::vector<face_info> > potential_pairs;
555 
556  void add_potential_contact_face(size_type ip, size_type ib, size_type ie,
557  short_type iff);
558  public:
559 
560  // stored information for contact pair
561  struct contact_pair {
562 
563  base_node slave_point; // The transformed slave point
564  base_small_vector slave_n; // Normal unit vector to slave surface
565  size_type slave_ind_boundary; // Boundary number
566  size_type slave_ind_element; // Element number
567  short_type slave_ind_face; // Face number in element
568  size_type slave_ind_pt; // Dof number for fem nodes or point number
569  // of integration method (depending on nodes_mode)
570 
571  base_node master_point_ref; // The master point on ref element
572  base_node master_point; // The transformed master point
573  base_small_vector master_n; // Normal unit vector to master surface
574  size_type master_ind_boundary; // Boundary number
575  size_type master_ind_element; // Element number
576  short_type master_ind_face; // Face number in element
577 
578  scalar_type signed_dist;
579 
580  size_type irigid_obstacle;
581 
582  contact_pair(void) {}
583  contact_pair(const base_node &spt, const base_small_vector &nx,
584  const boundary_point &bp,
585  const base_node &mptr, const base_node &mpt,
586  const base_small_vector &ny,
587  const face_info &mfi, scalar_type sd)
588  : slave_point(spt), slave_n(nx),
589  slave_ind_boundary(bp.ind_boundary), slave_ind_element(bp.ind_element),
590  slave_ind_face(bp.ind_face), slave_ind_pt(bp.ind_pt),
591  master_point_ref(mptr), master_point(mpt), master_n(ny),
592  master_ind_boundary(mfi.ind_boundary), master_ind_element(mfi.ind_element),
593  master_ind_face(mfi.ind_face),
594  signed_dist(sd), irigid_obstacle(size_type(-1)) {}
595  contact_pair(const base_node &spt, const base_small_vector &nx,
596  const boundary_point &bp,
597  const base_node &mpt, const base_small_vector &ny,
598  size_type ir, scalar_type sd)
599  : slave_point(spt), slave_n(nx), slave_ind_boundary(bp.ind_boundary),
600  slave_ind_element(bp.ind_element), slave_ind_face(bp.ind_face),
601  slave_ind_pt(bp.ind_pt), master_point(mpt), master_n(ny),
602  signed_dist(sd),
603  irigid_obstacle(ir) {}
604 
605  };
606 
607 
608  // Compute the influence boxes of master boundary elements. To be run
609  // before the detection of contact pairs. The influence box is the
610  // bounding box extended by a distance equal to the release distance.
611  void compute_influence_boxes(void);
612 
613  // For delaunay triangulation. Advantages compared to influence boxes:
614  // No degeneration of the algorithm complexity with refinement and
615  // more easy to extend to fictitious domain with contact.
616  // Stores all the boundary deformed points relatively to
617  // an integration method or to finite element nodes (depending on
618  // nodes_mode). Storing sufficient information to perform
619  // a Delaunay triangulation and to be able to recover the boundary
620  // number, element number, face number, unit normal vector ...
621  void compute_boundary_points(bool slave_only = false);
622  void compute_potential_contact_pairs_delaunay(void);
623  void compute_potential_contact_pairs_influence_boxes(void);
624 
625  protected:
626 
627  std::vector<contact_pair> contact_pairs;
628 
629  void clear_aux_info(void); // Delete auxiliary information
630 
631  public:
632 
633  size_type dim(void) const { return N; }
634  const std::vector<contact_pair> &ct_pairs(void) const
635  { return contact_pairs; }
636 
637 
638  const getfem::mesh_fem &mfdisp_of_boundary(size_type n) const
639  { return *(contact_boundaries[n].mfu); }
640  const getfem::mesh_fem &mfmult_of_boundary(size_type n) const
641  { return *(contact_boundaries[n].mflambda); }
642  const getfem::mesh_im &mim_of_boundary(size_type n) const
643  { return *(contact_boundaries[n].mim); }
644  size_type nb_variables(void) const { return Us.size(); }
645  size_type nb_multipliers(void) const { return lambdas.size(); }
646  const std::string &varname(size_type i) const { return Unames[i]; }
647  const std::string &multname(size_type i) const { return lambdanames[i]; }
648  const model_real_plain_vector &disp_of_boundary(size_type n) const
649  { return ext_Us[contact_boundaries[n].ind_U]; }
650  const model_real_plain_vector &disp0_of_boundary(size_type n) const
651  { return ext_Ws[contact_boundaries[n].ind_U]; }
652  const model_real_plain_vector &mult_of_boundary(size_type n) const
653  { return ext_lambdas[contact_boundaries[n].ind_lambda]; }
654  size_type region_of_boundary(size_type n) const
655  { return contact_boundaries[n].region; }
656  const std::string &varname_of_boundary(size_type n) const
657  { return Unames[contact_boundaries[n].ind_U]; }
658  size_type ind_varname_of_boundary(size_type n) const
659  { return contact_boundaries[n].ind_U; }
660  const std::string &multname_of_boundary(size_type n) const {
661  static const std::string vname;
662  size_type ind = contact_boundaries[n].ind_lambda;
663  return (ind == size_type(-1)) ? vname : lambdanames[ind];
664  }
665  size_type ind_multname_of_boundary(size_type n) const
666  { return contact_boundaries[n].ind_lambda; }
667  size_type nb_boundaries(void) const { return contact_boundaries.size(); }
668  bool is_self_contact(void) const { return self_contact; }
669  bool is_slave_boundary(size_type n) const { return contact_boundaries[n].slave; }
670  void set_raytrace(bool b) { raytrace = b; }
671  void set_nodes_mode(int m) { nodes_mode = m; }
672  size_type nb_contact_pairs(void) const { return contact_pairs.size(); }
673  const contact_pair &get_contact_pair(size_type i)
674  { return contact_pairs[i]; }
675 
676  multi_contact_frame(size_type NN, scalar_type r_dist,
677  bool dela = true, bool selfc = true,
678  scalar_type cut_a = 0.3, bool rayt = false,
679  int fem_nodes = 0, bool refc = false);
680  multi_contact_frame(const model &md, size_type NN, scalar_type r_dist,
681  bool dela = true, bool selfc = true,
682  scalar_type cut_a = 0.3, bool rayt = false,
683  int fem_nodes = 0, bool refc = false);
684 
685  size_type add_obstacle(const std::string &obs);
686 
687  size_type add_slave_boundary(const getfem::mesh_im &mim,
688  const getfem::mesh_fem *mfu,
689  const model_real_plain_vector *U,
690  size_type reg,
691  const getfem::mesh_fem *mflambda = 0,
692  const model_real_plain_vector *lambda = 0,
693  const model_real_plain_vector *w = 0,
694  const std::string &varname = std::string(),
695  const std::string &multname = std::string(),
696  const std::string &wname = std::string());
697 
698  size_type add_slave_boundary(const getfem::mesh_im &mim, size_type reg,
699  const std::string &varname,
700  const std::string &multname = std::string(),
701  const std::string &wname = std::string());
702 
703 
704  size_type add_master_boundary(const getfem::mesh_im &mim,
705  const getfem::mesh_fem *mfu,
706  const model_real_plain_vector *U,
707  size_type reg,
708  const getfem::mesh_fem *mflambda = 0,
709  const model_real_plain_vector *lambda = 0,
710  const model_real_plain_vector *w = 0,
711  const std::string &varname = std::string(),
712  const std::string &multname = std::string(),
713  const std::string &wname = std::string());
714 
715  size_type add_master_boundary(const getfem::mesh_im &mim, size_type reg,
716  const std::string &varname,
717  const std::string &multname = std::string(),
718  const std::string &wname = std::string());
719 
720 
721 
722  // The whole process of the computation of contact pairs
723  // Contact pairs are seached for a certain boundary (master or slave,
724  // depending on the contact algorithm) on the master ones. If contact pairs
725  // are searched for a master boundary, self-contact is taken into account
726  // if the flag 'self_contact' is set to 'true'. Self-contact is never taken
727  // into account for a slave boundary.
728  void compute_contact_pairs(void);
729 
730  };
731 
732 
733 
734  //=========================================================================
735  //
736  // Raytracing Interpolate transformation
737  //
738  //=========================================================================
739 
740  /** Add a raytracing interpolate transformation called 'transname' to a model
741  to be used by the generic assembly bricks.
742  */
743  void add_raytracing_transformation
744  (model &md, const std::string &transname, scalar_type release_distance);
745 
746  /** Add a raytracing interpolate transformation called 'transname' to a
747  workspace to be used by the generic assembly bricks.
748  */
749  void add_raytracing_transformation
750  (ga_workspace &workspace, const std::string &transname,
751  scalar_type release_distance);
752 
753  /** Add a master boundary with corresponding displacement variable
754  'dispname' on a specific boundary 'region' to an existing raytracing
755  interpolate transformation called 'transname'.
756  */
757  void add_master_contact_boundary_to_raytracing_transformation
758  (model &md, const std::string &transname,const mesh &m,
759  const std::string &dispname, size_type region);
760 
761  /** Add a slave boundary with corresponding displacement variable
762  'dispname' on a specific boundary 'region' to an existing raytracing
763  interpolate transformation called 'transname'.
764  */
765  void add_slave_contact_boundary_to_raytracing_transformation
766  (model &md, const std::string &transname, const mesh &m,
767  const std::string &dispname, size_type region);
768 
769  /** Add a master boundary with corresponding displacement variable
770  'dispname' on a specific boundary 'region' to an existing raytracing
771  interpolate transformation called 'transname'.
772  */
773  void add_master_contact_boundary_to_raytracing_transformation
774  (ga_workspace &workspace, const std::string &transname, const mesh &m,
775  const std::string &dispname, size_type region);
776 
777  /** Add a slave boundary with corresponding displacement variable
778  'dispname' on a specific boundary 'region' to an existing raytracing
779  interpolate transformation called 'transname'.
780  */
781  void add_slave_contact_boundary_to_raytracing_transformation
782  (ga_workspace &workspace, const std::string &transname, const mesh &m,
783  const std::string &dispname, size_type region);
784 
785  /** Add a rigid obstacle whose geometry corresponds to the zero level-set
786  of the high-level generic assembly expression `expr`
787  to an existing raytracing interpolate transformation called 'transname'.
788  In `expr`, the current position is denoted `X` with components
789  `X(1), X(2), ...`. It is also allowed to use `x` instead of `X(1)`,
790  `y` instead of `X(2)`, `z` instead of `X(3)` and `w` instead of `X(4)`.
791  */
792  void add_rigid_obstacle_to_raytracing_transformation
793  (model &md, const std::string &transname,
794  const std::string &expr, size_type N);
795 
796  /** Add a rigid obstacle whose geometry corresponds to the zero level-set
797  of the high-level generic assembly expression 'expr'
798  to an existing raytracing interpolate transformation called 'transname'.
799  */
800  void add_rigid_obstacle_to_raytracing_transformation
801  (ga_workspace &workspace, const std::string &transname,
802  const std::string &expr, size_type N);
803 
804  //=========================================================================
805  //
806  // Projection Interpolate transformation
807  //
808  //=========================================================================
809 
810  /** Add a projection interpolate transformation called 'transname' to a model
811  to be used by the generic assembly bricks.
812  */
813  void add_projection_transformation
814  (model &md, const std::string &transname, scalar_type release_distance);
815 
816  /** Add a projection interpolate transformation called 'transname' to a
817  workspace to be used by the generic assembly bricks.
818  */
819  void add_projection_transformation
820  (ga_workspace &workspace, const std::string &transname,
821  scalar_type release_distance);
822 
823  /** Add a master boundary with corresponding displacement variable
824  'dispname' on a specific boundary 'region' to an existing projection
825  interpolate transformation called 'transname'.
826  */
827  void add_master_contact_boundary_to_projection_transformation
828  (model &md, const std::string &transname,const mesh &m,
829  const std::string &dispname, size_type region);
830  /** Add a master boundary with corresponding displacement variable
831  'dispname' on a specific boundary 'region' to an existing projection
832  interpolate transformation called 'transname'.
833  */
834  void add_master_contact_boundary_to_projection_transformation
835  (ga_workspace &workspace, const std::string &transname, const mesh &m,
836  const std::string &dispname, size_type region);
837 
838  /** Add a slave boundary with corresponding displacement variable
839  'dispname' on a specific boundary 'region' to an existing projection
840  interpolate transformation called 'transname'.
841  */
842  void add_slave_contact_boundary_to_projection_transformation
843  (model &md, const std::string &transname, const mesh &m,
844  const std::string &dispname, size_type region);
845 
846  /** Add a slave boundary with corresponding displacement variable
847  'dispname' on a specific boundary 'region' to an existing projection
848  interpolate transformation called 'transname'.
849  */
850  void add_slave_contact_boundary_to_projection_transformation
851  (ga_workspace &workspace, const std::string &transname, const mesh &m,
852  const std::string &dispname, size_type region);
853 
854  /** Add a rigid obstacle whose geometry corresponds to the zero level-set
855  of the high-level generic assembly expression `expr`
856  to an existing projection interpolate transformation called 'transname'.
857  In `expr`, the current position is denoted `X` with components
858  `X(1), X(2), ...`. It is also allowed to use `x` instead of `X(1)`,
859  `y` instead of `X(2)`, `z` instead of `X(3)` and `w` instead of `X(4)`.
860  */
861  void add_rigid_obstacle_to_projection_transformation
862  (model &md, const std::string &transname,
863  const std::string &expr, size_type N);
864 
865  /** Add a rigid obstacle whose geometry corresponds to the zero level-set
866  of the high-level generic assembly expression 'expr'
867  to an existing projection interpolate transformation called 'transname'.
868  */
869  void add_rigid_obstacle_to_projection_transformation
870  (ga_workspace &workspace, const std::string &transname,
871  const std::string &expr, size_type N);
872 } /* end of namespace getfem. */
873 
874 
875 #endif /* GETFEM_CONTACT_AND_FRICTION_COMMON_H__ */
bgeot_rtree.h
region-tree for window/point search on a set of rectangles.
bgeot::rtree
Balanced tree of n-dimensional rectangles.
Definition: bgeot_rtree.h:98
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:152
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem_assembling_tensors.h
Generic assembly implementation.
getfem_mesher.h
An experimental mesher.
getfem_models.h
Model representation in Getfem.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:978
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:558
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:210
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:264
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1277
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:73
bgeot::compute_normal
base_small_vector compute_normal(const geotrans_interpolation_context &c, size_type face)
norm of returned vector is the ratio between the face surface on the real element and the face surfac...
Definition: bgeot_geometric_trans.cc:1082
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:47
getfem::add_rigid_obstacle_to_raytracing_transformation
void add_rigid_obstacle_to_raytracing_transformation(model &md, const std::string &transname, const std::string &expr, size_type N)
Add a rigid obstacle whose geometry corresponds to the zero level-set of the high-level generic assem...
Definition: getfem_contact_and_friction_common.cc:2482
getfem::add_master_contact_boundary_to_projection_transformation
void add_master_contact_boundary_to_projection_transformation(model &md, const std::string &transname, const mesh &m, const std::string &dispname, size_type region)
Add a master boundary with corresponding displacement variable 'dispname' on a specific boundary 'reg...
Definition: getfem_contact_and_friction_common.cc:2516
getfem::add_raytracing_transformation
void add_raytracing_transformation(model &md, const std::string &transname, scalar_type release_distance)
Add a raytracing interpolate transformation called 'transname' to a model to be used by the generic a...
Definition: getfem_contact_and_friction_common.cc:2428
getfem::add_rigid_obstacle_to_projection_transformation
void add_rigid_obstacle_to_projection_transformation(model &md, const std::string &transname, const std::string &expr, size_type N)
Add a rigid obstacle whose geometry corresponds to the zero level-set of the high-level generic assem...
Definition: getfem_contact_and_friction_common.cc:2556
getfem::add_master_contact_boundary_to_raytracing_transformation
void add_master_contact_boundary_to_raytracing_transformation(model &md, const std::string &transname, const mesh &m, const std::string &dispname, size_type region)
Add a master boundary with corresponding displacement variable 'dispname' on a specific boundary 'reg...
Definition: getfem_contact_and_friction_common.cc:2442
getfem::add_projection_transformation
void add_projection_transformation(model &md, const std::string &transname, scalar_type release_distance)
Add a projection interpolate transformation called 'transname' to a model to be used by the generic a...
Definition: getfem_contact_and_friction_common.cc:2502
getfem::add_slave_contact_boundary_to_projection_transformation
void add_slave_contact_boundary_to_projection_transformation(model &md, const std::string &transname, const mesh &m, const std::string &dispname, size_type region)
Add a slave boundary with corresponding displacement variable 'dispname' on a specific boundary 'regi...
Definition: getfem_contact_and_friction_common.cc:2526
getfem::add_slave_contact_boundary_to_raytracing_transformation
void add_slave_contact_boundary_to_raytracing_transformation(model &md, const std::string &transname, const mesh &m, const std::string &dispname, size_type region)
Add a slave boundary with corresponding displacement variable 'dispname' on a specific boundary 'regi...
Definition: getfem_contact_and_friction_common.cc:2452