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