GetFEM  5.5
getfem_contact_and_friction_integral.cc
1 /*===========================================================================
2 
3  Copyright (C) 2011-2026 Yves Renard, Konstantinos Poulios.
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, see https://www.gnu.org/licenses/.
18 
19 ===========================================================================*/
20 
21 #include "getfem/bgeot_rtree.h"
22 #include "getfem/getfem_contact_and_friction_integral.h"
23 #include "getfem/getfem_contact_and_friction_common.h"
24 #include "getfem/getfem_projected_fem.h"
25 #include "gmm/gmm_condition_number.h"
26 
27 namespace getfem {
28 
29 
30 
31  template <typename T> inline static T Heav(T a)
32  { return (a < T(0)) ? T(0) : T(1); }
33 
34 
35  //=========================================================================
36  //
37  // Basic non linear term used for contact bricks.
38  //
39  //=========================================================================
40 
41  void contact_nonlinear_term::adjust_tensor_size(void) {
42  sizes_.resize(1); sizes_[0] = 1;
43  switch (option) {
44  // one-dimensional tensors [N]
45  case RHS_U_V1: case RHS_U_V2: case RHS_U_V4:
46  case RHS_U_V5: case RHS_U_FRICT_V6: case RHS_U_FRICT_V7:
47  case RHS_U_FRICT_V8: case RHS_U_FRICT_V1:
48  case RHS_U_FRICT_V4: case RHS_U_FRICT_V5:
49  case RHS_L_FRICT_V1: case RHS_L_FRICT_V2: case RHS_L_FRICT_V4:
50  case K_UL_V1: case K_UL_V2: case K_UL_V3:
51  case UZAWA_PROJ_FRICT: case UZAWA_PROJ_FRICT_SAXCE:
52  sizes_[0] = N; break;
53  // two-dimensional tensors [N x N]
54  case K_UU_V1: case K_UU_V2:
55  case K_UL_FRICT_V1: case K_UL_FRICT_V2: case K_UL_FRICT_V3:
56  case K_UL_FRICT_V4: case K_UL_FRICT_V5:
57  case K_UL_FRICT_V7: case K_UL_FRICT_V8:
58  case K_LL_FRICT_V1: case K_LL_FRICT_V2: case K_LL_FRICT_V4:
59  case K_UU_FRICT_V1: case K_UU_FRICT_V2:
60  case K_UU_FRICT_V3: case K_UU_FRICT_V4: case K_UU_FRICT_V5:
61  sizes_.resize(2); sizes_[0] = sizes_[1] = N; break;
62  }
63 
64  // adjust temporary variables sizes
65  lnt.resize(N); lt.resize(N); zt.resize(N); no.resize(N);
66  aux1.resize(1); auxN.resize(N); V.resize(N);
67  gmm::resize(GP, N, N);
68  }
69 
70  void contact_nonlinear_term::friction_law
71  (scalar_type p, scalar_type &tau) {
72  tau = (p > scalar_type(0)) ? tau_adh + f_coeff * p : scalar_type(0);
73  if (tau > tresca_lim) tau = tresca_lim;
74  }
75 
76  void contact_nonlinear_term::friction_law
77  (scalar_type p, scalar_type &tau, scalar_type &tau_grad) {
78  if (p <= scalar_type(0)) {
79  tau = scalar_type(0);
80  tau_grad = scalar_type(0);
81  }
82  else {
83  tau = tau_adh + f_coeff * p;
84  if (tau > tresca_lim) {
85  tau = tresca_lim;
86  tau_grad = scalar_type(0);
87  }
88  else
89  tau_grad = f_coeff;
90  }
91  }
92 
93  void contact_nonlinear_term::compute
94  (fem_interpolation_context &/* ctx */, bgeot::base_tensor &t) {
95 
96  t.adjust_sizes(sizes_);
97  scalar_type e, augm_ln, rho, rho_grad;
98  dim_type i, j;
99  bool coulomb;
100 
101  switch (option) {
102 
103  // scalar tensors [1]
104 
105  case RHS_L_V1:
106  t[0] = (ln+gmm::neg(ln-r*(un - g)))/r; break;
107  case RHS_L_V2:
108  t[0] = (un-g) + gmm::pos(ln)/r; break;
109 
110  case K_LL_V1:
111  t[0] = (Heav(r*(un-g)-ln) - scalar_type(1))/r; break;
112  case K_LL_V2:
113  t[0] = -Heav(ln)/r; break;
114 
115  case UZAWA_PROJ:
116  t[0] = -gmm::neg(ln - r*(un - g)); break;
117 
118  case CONTACT_FLAG:
119  // here ln is expected to be a threshold value expressing a penetration
120  // (positive value) or separation (negative value) distance
121  t[0] = Heav(un-g - ln); break;
122  case CONTACT_PRESSURE:
123  t[0] = -ln; break;
124 
125  // one-dimensional tensors [N]
126 
127  case RHS_U_V1:
128  for (i=0; i<N; ++i) t[i] = ln * no[i];
129  break;
130  case RHS_U_V2:
131  e = -gmm::neg(ln-r*(un - g));
132  for (i=0; i<N; ++i) t[i] = e * no[i];
133  break;
134  case RHS_U_V4:
135  e = -gmm::neg(ln);
136  for (i=0; i<N; ++i) t[i] = e * no[i];
137  break;
138  case RHS_U_V5:
139  e = - gmm::pos(un-g) * r;
140  for (i=0; i<N; ++i) t[i] = e * no[i];
141  break;
142  case RHS_U_FRICT_V6:
143  e = gmm::neg(ln-r*(un - g));
144  friction_law(e, rho);
145  auxN = lt - zt; ball_projection(auxN, rho);
146  for (i=0; i<N; ++i) t[i] = auxN[i] - e*no[i];
147  break;
148  case RHS_U_FRICT_V7:
149  e = gmm::neg(-r*(un - g));
150  friction_law(e, rho);
151  auxN = - zt; ball_projection(auxN, rho);
152  for (i=0; i<N; ++i) t[i] = auxN[i] - e*no[i];
153  break;
154  case RHS_U_FRICT_V8: // ignores friction_law, assumes pure Coulomb friction
155  auxN = lnt - (r*(un-g) - f_coeff * gmm::vect_norm2(zt)) * no - zt;
156  De_Saxce_projection(auxN, no, f_coeff);
157  for (i=0; i<N; ++i) t[i] = auxN[i];
158  break;
159  case RHS_U_FRICT_V1:
160  for (i=0; i<N; ++i) t[i] = lnt[i];
161  break;
162  case RHS_U_FRICT_V4:
163  e = gmm::neg(ln);
164  // if (e > 0. && ctx.xreal()[1] > 1.)
165  // cout << "x = " << ctx.xreal() << " e = " << e << endl;
166  friction_law(e, rho);
167  auxN = lt; ball_projection(auxN, rho);
168  // if (gmm::vect_norm2(auxN) > 0. && ctx.xreal()[1] > 1.)
169  // cout << "x = " << ctx.xreal() << " auxN = " << auxN << endl;
170  for (i=0; i<N; ++i) t[i] = auxN[i] - e*no[i];
171  break;
172  case RHS_U_FRICT_V5: // ignores friction_law, assumes pure Coulomb friction
173  auxN = lnt; De_Saxce_projection(auxN, no, f_coeff);
174  for (i=0; i<N; ++i) t[i] = auxN[i];
175  break;
176  case RHS_L_FRICT_V1:
177  e = gmm::neg(ln-r*(un-g));
178  friction_law(e, rho);
179  auxN = zt - lt; ball_projection(auxN, rho); auxN += lt;
180  for (i=0; i<N; ++i) t[i] = ((e+ln)*no[i] + auxN[i])/ r;
181  break;
182  case RHS_L_FRICT_V2:
183  e = r*(un-g) + gmm::pos(ln);
184  friction_law(gmm::neg(ln), rho);
185  auxN = lt; ball_projection(auxN, rho);
186  for (i=0; i<N; ++i) t[i] = (no[i]*e + zt[i] + lt[i] - auxN[i])/r;
187  break;
188  case RHS_L_FRICT_V4: // ignores friction_law, assumes pure Coulomb friction
189  auxN = lnt;
190  De_Saxce_projection(auxN, no, f_coeff);
191  auxN -= lnt + (r*(un-g) - f_coeff * gmm::vect_norm2(zt)) * no + zt;
192  for (i=0; i<N; ++i) t[i] = -auxN[i]/r;
193  break;
194  case K_UL_V1:
195  for (i=0; i<N; ++i) t[i] = -no[i];
196  break;
197  case K_UL_V2 :
198  e = -Heav(-ln); //Heav(ln)-scalar_type(1);
199  for (i=0; i<N; ++i) t[i] = e*no[i];
200  break;
201  case K_UL_V3:
202  e = -Heav(r*(un-g)-ln);
203  for (i=0; i<N; ++i) t[i] = e*no[i];
204  break;
205  case UZAWA_PROJ_FRICT:
206  e = gmm::neg(ln - r*(un - g));
207  friction_law(e, rho);
208  auxN = lt - zt; ball_projection(auxN, rho);
209  for (i=0; i<N; ++i) t[i] = auxN[i] - e*no[i];
210  break;
211  case UZAWA_PROJ_FRICT_SAXCE: // ignores friction_law, assumes pure Coulomb friction
212  auxN = lnt - (r*(un-g) - f_coeff * gmm::vect_norm2(zt)) * no - zt;
213  De_Saxce_projection(auxN, no, f_coeff);
214  for (i=0; i<N; ++i) t[i] = auxN[i];
215  break;
216 
217  // two-dimensional tensors [N x N]
218 
219  case K_UU_V1:
220  e = r * Heav(un - g);
221  for (i=0; i<N; ++i) for (j=0; j<N; ++j) t[i*N+j] = e * no[i] * no[j];
222  break;
223  case K_UU_V2:
224  e = r * Heav(r*(un - g)-ln);
225  for (i=0; i<N; ++i) for (j=0; j<N; ++j) t[i*N+j] = e * no[i] * no[j];
226  break;
227 
228  case K_UL_FRICT_V1:
229  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
230  t[i*N+j] = ((i == j) ? -scalar_type(1) : scalar_type(0));
231  break;
232  case K_UL_FRICT_V2:
233  friction_law(gmm::neg(ln), rho, rho_grad);
234  ball_projection_grad(lt, rho, GP);
235  e = gmm::vect_sp(GP, no, no) - Heav(-ln);
236  coulomb = (rho_grad > 0) && bool(Heav(-ln));
237  if (coulomb) ball_projection_grad_r(lt, rho, V);
238  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
239  t[i*N+j] = no[i]*no[j]*e - GP(i,j) +
240  (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0));
241  break;
242  case K_UL_FRICT_V3:
243  augm_ln = ln - r*(un-g);
244  friction_law(gmm::neg(augm_ln), rho, rho_grad);
245  auxN = lt - zt;
246  ball_projection_grad(auxN, rho, GP);
247  e = gmm::vect_sp(GP, no, no) - Heav(-augm_ln);
248  coulomb = (rho_grad > 0) && bool(Heav(-augm_ln));
249  if (coulomb) ball_projection_grad_r(auxN, rho, V);
250  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
251  t[i*N+j] = no[i]*no[j]*e - GP(i,j) +
252  (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0));
253  break;
254  case K_UL_FRICT_V4:
255  augm_ln = ln - r*(un-g);
256  friction_law(gmm::neg(augm_ln), rho, rho_grad);
257  auxN = lt - zt;
258  ball_projection_grad(auxN, rho, GP); gmm::scale(GP, alpha);
259  e = gmm::vect_sp(GP, no, no) - Heav(-augm_ln);
260  coulomb = (rho_grad > 0) && bool(Heav(-augm_ln));
261  if (coulomb) ball_projection_grad_r(auxN, rho, V);
262  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
263  t[i*N+j] = no[i]*no[j]*e - GP(i,j) +
264  (coulomb ? rho_grad*V[i]*no[j] : scalar_type(0));
265  break;
266  case K_UL_FRICT_V5:
267  e = alpha - scalar_type(1);
268  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
269  t[i*N+j] = no[i]*no[j]*e - ((i == j) ? alpha : scalar_type(0));
270  break;
271  case K_UL_FRICT_V7: // ignores friction_law, assumes pure Coulomb friction
272  De_Saxce_projection_grad(lnt, no, f_coeff, GP);
273  for (i=0; i<N; ++i) for (j=0; j<N; ++j) t[i*N+j] = -GP(j,i);
274  break;
275  case K_UL_FRICT_V8: // ignores friction_law, assumes pure Coulomb friction
276  {
277  scalar_type nzt = gmm::vect_norm2(zt);
278  gmm::copy(gmm::identity_matrix(), GP); gmm::scale(GP, alpha);
279  gmm::rank_one_update(GP, gmm::scaled(no, scalar_type(1)-alpha), no);
280  if (nzt != scalar_type(0))
281  gmm::rank_one_update(GP, gmm::scaled(no, -f_coeff*alpha/nzt), zt);
282  for (i=0; i<N; ++i) for (j=0; j<N; ++j) t[i*N+j] = - GP(i,j);
283  }
284  break;
285  case K_LL_FRICT_V1:
286  augm_ln = ln - r*(un-g);
287  friction_law(gmm::neg(augm_ln), rho, rho_grad);
288  auxN = lt - zt;
289  ball_projection_grad(auxN, rho, GP);
290  e = Heav(-augm_ln) - gmm::vect_sp(GP, no, no);
291  coulomb = (rho_grad > 0) && bool(Heav(-augm_ln));
292  if (coulomb) ball_projection_grad_r(auxN, rho, V);
293  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
294  t[i*N+j] = (no[i]*no[j]*e + GP(i,j)
295  - ((i == j) ? scalar_type(1) : scalar_type(0))
296  - (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0))) / r;
297  break;
298  case K_LL_FRICT_V2:
299  friction_law(gmm::neg(ln), rho, rho_grad);
300  ball_projection_grad(lt, rho, GP);
301  e = Heav(-ln) - gmm::vect_sp(GP, no, no);
302  coulomb = (rho_grad > 0) && bool(Heav(-ln));
303  if (coulomb) ball_projection_grad_r(lt, rho, V);
304  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
305  t[i*N+j] = (no[i]*no[j]*e + GP(i,j)
306  - ((i == j) ? scalar_type(1) : scalar_type(0))
307  - (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0))) / r;
308  break;
309  case K_LL_FRICT_V4: // ignores friction_law, assumes pure Coulomb friction
310  De_Saxce_projection_grad(lnt, no, f_coeff, GP);
311  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
312  t[i*N+j] = (GP(i,j) - ((i == j) ? scalar_type(1) : scalar_type(0)))/r;
313  break;
314  case K_UU_FRICT_V1:
315  e = r * Heav(r*(un-g)-ln);
316  for (i=0; i<N; ++i) for (j=0; j<N; ++j) t[i*N+j] = no[i]*no[j]*e;
317  break;
318  case K_UU_FRICT_V2:
319  friction_law(-ln, rho, rho_grad);
320  auxN = lt - zt;
321  ball_projection_grad(auxN, rho, GP); gmm::scale(GP, alpha);
322  e = Heav(r*(un-g)-ln) - gmm::vect_sp(GP, no, no);
323  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
324  t[i*N+j] = r*(no[i]*no[j]*e + GP(i,j));
325  break;
326  case K_UU_FRICT_V3:
327  case K_UU_FRICT_V4:
328  augm_ln = (option == K_UU_FRICT_V3) ? ln - r*(un-g) : - r*(un-g);
329  auxN = (option == K_UU_FRICT_V3) ? lt - zt : -zt;
330  friction_law(gmm::neg(augm_ln), rho, rho_grad);
331  ball_projection_grad(auxN, rho, GP); gmm::scale(GP, alpha);
332  e = Heav(-augm_ln) - gmm::vect_sp(GP, no, no);
333  coulomb = (rho_grad > 0) && bool(Heav(-augm_ln));
334  if (coulomb) ball_projection_grad_r(auxN, rho, V);
335  for (i=0; i<N; ++i) for (j=0; j<N; ++j)
336  t[i*N+j] = r*(no[i]*no[j]*e + GP(i,j)
337  - (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0)));
338  break;
339  case K_UU_FRICT_V5: // ignores friction_law, assumes pure Coulomb friction
340  {
341  scalar_type nzt = gmm::vect_norm2(zt);
342  auxN = lnt - (r*(un-g) - f_coeff * nzt) * no - zt;
343  base_matrix A(N, N), B(N, N);
344  De_Saxce_projection_grad(auxN, no, f_coeff, A);
345  gmm::copy(gmm::identity_matrix(), B); gmm::scale(B, alpha);
346  gmm::rank_one_update(B, gmm::scaled(no, scalar_type(1)-alpha), no);
347  if (nzt != scalar_type(0))
348  gmm::rank_one_update(B, gmm::scaled(no, -f_coeff*alpha/nzt), zt);
349  gmm::mult(A, B, GP);
350  for (i=0; i<N; ++i) for (j=0; j<N; ++j) t[i*N+j] = r*GP(j,i);
351  }
352  break;
353  default : GMM_ASSERT1(false, "Invalid option");
354  }
355  }
356 
357 
358  //=========================================================================
359  //
360  // Non linear term used for contact with rigid obstacle bricks.
361  //
362  //=========================================================================
363 
364  void contact_rigid_obstacle_nonlinear_term::prepare
365  (fem_interpolation_context& ctx, size_type nb) {
366  size_type cv = ctx.convex_num();
367 
368  switch (nb) { // last is computed first
369  case 1 : // calculate [un] and [zt] interpolating [U],[WT],[VT] on [mf_u]
370  slice_vector_on_basic_dof_of_element(mf_u, U, cv, coeff);
371  ctx.pf()->interpolation(ctx, coeff, V, N);
372  un = gmm::vect_sp(V, no);
373  if (!contact_only) {
374  if (gmm::vect_size(WT) == gmm::vect_size(U)) {
375  slice_vector_on_basic_dof_of_element(mf_u, WT, cv, coeff);
376  ctx.pf()->interpolation(ctx, coeff, auxN, N);
377  auxN -= gmm::vect_sp(auxN, no) * no;
378  if (gmm::vect_size(VT) == gmm::vect_size(U)) {
379  slice_vector_on_basic_dof_of_element(mf_u, VT, cv, coeff);
380  ctx.pf()->interpolation(ctx, coeff, vt, N);
381  vt -= gmm::vect_sp(vt, no) * no;
382  // zt = r*(alpha*(u_T-w_T) + (1-gamma)*v_T)
383  zt = (((V - un * no) - auxN) * alpha + vt * (1 - gamma)) * r;
384  } else {
385  // zt = r*alpha*(u_T-w_T)
386  zt = ((V - un * no) - auxN) * (r * alpha);
387  }
388  } else {
389  // zt = r*alpha*u_T
390  zt = (V - un * no) * (r * alpha);
391  }
392  }
393  break;
394 
395  case 2 : // calculate [g] and [no] interpolating [obs] on [mf_obs]
396  // calculate [ln] and [lt] from [lnt] and [no]
397  slice_vector_on_basic_dof_of_element(mf_obs, obs, cv, coeff);
398  ctx.pf()->interpolation_grad(ctx, coeff, grad, 1);
399  gmm::copy(gmm::mat_row(grad, 0), no);
400  no /= -gmm::vect_norm2(no);
401  ctx.pf()->interpolation(ctx, coeff, aux1, 1);
402  g = aux1[0];
403 
404  if (!contact_only && pmf_lambda) {
405  ln = gmm::vect_sp(lnt, no);
406  lt = lnt - ln * no;
407  }
408 
409  break;
410 
411  case 3 : // calculate [ln] or [lnt] interpolating [lambda] on [mf_lambda]
412  if (pmf_lambda) {
413  slice_vector_on_basic_dof_of_element(*pmf_lambda, lambda, cv, coeff);
414  if (contact_only) {
415  ctx.pf()->interpolation(ctx, coeff, aux1, 1);
416  ln = aux1[0];
417  } else {
418  ctx.pf()->interpolation(ctx, coeff, lnt, N);
419  }
420  }
421  break;
422 
423  case 4 : // calculate [f_coeff] interpolating [friction_coeff] on [mf_coeff]
424  // calculate [tau_adh] interpolating [tau_adhesion] on [mf_coeff]
425  // calculate [tresca_lim] interpolating [tresca_limit] on [mf_coeff]
426  GMM_ASSERT1(!contact_only, "Invalid friction option");
427  if (pmf_coeff) {
428  slice_vector_on_basic_dof_of_element(*pmf_coeff, friction_coeff, cv, coeff);
429  ctx.pf()->interpolation(ctx, coeff, aux1, 1);
430  f_coeff = aux1[0];
431  if (gmm::vect_size(tau_adhesion)) {
432  slice_vector_on_basic_dof_of_element(*pmf_coeff, tau_adhesion, cv, coeff);
433  ctx.pf()->interpolation(ctx, coeff, aux1, 1);
434  tau_adh = aux1[0];
435  if (gmm::vect_size(tresca_limit)) {
436  slice_vector_on_basic_dof_of_element(*pmf_coeff, tresca_limit, cv, coeff);
437  ctx.pf()->interpolation(ctx, coeff, aux1, 1);
438  tresca_lim = aux1[0];
439  }
440  }
441  }
442  break;
443 
444  default : GMM_ASSERT1(false, "Invalid option");
445  }
446 
447  }
448 
449 
450  //=========================================================================
451  //
452  // Non linear term used for contact between non-matching meshes bricks.
453  //
454  //=========================================================================
455 
456  void contact_nonmatching_meshes_nonlinear_term::prepare
457  (fem_interpolation_context& ctx, size_type nb) {
458 
459  size_type cv = ctx.convex_num();
460 
461  // - this method is called for nb=4,3,2,1 corresponding to
462  // NonLin(#i0,#i1,#i2,#i3,#i4) before the compute() method is called
463  // for i0
464  // - in each call ctx.pf() corresponds to the fem of the cv convex
465  // on the i4,i3,i2 and i1 mesh_fem respectively
466  // - this method expects that i1,i2,i3 and i4 mesh_fems will correspond
467  // to mf_u1, mf_u2, mf_lambda and mf_coeff
468 
469  switch (nb) { // last is computed first
470  case 1 : // calculate [un] and [zt] interpolating [U1],[WT1] on [mf_u1]
471  // and subtracting [un] and [zt] calculated on [mf_u2]
472  slice_vector_on_basic_dof_of_element(mf_u1, U1, cv, coeff);
473  ctx.pf()->interpolation(ctx, coeff, V, N);
474  {
475  scalar_type un1 = gmm::vect_sp(V, no);
476  if (!contact_only) {
477  if (gmm::vect_size(WT1) == gmm::vect_size(U1)) {
478  slice_vector_on_basic_dof_of_element(mf_u1, WT1, cv, coeff);
479  ctx.pf()->interpolation(ctx, coeff, auxN, N);
480  auxN -= gmm::vect_sp(auxN, no) * no;
481  zt = ((V - un1 * no) - auxN) * (r * alpha) - zt; // = zt1 - zt2 , with zt = r*alpha*(u_T-w_T)
482  } else {
483  zt = (V - un1 * no) * (r * alpha) - zt; // = zt1 - zt2 , with zt = r*alpha*u_T
484  }
485  }
486  un = un1 - un; // = un1 - un2
487  }
488  break;
489 
490  case 2 : // calculate [g] and [no]
491  // calculate [ln] and [lt] from [lnt] and [no]
492  // calculate [un] and [zt] interpolating [U2],[WT2] on [mf_u2]
493  {
494  const projected_fem &pfe = dynamic_cast<const projected_fem&>(*ctx.pf());
495  pfe.projection_data(ctx, no, g);
496  gmm::scale(no, scalar_type(-1)); // pointing outwards from mf_u1
497  }
498 
499  if (!contact_only && pmf_lambda) {
500  ln = gmm::vect_sp(lnt, no);
501  lt = lnt - ln * no;
502  }
503 
504  slice_vector_on_basic_dof_of_element(mf_u2, U2, cv, coeff);
505  ctx.pf()->interpolation(ctx, coeff, V, N);
506  un = gmm::vect_sp(V, no);
507  if (!contact_only) {
508  if (gmm::vect_size(WT2) == gmm::vect_size(U2)) {
509  slice_vector_on_basic_dof_of_element(mf_u2, WT2, cv, coeff);
510  ctx.pf()->interpolation(ctx, coeff, auxN, N);
511  auxN -= gmm::vect_sp(auxN, no) * no;
512  zt = ((V - un * no) - auxN) * (r * alpha); // zt = r*alpha*(u_T-w_T)
513  } else {
514  zt = (V - un * no) * (r * alpha); // zt = r*alpha*u_T
515  }
516  }
517  break;
518 
519  case 3 : // calculate [ln] or [lnt] interpolating [lambda] on [mf_lambda]
520  if (pmf_lambda) {
521  slice_vector_on_basic_dof_of_element(*pmf_lambda, lambda, cv, coeff);
522  if (contact_only) {
523  ctx.pf()->interpolation(ctx, coeff, aux1, 1);
524  ln = aux1[0];
525  } else {
526  ctx.pf()->interpolation(ctx, coeff, lnt, N);
527  }
528  }
529  break;
530 
531  case 4 : // calculate [f_coeff] interpolating [friction_coeff] on [mf_coeff]
532  // calculate [tau_adh] interpolating [tau_adhesion] on [mf_coeff]
533  // calculate [tresca_lim] interpolating [tresca_limit] on [mf_coeff]
534  GMM_ASSERT1(!contact_only, "Invalid friction option");
535  if (pmf_coeff) {
536  slice_vector_on_basic_dof_of_element(*pmf_coeff, friction_coeff, cv, coeff);
537  ctx.pf()->interpolation(ctx, coeff, aux1, 1);
538  f_coeff = aux1[0];
539  if (gmm::vect_size(tau_adhesion)) {
540  slice_vector_on_basic_dof_of_element(*pmf_coeff, tau_adhesion, cv, coeff);
541  ctx.pf()->interpolation(ctx, coeff, aux1, 1);
542  tau_adh = aux1[0];
543  if (gmm::vect_size(tresca_limit)) {
544  slice_vector_on_basic_dof_of_element(*pmf_coeff, tresca_limit, cv, coeff);
545  ctx.pf()->interpolation(ctx, coeff, aux1, 1);
546  tresca_lim = aux1[0];
547  }
548  }
549  }
550  break;
551 
552  default : GMM_ASSERT1(false, "Invalid option");
553  }
554 
555  }
556 
557 
558  //=========================================================================
559  //
560  // Integral augmented Lagrangian brick (given obstacle, u, lambda).
561  //
562  //=========================================================================
563 
564  template<typename MAT, typename VECT1>
565  void asm_Alart_Curnier_contact_rigid_obstacle_tangent_matrix // frictionless
566  (MAT &Kul, MAT &Klu, MAT &Kll, MAT &Kuu,
567  const mesh_im &mim,
568  const getfem::mesh_fem &mf_u, const VECT1 &U,
569  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
570  const getfem::mesh_fem &mf_lambda, const VECT1 &lambda,
571  scalar_type r, const mesh_region &rg, int option = 1) {
572 
573  size_type subterm1 = (option == 3) ? K_UL_V2 : K_UL_V1;
574  size_type subterm2 = (option == 3) ? K_UL_V1 : K_UL_V3;
575  size_type subterm3 = (option == 3) ? K_LL_V2 : K_LL_V1;
576  size_type subterm4 = (option == 2) ? K_UU_V2 : K_UU_V1;
577 
578  contact_rigid_obstacle_nonlinear_term
579  nterm1(subterm1, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda),
580  nterm2(subterm2, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda),
581  nterm3(subterm3, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda),
582  nterm4(subterm4, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda);
583 
584  getfem::generic_assembly assem;
585  switch (option) {
586  case 1: case 3:
587  assem.set
588  ("M$1(#1,#3)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#1).Base(#3))(i,:,i,:); " // UL
589  "M$2(#3,#1)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3).vBase(#1))(i,:,:,i); " // LU
590  "M$3(#3,#3)+=comp(NonLin$3(#1,#1,#2,#3).Base(#3).Base(#3))(i,:,:)"); // LL
591  break;
592  case 2:
593  assem.set
594  ("M$1(#1,#3)+=comp(NonLin$2(#1,#1,#2,#3).vBase(#1).Base(#3))(i,:,i,:); " // UL
595  "M$3(#3,#3)+=comp(NonLin$3(#1,#1,#2,#3).Base(#3).Base(#3))(i,:,:);" // LL
596  "M$4(#1,#1)+=comp(NonLin$4(#1,#1,#2,#3).vBase(#1).vBase(#1))(i,j,:,i,:,j)"); // UU
597  break;
598  }
599  assem.push_mi(mim);
600  assem.push_mf(mf_u);
601  assem.push_mf(mf_obs);
602  assem.push_mf(mf_lambda);
603  assem.push_nonlinear_term(&nterm1);
604  assem.push_nonlinear_term(&nterm2);
605  assem.push_nonlinear_term(&nterm3);
606  assem.push_nonlinear_term(&nterm4);
607  assem.push_mat(Kul);
608  assem.push_mat(Klu);
609  assem.push_mat(Kll);
610  assem.push_mat(Kuu);
611  assem.assembly(rg);
612  }
613 
614  template<typename MAT, typename VECT1>
615  void asm_Alart_Curnier_contact_rigid_obstacle_tangent_matrix // with friction
616  (MAT &Kul, MAT &Klu, MAT &Kll, MAT &Kuu,
617  const mesh_im &mim,
618  const getfem::mesh_fem &mf_u, const VECT1 &U,
619  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
620  const getfem::mesh_fem &mf_lambda, const VECT1 &lambda,
621  const getfem::mesh_fem *pmf_coeff, const VECT1 *f_coeffs, scalar_type r,
622  scalar_type alpha, const VECT1 *WT,
623  scalar_type gamma, const VECT1 *VT,
624  const mesh_region &rg, int option = 1) {
625 
626  size_type subterm1, subterm2, subterm3;
627  switch (option) {
628  case 1 : subterm1 = K_UL_FRICT_V1; subterm2 = K_UL_FRICT_V4;
629  subterm3 = K_LL_FRICT_V1; break;
630  case 2 : subterm1 = K_UL_FRICT_V3; subterm2 = K_UL_FRICT_V4;
631  subterm3 = K_LL_FRICT_V1; break;
632  case 3 : subterm1 = K_UL_FRICT_V2; subterm2 = K_UL_FRICT_V5;
633  subterm3 = K_LL_FRICT_V2; break;
634  case 4 : subterm1 = K_UL_FRICT_V7; subterm2 = K_UL_FRICT_V8;
635  subterm3 = K_LL_FRICT_V4; break;
636  default : GMM_ASSERT1(false, "Incorrect option");
637  }
638 
639  size_type subterm4 = K_UU_FRICT_V3;
640 
641  contact_rigid_obstacle_nonlinear_term
642  nterm1(subterm1, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
643  pmf_coeff, f_coeffs, alpha, WT, gamma, VT),
644  nterm2(subterm2, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
645  pmf_coeff, f_coeffs, alpha, WT, gamma, VT),
646  nterm3(subterm3, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
647  pmf_coeff, f_coeffs, alpha, WT, gamma, VT),
648  nterm4(subterm4, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
649  pmf_coeff, f_coeffs, alpha, WT, gamma, VT);
650 
651  const std::string aux_fems = pmf_coeff ? "#1,#2,#3,#4" : "#1,#2,#3";
652 
653  getfem::generic_assembly assem;
654  switch (option) {
655  case 1: case 3: case 4:
656  assem.set
657  ("M$1(#1,#3)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1).vBase(#3))(i,j,:,i,:,j); " // UL
658  "M$2(#3,#1)+=comp(NonLin$2(#1," + aux_fems + ").vBase(#3).vBase(#1))(i,j,:,j,:,i); " // LU
659  "M$3(#3,#3)+=comp(NonLin$3(#1," + aux_fems + ").vBase(#3).vBase(#3))(i,j,:,i,:,j)"); // LL
660  break;
661  case 2:
662  assem.set
663  ("M$1(#1,#3)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1).vBase(#3))(i,j,:,i,:,j); " // UL
664  "M$2(#3,#1)+=comp(NonLin$2(#1," + aux_fems + ").vBase(#3).vBase(#1))(i,j,:,j,:,i); " // LU
665  "M$3(#3,#3)+=comp(NonLin$3(#1," + aux_fems + ").vBase(#3).vBase(#3))(i,j,:,i,:,j);" // LL
666  "M$4(#1,#1)+=comp(NonLin$4(#1," + aux_fems + ").vBase(#1).vBase(#1))(i,j,:,i,:,j)"); // UU
667  break;
668  }
669  assem.push_mi(mim);
670  assem.push_mf(mf_u);
671  assem.push_mf(mf_obs);
672  assem.push_mf(mf_lambda);
673  if (pmf_coeff)
674  assem.push_mf(*pmf_coeff);
675  assem.push_nonlinear_term(&nterm1);
676  assem.push_nonlinear_term(&nterm2);
677  assem.push_nonlinear_term(&nterm3);
678  assem.push_nonlinear_term(&nterm4);
679  assem.push_mat(Kul);
680  assem.push_mat(Klu);
681  assem.push_mat(Kll);
682  assem.push_mat(Kuu);
683  assem.assembly(rg);
684  }
685 
686  template<typename VECT1>
687  void asm_Alart_Curnier_contact_rigid_obstacle_rhs // frictionless
688  (VECT1 &Ru, VECT1 &Rl,
689  const mesh_im &mim,
690  const getfem::mesh_fem &mf_u, const VECT1 &U,
691  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
692  const getfem::mesh_fem &mf_lambda, const VECT1 &lambda,
693  scalar_type r, const mesh_region &rg, int option = 1) {
694 
695  size_type subterm1;
696  switch (option) {
697  case 1 : subterm1 = RHS_U_V1; break;
698  case 2 : subterm1 = RHS_U_V2; break;
699  case 3 : subterm1 = RHS_U_V4; break;
700  default : GMM_ASSERT1(false, "Incorrect option");
701  }
702  size_type subterm2 = (option == 3) ? RHS_L_V2 : RHS_L_V1;
703 
704  contact_rigid_obstacle_nonlinear_term
705  nterm1(subterm1, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda),
706  nterm2(subterm2, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda);
707 
708  getfem::generic_assembly assem;
709  assem.set("V$1(#1)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#1))(i,:,i); "
710  "V$2(#3)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3))(i,:)");
711  assem.push_mi(mim);
712  assem.push_mf(mf_u);
713  assem.push_mf(mf_obs);
714  assem.push_mf(mf_lambda);
715  assem.push_nonlinear_term(&nterm1);
716  assem.push_nonlinear_term(&nterm2);
717  assem.push_vec(Ru);
718  assem.push_vec(Rl);
719  assem.assembly(rg);
720 
721  }
722 
723  template<typename VECT1>
724  void asm_Alart_Curnier_contact_rigid_obstacle_rhs // with friction
725  (VECT1 &Ru, VECT1 &Rl,
726  const mesh_im &mim,
727  const getfem::mesh_fem &mf_u, const VECT1 &U,
728  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
729  const getfem::mesh_fem &mf_lambda, const VECT1 &lambda,
730  const getfem::mesh_fem *pmf_coeff, const VECT1 *f_coeffs, scalar_type r,
731  scalar_type alpha, const VECT1 *WT,
732  scalar_type gamma, const VECT1 *VT,
733  const mesh_region &rg, int option = 1) {
734 
735  size_type subterm1, subterm2;
736  switch (option) {
737  case 1 : subterm1 = RHS_U_FRICT_V1; subterm2 = RHS_L_FRICT_V1; break;
738  case 2 : subterm1 = RHS_U_FRICT_V6; subterm2 = RHS_L_FRICT_V1; break;
739  case 3 : subterm1 = RHS_U_FRICT_V4; subterm2 = RHS_L_FRICT_V2; break;
740  case 4 : subterm1 = RHS_U_FRICT_V5; subterm2 = RHS_L_FRICT_V4; break;
741  default : GMM_ASSERT1(false, "Incorrect option");
742  }
743 
744  contact_rigid_obstacle_nonlinear_term
745  nterm1(subterm1, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
746  pmf_coeff, f_coeffs, alpha, WT, gamma, VT),
747  nterm2(subterm2, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
748  pmf_coeff, f_coeffs, alpha, WT, gamma, VT);
749 
750  const std::string aux_fems = pmf_coeff ? "#1,#2,#3,#4" : "#1,#2,#3";
751 
752  getfem::generic_assembly assem;
753  assem.set("V$1(#1)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1))(i,:,i); "
754  "V$2(#3)+=comp(NonLin$2(#1," + aux_fems + ").vBase(#3))(i,:,i)");
755  assem.push_mi(mim);
756  assem.push_mf(mf_u);
757  assem.push_mf(mf_obs);
758  assem.push_mf(mf_lambda);
759  if (pmf_coeff)
760  assem.push_mf(*pmf_coeff);
761  assem.push_nonlinear_term(&nterm1);
762  assem.push_nonlinear_term(&nterm2);
763  assem.push_vec(Ru);
764  assem.push_vec(Rl);
765  assem.assembly(rg);
766  }
767 
768  struct integral_contact_rigid_obstacle_brick : public virtual_brick {
769 
770  bool contact_only;
771  int option;
772 
773  // option = 1 : Alart-Curnier
774  // option = 2 : symmetric Alart-Curnier (with friction, almost symmetric),
775  // option = 3 : Unsymmetric method based on augmented multipliers
776  // option = 4 : Unsymmetric method based on augmented multipliers
777  // with De-Saxce projection.
778 
779  virtual void asm_real_tangent_terms(const model &md, size_type /* ib */,
780  const model::varnamelist &vl,
781  const model::varnamelist &dl,
782  const model::mimlist &mims,
783  model::real_matlist &matl,
784  model::real_veclist &vecl,
785  model::real_veclist &,
786  size_type region,
787  build_version version) const {
788  GMM_ASSERT1(mims.size() == 1,
789  "Integral contact with rigid obstacle bricks need a single mesh_im");
790  GMM_ASSERT1(vl.size() == 2,
791  "Integral contact with rigid obstacle bricks need two variables");
792  GMM_ASSERT1(dl.size() >= 2 && dl.size() <= 7,
793  "Wrong number of data for integral contact with rigid obstacle "
794  << "brick, " << dl.size() << " should be between 2 and 7.");
795  GMM_ASSERT1(matl.size() == size_type(3 + (option == 2 && !contact_only)),
796  "Wrong number of terms for "
797  "integral contact with rigid obstacle brick");
798 
799  // variables : u, lambda. The variable lambda should be scalar in the
800  // frictionless case and vector valued in the case with
801  // friction.
802  // data : obstacle, r for the version without friction
803  // : obstacle, r, friction_coeffs, alpha, w_t, gamma, v_t for
804  // the version with friction. alpha, w_t , gamma and v_t
805  // are optional and equal to 1, 0, 1 and 0 by default,
806  // respectively.
807 
808  const model_real_plain_vector &u = md.real_variable(vl[0]);
809  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
810  const model_real_plain_vector &lambda = md.real_variable(vl[1]);
811  const mesh_fem &mf_lambda = md.mesh_fem_of_variable(vl[1]);
812  GMM_ASSERT1(mf_lambda.get_qdim() == (contact_only ? 1 : mf_u.get_qdim()),
813  "The contact stress has not the right dimension");
814  const model_real_plain_vector &obstacle = md.real_variable(dl[0]);
815  const mesh_fem &mf_obstacle = md.mesh_fem_of_variable(dl[0]);
816  size_type sl = gmm::vect_size(obstacle) * mf_obstacle.get_qdim()
817  / mf_obstacle.nb_dof();
818  GMM_ASSERT1(sl == 1, "the data corresponding to the obstacle has not "
819  "the right format");
820 
821  const model_real_plain_vector &vr = md.real_variable(dl[1]);
822  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Parameter r should be a scalar");
823  const mesh_im &mim = *mims[0];
824 
825  const model_real_plain_vector dummy_vec(0);
826  const model_real_plain_vector &friction_coeffs = contact_only
827  ? dummy_vec : md.real_variable(dl[2]);
828  const mesh_fem *pmf_coeff = contact_only ? 0 : md.pmesh_fem_of_variable(dl[2]);
829  sl = gmm::vect_size(friction_coeffs);
830  if (pmf_coeff) { sl *= pmf_coeff->get_qdim(); sl /= pmf_coeff->nb_dof(); }
831  GMM_ASSERT1(sl == 1 || sl == 2 || sl == 3 || contact_only,
832  "the data corresponding to the friction coefficient "
833  "has not the right format");
834 
835  scalar_type alpha = 1;
836  if (!contact_only && dl.size() >= 4) {
837  alpha = md.real_variable(dl[3])[0];
838  GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[3])) == 1,
839  "Parameter alpha should be a scalar");
840  }
841 
842  const model_real_plain_vector *WT = 0;
843  if (!contact_only && dl.size() >= 5) {
844  if (dl[4].compare(vl[0]) != 0)
845  WT = &(md.real_variable(dl[4]));
846  else if (md.n_iter_of_variable(vl[0]) > 1)
847  WT = &(md.real_variable(vl[0],1));
848  }
849 
850  scalar_type gamma = 1;
851  if (!contact_only && dl.size() >= 6) {
852  GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[5])) == 1,
853  "Parameter gamma should be a scalar");
854  gamma = md.real_variable(dl[5])[0];
855  }
856 
857  const model_real_plain_vector *VT
858  = (!contact_only && dl.size()>=7) ? &(md.real_variable(dl[6])) : 0;
859 
860  mesh_region rg(region);
861  mf_u.linked_mesh().intersect_with_mpi_region(rg);
862 
863  if (version & model::BUILD_MATRIX) {
864  GMM_TRACE2("Integral contact with rigid obstacle friction tangent term");
865  gmm::clear(matl[0]); gmm::clear(matl[1]); gmm::clear(matl[2]);
866  if (matl.size() >= 4) gmm::clear(matl[3]);
867  size_type fourthmat = (matl.size() >= 4) ? 3 : 1;
868  if (contact_only)
869  asm_Alart_Curnier_contact_rigid_obstacle_tangent_matrix
870  (matl[0], matl[1], matl[2], matl[fourthmat], mim,
871  mf_u, u, mf_obstacle, obstacle, mf_lambda, lambda, vr[0],
872  rg, option);
873  else
874  asm_Alart_Curnier_contact_rigid_obstacle_tangent_matrix
875  (matl[0], matl[1], matl[2], matl[fourthmat], mim,
876  mf_u, u, mf_obstacle, obstacle, mf_lambda, lambda,
877  pmf_coeff, &friction_coeffs, vr[0], alpha, WT, gamma, VT,
878  rg, option);
879  }
880 
881  if (version & model::BUILD_RHS) {
882  gmm::clear(vecl[0]); gmm::clear(vecl[1]); gmm::clear(vecl[2]);
883  if (matl.size() >= 4) gmm::clear(vecl[3]);
884 
885  if (contact_only)
886  asm_Alart_Curnier_contact_rigid_obstacle_rhs
887  (vecl[0], vecl[2], mim,
888  mf_u, u, mf_obstacle, obstacle, mf_lambda, lambda, vr[0],
889  rg, option);
890  else
891  asm_Alart_Curnier_contact_rigid_obstacle_rhs
892  (vecl[0], vecl[2], mim,
893  mf_u, u, mf_obstacle, obstacle, mf_lambda, lambda,
894  pmf_coeff, &friction_coeffs, vr[0], alpha, WT, gamma, VT,
895  rg, option);
896  }
897 
898  }
899 
900  integral_contact_rigid_obstacle_brick(bool contact_only_, int option_) {
901  option = option_;
902  contact_only = contact_only_;
903  set_flags(contact_only
904  ? "Integral contact with rigid obstacle brick"
905  : "Integral contact and friction with rigid obstacle brick",
906  false /* is linear*/,
907  (option==2) && contact_only /* is symmetric */,
908  false /* is coercive */,
909  true /* is real */, false /* is complex */);
910  }
911 
912  };
913 
914 
915  //=========================================================================
916  // Add a frictionless contact condition with a rigid obstacle given
917  // by a level set.
918  //=========================================================================
919 
920  size_type add_integral_contact_with_rigid_obstacle_brick
921  (model &md, const mesh_im &mim, const std::string &varname_u,
922  const std::string &multname_n, const std::string &dataname_obs,
923  const std::string &dataname_r, size_type region, int option) {
924 
925  pbrick pbr
926  = std::make_shared<integral_contact_rigid_obstacle_brick>(true, option);
927 
928  model::termlist tl;
929 
930  switch (option) {
931  case 1 : case 3 :
932  tl.push_back(model::term_description(varname_u, multname_n, false)); // UL
933  tl.push_back(model::term_description(multname_n, varname_u, false)); // LU
934  tl.push_back(model::term_description(multname_n, multname_n, true)); // LL
935  break;
936  case 2 :
937  tl.push_back(model::term_description(varname_u, multname_n, true)); // UL
938  tl.push_back(model::term_description(varname_u, varname_u, true)); // UU (fourthmat == 1)
939  tl.push_back(model::term_description(multname_n, multname_n, true)); // LL
940  break;
941  default :GMM_ASSERT1(false,
942  "Incorrect option for integral contact brick");
943 
944  }
945  model::varnamelist dl(1, dataname_obs);
946  dl.push_back(dataname_r);
947 
948  model::varnamelist vl(1, varname_u);
949  vl.push_back(multname_n);
950 
951  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
952  }
953 
954 
955  //=========================================================================
956  // Add a contact condition with Coulomb friction with a rigid obstacle
957  // given by a level set.
958  //=========================================================================
959 
960  size_type add_integral_contact_with_rigid_obstacle_brick
961  (model &md, const mesh_im &mim, const std::string &varname_u,
962  const std::string &multname, const std::string &dataname_obs,
963  const std::string &dataname_r, const std::string &dataname_friction_coeffs,
964  size_type region, int option,
965  const std::string &dataname_alpha, const std::string &dataname_wt,
966  const std::string &dataname_gamma, const std::string &dataname_vt) {
967 
968  pbrick pbr = std::make_shared<integral_contact_rigid_obstacle_brick>
969  (false, option);
970 
971  model::termlist tl;
972 
973  switch (option) {
974  case 1: case 3: case 4:
975  tl.push_back(model::term_description(varname_u, multname, false)); // UL
976  tl.push_back(model::term_description(multname, varname_u, false)); // LU
977  tl.push_back(model::term_description(multname, multname, true)); // LL
978  break;
979  case 2:
980  tl.push_back(model::term_description(varname_u, multname, false)); // UL
981  tl.push_back(model::term_description(multname, varname_u, false)); // LU
982  tl.push_back(model::term_description(multname, multname, true)); // LL
983  tl.push_back(model::term_description(varname_u, varname_u, true)); // UU
984  break;
985  default :GMM_ASSERT1(false,
986  "Incorrect option for integral contact brick");
987  }
988  model::varnamelist dl(1, dataname_obs);
989  dl.push_back(dataname_r);
990  dl.push_back(dataname_friction_coeffs);
991  if (dataname_alpha.size()) {
992  dl.push_back(dataname_alpha);
993  if (dataname_wt.size()) {
994  dl.push_back(dataname_wt);
995  if (dataname_gamma.size()) {
996  dl.push_back(dataname_gamma);
997  if (dataname_vt.size()) dl.push_back(dataname_vt);
998  }
999  }
1000  }
1001 
1002  model::varnamelist vl(1, varname_u);
1003  vl.push_back(multname);
1004 
1005  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1006  }
1007 
1008 
1009  //=========================================================================
1010  //
1011  // Integral penalized contact with friction (given obstacle, u, lambda).
1012  //
1013  //=========================================================================
1014 
1015  template<typename MAT, typename VECT1>
1016  void asm_penalized_contact_rigid_obstacle_tangent_matrix // frictionless
1017  (MAT &Kuu,
1018  const mesh_im &mim,
1019  const getfem::mesh_fem &mf_u, const VECT1 &U,
1020  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
1021  const getfem::mesh_fem *pmf_lambda, const VECT1 *lambda,
1022  scalar_type r, const mesh_region &rg, int option = 1) {
1023 
1024  contact_rigid_obstacle_nonlinear_term
1025  nterm((option == 1) ? K_UU_V1 : K_UU_V2, r,
1026  mf_u, U, mf_obs, obs, pmf_lambda, lambda);
1027 
1028  const std::string aux_fems = pmf_lambda ? "#1,#2,#3": "#1,#2";
1029  getfem::generic_assembly assem;
1030  assem.set("M(#1,#1)+=comp(NonLin(#1," + aux_fems + ").vBase(#1).vBase(#1))(i,j,:,i,:,j)");
1031  assem.push_mi(mim);
1032  assem.push_mf(mf_u);
1033  assem.push_mf(mf_obs);
1034  if (pmf_lambda)
1035  assem.push_mf(*pmf_lambda);
1036  assem.push_nonlinear_term(&nterm);
1037  assem.push_mat(Kuu);
1038  assem.assembly(rg);
1039  }
1040 
1041 
1042  template<typename VECT1>
1043  void asm_penalized_contact_rigid_obstacle_rhs // frictionless
1044  (VECT1 &Ru,
1045  const mesh_im &mim,
1046  const getfem::mesh_fem &mf_u, const VECT1 &U,
1047  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
1048  const getfem::mesh_fem *pmf_lambda, const VECT1 *lambda,
1049  scalar_type r, const mesh_region &rg, int option = 1) {
1050 
1051  contact_rigid_obstacle_nonlinear_term
1052  nterm((option == 1) ? RHS_U_V5 : RHS_U_V2, r,
1053  mf_u, U, mf_obs, obs, pmf_lambda, lambda);
1054 
1055  const std::string aux_fems = pmf_lambda ? "#1,#2,#3": "#1,#2";
1056  getfem::generic_assembly assem;
1057  assem.set("V(#1)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1))(i,:,i); ");
1058  assem.push_mi(mim);
1059  assem.push_mf(mf_u);
1060  assem.push_mf(mf_obs);
1061  if (pmf_lambda)
1062  assem.push_mf(*pmf_lambda);
1063  assem.push_nonlinear_term(&nterm);
1064  assem.push_vec(Ru);
1065  assem.assembly(rg);
1066  }
1067 
1068  template<typename MAT, typename VECT1>
1069  void asm_penalized_contact_rigid_obstacle_tangent_matrix // with friction
1070  (MAT &Kuu,
1071  const mesh_im &mim,
1072  const getfem::mesh_fem &mf_u, const VECT1 &U,
1073  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
1074  const getfem::mesh_fem *pmf_lambda, const VECT1 *lambda,
1075  const getfem::mesh_fem *pmf_coeff, const VECT1 *f_coeffs, scalar_type r,
1076  scalar_type alpha, const VECT1 *WT,
1077  const mesh_region &rg, int option = 1) {
1078 
1079  size_type subterm = 0;
1080  switch (option) {
1081  case 1 : subterm = K_UU_FRICT_V4; break;
1082  case 2 : subterm = K_UU_FRICT_V3; break;
1083  case 3 : subterm = K_UU_FRICT_V5; break;
1084  }
1085 
1086  contact_rigid_obstacle_nonlinear_term
1087  nterm(subterm, r, mf_u, U, mf_obs, obs, pmf_lambda, lambda,
1088  pmf_coeff, f_coeffs, alpha, WT);
1089 
1090  const std::string aux_fems = pmf_coeff ? "#1,#2,#3,#4"
1091  : (pmf_lambda ? "#1,#2,#3": "#1,#2");
1092  getfem::generic_assembly assem;
1093  assem.set("M(#1,#1)+=comp(NonLin(#1," + aux_fems + ").vBase(#1).vBase(#1))(i,j,:,i,:,j)");
1094  assem.push_mi(mim);
1095  assem.push_mf(mf_u);
1096  assem.push_mf(mf_obs);
1097 
1098  if (pmf_lambda)
1099  assem.push_mf(*pmf_lambda);
1100  else if (pmf_coeff)
1101  assem.push_mf(*pmf_coeff); // dummy
1102 
1103  if (pmf_coeff)
1104  assem.push_mf(*pmf_coeff);
1105 
1106  assem.push_nonlinear_term(&nterm);
1107  assem.push_mat(Kuu);
1108  assem.assembly(rg);
1109  }
1110 
1111 
1112  template<typename VECT1>
1113  void asm_penalized_contact_rigid_obstacle_rhs // with friction
1114  (VECT1 &Ru,
1115  const mesh_im &mim,
1116  const getfem::mesh_fem &mf_u, const VECT1 &U,
1117  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
1118  const getfem::mesh_fem *pmf_lambda, const VECT1 *lambda,
1119  const getfem::mesh_fem *pmf_coeff, const VECT1 *f_coeffs, scalar_type r,
1120  scalar_type alpha, const VECT1 *WT,
1121  const mesh_region &rg, int option = 1) {
1122 
1123  size_type subterm = 0;
1124  switch (option) {
1125  case 1 : subterm = RHS_U_FRICT_V7; break;
1126  case 2 : subterm = RHS_U_FRICT_V6; break;
1127  case 3 : subterm = RHS_U_FRICT_V8; break;
1128  }
1129 
1130  contact_rigid_obstacle_nonlinear_term
1131  nterm(subterm, r, mf_u, U, mf_obs, obs, pmf_lambda, lambda,
1132  pmf_coeff, f_coeffs, alpha, WT);
1133 
1134  const std::string aux_fems = pmf_coeff ? "#1,#2,#3,#4"
1135  : (pmf_lambda ? "#1,#2,#3": "#1,#2");
1136  getfem::generic_assembly assem;
1137  assem.set("V(#1)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1))(i,:,i); ");
1138  assem.push_mi(mim);
1139  assem.push_mf(mf_u);
1140  assem.push_mf(mf_obs);
1141 
1142  if (pmf_lambda)
1143  assem.push_mf(*pmf_lambda);
1144  else if (pmf_coeff)
1145  assem.push_mf(*pmf_coeff); // dummy
1146 
1147  if (pmf_coeff)
1148  assem.push_mf(*pmf_coeff);
1149 
1150  assem.push_nonlinear_term(&nterm);
1151  assem.push_vec(Ru);
1152  assem.assembly(rg);
1153  }
1154 
1155 
1156  struct penalized_contact_rigid_obstacle_brick : public virtual_brick {
1157 
1158  bool contact_only;
1159  int option;
1160 
1161  virtual void asm_real_tangent_terms(const model &md, size_type /* ib */,
1162  const model::varnamelist &vl,
1163  const model::varnamelist &dl,
1164  const model::mimlist &mims,
1165  model::real_matlist &matl,
1166  model::real_veclist &vecl,
1167  model::real_veclist &,
1168  size_type region,
1169  build_version version) const {
1170  // Integration method
1171  GMM_ASSERT1(mims.size() == 1,
1172  "Penalized contact with rigid obstacle bricks need a single mesh_im");
1173  const mesh_im &mim = *mims[0];
1174 
1175  // Variables : u
1176  GMM_ASSERT1(vl.size() == 1,
1177  "Penalized contact with rigid obstacle bricks need a single variable");
1178  const model_real_plain_vector &u = md.real_variable(vl[0]);
1179  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
1180 
1181  size_type N = mf_u.linked_mesh().dim();
1182 
1183  // Data : obs, r, [lambda,] [friction_coeffs,] [alpha,] [WT]
1184  size_type nb_data_1 = ((option == 1) ? 2 : 3) + (contact_only ? 0 : 1);
1185  size_type nb_data_2 = nb_data_1 + (contact_only ? 0 : 2);
1186  GMM_ASSERT1(dl.size() >= nb_data_1 && dl.size() <= nb_data_2,
1187  "Wrong number of data for penalized contact with rigid obstacle "
1188  << "brick, " << dl.size() << " should be between "
1189  << nb_data_1 << " and " << nb_data_2 << ".");
1190 
1191  size_type nd = 0;
1192  const model_real_plain_vector &obs = md.real_variable(dl[nd]);
1193  const mesh_fem &mf_obs = md.mesh_fem_of_variable(dl[nd]);
1194  size_type sl = gmm::vect_size(obs) * mf_obs.get_qdim() / mf_obs.nb_dof();
1195  GMM_ASSERT1(sl == 1, "the data corresponding to the obstacle has not "
1196  "the right format");
1197 
1198  nd++;
1199  const model_real_plain_vector &vr = md.real_variable(dl[nd]);
1200  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Parameter r should be a scalar");
1201 
1202  const model_real_plain_vector *lambda = 0;
1203  const mesh_fem *pmf_lambda = 0;
1204  if (option != 1) {
1205  nd++;
1206  lambda = &(md.real_variable(dl[nd]));
1207  pmf_lambda = md.pmesh_fem_of_variable(dl[nd]);
1208  sl = gmm::vect_size(*lambda) * pmf_lambda->get_qdim() / pmf_lambda->nb_dof();
1209  GMM_ASSERT1(sl == (contact_only ? 1 : N),
1210  "the data corresponding to the contact stress "
1211  "has not the right format");
1212  }
1213 
1214  const model_real_plain_vector *f_coeffs = 0;
1215  const mesh_fem *pmf_coeff = 0;
1216  scalar_type alpha = 1;
1217  const model_real_plain_vector *WT = 0;
1218  if (!contact_only) {
1219  nd++;
1220  f_coeffs = &(md.real_variable(dl[nd]));
1221  pmf_coeff = md.pmesh_fem_of_variable(dl[nd]);
1222  sl = gmm::vect_size(*f_coeffs);
1223  if (pmf_coeff) { sl *= pmf_coeff->get_qdim(); sl /= pmf_coeff->nb_dof(); }
1224  GMM_ASSERT1(sl == 1 || sl == 2 || sl == 3,
1225  "the data corresponding to the friction coefficient "
1226  "has not the right format");
1227  if (dl.size() > nd+1) {
1228  nd++;
1229  alpha = md.real_variable(dl[nd])[0];
1230  GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[nd])) == 1,
1231  "Parameter alpha should be a scalar");
1232  }
1233 
1234  if (dl.size() > nd+1) {
1235  nd++;
1236  if (dl[nd].compare(vl[0]) != 0)
1237  WT = &(md.real_variable(dl[nd]));
1238  else if (md.n_iter_of_variable(vl[0]) > 1)
1239  WT = &(md.real_variable(vl[0],1));
1240  }
1241  }
1242 
1243  GMM_ASSERT1(matl.size() == 1, "Wrong number of terms for "
1244  "penalized contact with rigid obstacle brick");
1245 
1246  mesh_region rg(region);
1247  mf_u.linked_mesh().intersect_with_mpi_region(rg);
1248 
1249  if (version & model::BUILD_MATRIX) {
1250  GMM_TRACE2("Penalized contact with rigid obstacle tangent term");
1251  gmm::clear(matl[0]);
1252  if (contact_only)
1253  asm_penalized_contact_rigid_obstacle_tangent_matrix
1254  (matl[0], mim, mf_u, u, mf_obs, obs, pmf_lambda, lambda,
1255  vr[0], rg, option);
1256  else
1257  asm_penalized_contact_rigid_obstacle_tangent_matrix
1258  (matl[0], mim, mf_u, u, mf_obs, obs, pmf_lambda, lambda,
1259  pmf_coeff, f_coeffs, vr[0], alpha, WT, rg, option);
1260  }
1261 
1262  if (version & model::BUILD_RHS) {
1263  gmm::clear(vecl[0]);
1264  if (contact_only)
1265  asm_penalized_contact_rigid_obstacle_rhs
1266  (vecl[0], mim, mf_u, u, mf_obs, obs, pmf_lambda, lambda,
1267  vr[0], rg, option);
1268  else
1269  asm_penalized_contact_rigid_obstacle_rhs
1270  (vecl[0], mim, mf_u, u, mf_obs, obs, pmf_lambda, lambda,
1271  pmf_coeff, f_coeffs, vr[0], alpha, WT, rg, option);
1272  }
1273 
1274  }
1275 
1276  penalized_contact_rigid_obstacle_brick(bool contact_only_, int option_) {
1277  contact_only = contact_only_;
1278  option = option_;
1279  set_flags(contact_only
1280  ? "Integral penalized contact with rigid obstacle brick"
1281  : "Integral penalized contact and friction with rigid obstacle brick",
1282  false /* is linear*/, contact_only /* is symmetric */,
1283  true /* is coercive */, true /* is real */,
1284  false /* is complex */);
1285  }
1286 
1287  };
1288 
1289 
1290  //=========================================================================
1291  // Add a frictionless contact condition with a rigid obstacle given
1292  // by a level set.
1293  //=========================================================================
1294 
1295  size_type add_penalized_contact_with_rigid_obstacle_brick
1296  (model &md, const mesh_im &mim, const std::string &varname_u,
1297  const std::string &dataname_obs, const std::string &dataname_r,
1298  size_type region, int option, const std::string &dataname_n) {
1299 
1300  pbrick pbr = std::make_shared<penalized_contact_rigid_obstacle_brick>
1301  (true, option);
1302 
1303  model::termlist tl;
1304  tl.push_back(model::term_description(varname_u, varname_u, true));
1305 
1306  model::varnamelist dl(1, dataname_obs);
1307  dl.push_back(dataname_r);
1308  switch (option) {
1309  case 1: break;
1310  case 2: dl.push_back(dataname_n); break;
1311  default: GMM_ASSERT1(false, "Penalized contact brick : invalid option");
1312  }
1313 
1314  model::varnamelist vl(1, varname_u);
1315 
1316  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1317  }
1318 
1319  //=========================================================================
1320  // Add a contact condition with friction with a rigid obstacle given
1321  // by a level set.
1322  //=========================================================================
1323 
1324  size_type add_penalized_contact_with_rigid_obstacle_brick
1325  (model &md, const mesh_im &mim, const std::string &varname_u,
1326  const std::string &dataname_obs, const std::string &dataname_r,
1327  const std::string &dataname_friction_coeffs,
1328  size_type region, int option, const std::string &dataname_lambda,
1329  const std::string &dataname_alpha, const std::string &dataname_wt) {
1330 
1331  pbrick pbr = std::make_shared<penalized_contact_rigid_obstacle_brick>
1332  (false, option);
1333 
1334  model::termlist tl;
1335  tl.push_back(model::term_description(varname_u, varname_u, false));
1336 
1337  model::varnamelist dl(1, dataname_obs);
1338  dl.push_back(dataname_r);
1339  switch (option) {
1340  case 1: break;
1341  case 2: case 3: dl.push_back(dataname_lambda); break;
1342  default: GMM_ASSERT1(false, "Penalized contact brick : invalid option");
1343  }
1344  dl.push_back(dataname_friction_coeffs);
1345  if (dataname_alpha.size() > 0) {
1346  dl.push_back(dataname_alpha);
1347  if (dataname_wt.size() > 0) dl.push_back(dataname_wt);
1348  }
1349 
1350  model::varnamelist vl(1, varname_u);
1351 
1352  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1353  }
1354 
1355 
1356  //=========================================================================
1357  //
1358  // Integral contact (with friction) between non-matching meshes.
1359  //
1360  //=========================================================================
1361 
1362  template<typename MAT, typename VEC>
1363  void asm_Alart_Curnier_contact_nonmatching_meshes_tangent_matrix // frictionless
1364  (MAT &Ku1l, MAT &Klu1, MAT &Ku2l, MAT &Klu2, MAT &Kll,
1365  MAT &Ku1u1, MAT &Ku2u2, MAT &Ku1u2,
1366  const mesh_im &mim,
1367  const getfem::mesh_fem &mf_u1, const VEC &U1,
1368  const getfem::mesh_fem &mf_u2, const VEC &U2,
1369  const getfem::mesh_fem &mf_lambda, const VEC &lambda,
1370  scalar_type r, const mesh_region &rg, int option = 1) {
1371 
1372  size_type subterm1 = (option == 3) ? K_UL_V2 : K_UL_V1;
1373  size_type subterm2 = (option == 3) ? K_UL_V1 : K_UL_V3;
1374  size_type subterm3 = (option == 3) ? K_LL_V2 : K_LL_V1;
1375  size_type subterm4 = (option == 2) ? K_UU_V2 : K_UU_V1;
1376 
1377  contact_nonmatching_meshes_nonlinear_term
1378  nterm1(subterm1, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda),
1379  nterm2(subterm2, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda),
1380  nterm3(subterm3, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda),
1381  nterm4(subterm4, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda);
1382 
1383  getfem::generic_assembly assem;
1384  switch (option) {
1385  case 1: case 3:
1386  assem.set
1387  ("M$1(#1,#3)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#1).Base(#3))(i,:,i,:); " // U1L
1388  "M$2(#3,#1)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3).vBase(#1))(i,:,:,i); " // LU1
1389  "M$3(#2,#3)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#2).Base(#3))(i,:,i,:); " // U2L
1390  "M$4(#3,#2)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3).vBase(#2))(i,:,:,i); " // LU2
1391  "M$5(#3,#3)+=comp(NonLin$3(#1,#1,#2,#3).Base(#3).Base(#3))(i,:,:)"); // LL
1392  break;
1393  case 2:
1394  assem.set
1395  ("M$1(#1,#3)+=comp(NonLin$2(#1,#1,#2,#3).vBase(#1).Base(#3))(i,:,i,:); " // U1L
1396  "M$3(#2,#3)+=comp(NonLin$2(#1,#1,#2,#3).vBase(#2).Base(#3))(i,:,i,:); " // U2L
1397  "M$5(#3,#3)+=comp(NonLin$3(#1,#1,#2,#3).Base(#3).Base(#3))(i,:,:); " // LL
1398  "M$6(#1,#1)+=comp(NonLin$4(#1,#1,#2,#3).vBase(#1).vBase(#1))(i,j,:,i,:,j); " // U1U1
1399  "M$7(#2,#2)+=comp(NonLin$4(#1,#1,#2,#3).vBase(#2).vBase(#2))(i,j,:,i,:,j); " // U2U2
1400  "M$8(#1,#2)+=comp(NonLin$4(#1,#1,#2,#3).vBase(#1).vBase(#2))(i,j,:,i,:,j)"); // U1U2
1401  break;
1402  }
1403  assem.push_mi(mim);
1404  assem.push_mf(mf_u1);
1405  assem.push_mf(mf_u2);
1406  assem.push_mf(mf_lambda);
1407  assem.push_nonlinear_term(&nterm1);
1408  assem.push_nonlinear_term(&nterm2);
1409  assem.push_nonlinear_term(&nterm3);
1410  assem.push_nonlinear_term(&nterm4);
1411  assem.push_mat(Ku1l);
1412  assem.push_mat(Klu1);
1413  assem.push_mat(Ku2l);
1414  assem.push_mat(Klu2);
1415  assem.push_mat(Kll);
1416  assem.push_mat(Ku1u1);
1417  assem.push_mat(Ku2u2);
1418  assem.push_mat(Ku1u2);
1419  assem.assembly(rg);
1420 
1421  gmm::scale(Ku2l, scalar_type(-1));
1422  if (option != 2) // Klu2 was calculated
1423  gmm::scale(Klu2, scalar_type(-1));
1424  gmm::scale(Ku1u2, scalar_type(-1));
1425  }
1426 
1427  template<typename MAT, typename VEC>
1428  void asm_Alart_Curnier_contact_nonmatching_meshes_tangent_matrix // with friction
1429  (MAT &Ku1l, MAT &Klu1, MAT &Ku2l, MAT &Klu2, MAT &Kll,
1430  MAT &Ku1u1, MAT &Ku2u2, MAT &Ku1u2, MAT &Ku2u1,
1431  const mesh_im &mim,
1432  const getfem::mesh_fem &mf_u1, const VEC &U1,
1433  const getfem::mesh_fem &mf_u2, const VEC &U2,
1434  const getfem::mesh_fem &mf_lambda, const VEC &lambda,
1435  const getfem::mesh_fem *pmf_coeff, const VEC *f_coeffs,
1436  scalar_type r, scalar_type alpha,
1437  const VEC *WT1, const VEC *WT2,
1438  const mesh_region &rg, int option = 1) {
1439 
1440  size_type subterm1, subterm2, subterm3;
1441  switch (option) {
1442  case 1 :
1443  subterm1 = K_UL_FRICT_V1; subterm2 = K_UL_FRICT_V4; subterm3 = K_LL_FRICT_V1;
1444  break;
1445  case 2 :
1446  subterm1 = K_UL_FRICT_V3; subterm2 = K_UL_FRICT_V4; subterm3 = K_LL_FRICT_V1;
1447  break;
1448  case 3 :
1449  subterm1 = K_UL_FRICT_V2; subterm2 = K_UL_FRICT_V5; subterm3 = K_LL_FRICT_V2;
1450  break;
1451  case 4 :
1452  subterm1 = K_UL_FRICT_V7; subterm2 = K_UL_FRICT_V8; subterm3 = K_LL_FRICT_V4;
1453  break;
1454  default : GMM_ASSERT1(false, "Incorrect option");
1455  }
1456 
1457  size_type subterm4 = K_UU_FRICT_V3;
1458 
1459  contact_nonmatching_meshes_nonlinear_term
1460  nterm1(subterm1, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1461  pmf_coeff, f_coeffs, alpha, WT1, WT2),
1462  nterm2(subterm2, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1463  pmf_coeff, f_coeffs, alpha, WT1, WT2),
1464  nterm3(subterm3, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1465  pmf_coeff, f_coeffs, alpha, WT1, WT2),
1466  nterm4(subterm4, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1467  pmf_coeff, f_coeffs, alpha, WT1, WT2);
1468 
1469  const std::string aux_fems = pmf_coeff ? "#1,#2,#3,#4" : "#1,#2,#3";
1470 
1471  getfem::generic_assembly assem;
1472  switch (option) {
1473  case 1: case 3: case 4:
1474  assem.set
1475  ("M$1(#1,#3)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1).vBase(#3))(i,j,:,i,:,j); " // U1L
1476  "M$2(#3,#1)+=comp(NonLin$2(#1," + aux_fems + ").vBase(#3).vBase(#1))(i,j,:,j,:,i); " // LU1
1477  "M$3(#2,#3)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#2).vBase(#3))(i,j,:,i,:,j); " // U2L
1478  "M$4(#3,#2)+=comp(NonLin$2(#1," + aux_fems + ").vBase(#3).vBase(#2))(i,j,:,j,:,i); " // LU2
1479  "M$5(#3,#3)+=comp(NonLin$3(#1," + aux_fems + ").vBase(#3).vBase(#3))(i,j,:,i,:,j)"); // LL
1480  break;
1481  case 2:
1482  assem.set
1483  ("M$1(#1,#3)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1).vBase(#3))(i,j,:,i,:,j); " // U1L
1484  "M$2(#3,#1)+=comp(NonLin$2(#1," + aux_fems + ").vBase(#3).vBase(#1))(i,j,:,j,:,i); " // LU1
1485  "M$3(#2,#3)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#2).vBase(#3))(i,j,:,i,:,j); " // U2L
1486  "M$4(#3,#2)+=comp(NonLin$2(#1," + aux_fems + ").vBase(#3).vBase(#2))(i,j,:,j,:,i); " // LU2
1487  "M$5(#3,#3)+=comp(NonLin$3(#1," + aux_fems + ").vBase(#3).vBase(#3))(i,j,:,i,:,j); " // LL
1488  "M$6(#1,#1)+=comp(NonLin$4(#1," + aux_fems + ").vBase(#1).vBase(#1))(i,j,:,i,:,j); " // U1U1
1489  "M$7(#2,#2)+=comp(NonLin$4(#1," + aux_fems + ").vBase(#2).vBase(#2))(i,j,:,i,:,j); " // U2U2
1490  "M$8(#1,#2)+=comp(NonLin$4(#1," + aux_fems + ").vBase(#1).vBase(#2))(i,j,:,i,:,j); " // U1U2
1491  "M$9(#2,#1)+=comp(NonLin$4(#1," + aux_fems + ").vBase(#2).vBase(#1))(i,j,:,i,:,j)"); // U2U1
1492  break;
1493  }
1494  assem.push_mi(mim);
1495  assem.push_mf(mf_u1);
1496  assem.push_mf(mf_u2);
1497  assem.push_mf(mf_lambda);
1498  if (pmf_coeff)
1499  assem.push_mf(*pmf_coeff);
1500  assem.push_nonlinear_term(&nterm1);
1501  assem.push_nonlinear_term(&nterm2);
1502  assem.push_nonlinear_term(&nterm3);
1503  assem.push_nonlinear_term(&nterm4);
1504  assem.push_mat(Ku1l);
1505  assem.push_mat(Klu1);
1506  assem.push_mat(Ku2l);
1507  assem.push_mat(Klu2);
1508  assem.push_mat(Kll);
1509  assem.push_mat(Ku1u1);
1510  assem.push_mat(Ku2u2);
1511  assem.push_mat(Ku1u2);
1512  assem.push_mat(Ku2u1);
1513  assem.assembly(rg);
1514 
1515  gmm::scale(Ku2l, scalar_type(-1));
1516  gmm::scale(Klu2, scalar_type(-1));
1517  gmm::scale(Ku1u2, scalar_type(-1));
1518  }
1519 
1520  template<typename VECT1>
1521  void asm_Alart_Curnier_contact_nonmatching_meshes_rhs // frictionless
1522  (VECT1 &Ru1, VECT1 &Ru2, VECT1 &Rl,
1523  const mesh_im &mim,
1524  const getfem::mesh_fem &mf_u1, const VECT1 &U1,
1525  const getfem::mesh_fem &mf_u2, const VECT1 &U2,
1526  const getfem::mesh_fem &mf_lambda, const VECT1 &lambda,
1527  scalar_type r, const mesh_region &rg, int option = 1) {
1528 
1529  size_type subterm1;
1530  switch (option) {
1531  case 1 : subterm1 = RHS_U_V1; break;
1532  case 2 : subterm1 = RHS_U_V2; break;
1533  case 3 : subterm1 = RHS_U_V4; break;
1534  default : GMM_ASSERT1(false, "Incorrect option");
1535  }
1536  size_type subterm2 = (option == 3) ? RHS_L_V2 : RHS_L_V1;
1537 
1538  contact_nonmatching_meshes_nonlinear_term
1539  nterm1(subterm1, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda),
1540  nterm2(subterm2, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda);
1541 
1542  getfem::generic_assembly assem;
1543  assem.set("V$1(#1)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#1))(i,:,i); "
1544  "V$2(#2)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#2))(i,:,i); "
1545  "V$3(#3)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3))(i,:)");
1546  assem.push_mi(mim);
1547  assem.push_mf(mf_u1);
1548  assem.push_mf(mf_u2);
1549  assem.push_mf(mf_lambda);
1550  assem.push_nonlinear_term(&nterm1);
1551  assem.push_nonlinear_term(&nterm2);
1552  assem.push_vec(Ru1);
1553  assem.push_vec(Ru2);
1554  assem.push_vec(Rl);
1555  assem.assembly(rg);
1556 
1557  gmm::scale(Ru2, scalar_type(-1));
1558  }
1559 
1560  template<typename VECT1>
1561  void asm_Alart_Curnier_contact_nonmatching_meshes_rhs // with friction
1562  (VECT1 &Ru1, VECT1 &Ru2, VECT1 &Rl,
1563  const mesh_im &mim,
1564  const getfem::mesh_fem &mf_u1, const VECT1 &U1,
1565  const getfem::mesh_fem &mf_u2, const VECT1 &U2,
1566  const getfem::mesh_fem &mf_lambda, const VECT1 &lambda,
1567  const getfem::mesh_fem *pmf_coeff, const VECT1 *f_coeffs,
1568  scalar_type r, scalar_type alpha,
1569  const VECT1 *WT1, const VECT1 *WT2,
1570  const mesh_region &rg, int option = 1) {
1571 
1572  size_type subterm1, subterm2;
1573  switch (option) {
1574  case 1 : subterm1 = RHS_U_FRICT_V1; subterm2 = RHS_L_FRICT_V1; break;
1575  case 2 : subterm1 = RHS_U_FRICT_V6; subterm2 = RHS_L_FRICT_V1; break;
1576  case 3 : subterm1 = RHS_U_FRICT_V4; subterm2 = RHS_L_FRICT_V2; break;
1577  case 4 : subterm1 = RHS_U_FRICT_V5; subterm2 = RHS_L_FRICT_V4; break;
1578  default : GMM_ASSERT1(false, "Incorrect option");
1579  }
1580 
1581  contact_nonmatching_meshes_nonlinear_term
1582  nterm1(subterm1, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1583  pmf_coeff, f_coeffs, alpha, WT1, WT2),
1584  nterm2(subterm2, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1585  pmf_coeff, f_coeffs, alpha, WT1, WT2);
1586 
1587  const std::string aux_fems = pmf_coeff ? "#1,#2,#3,#4" : "#1,#2,#3";
1588 
1589  getfem::generic_assembly assem;
1590  assem.set("V$1(#1)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1))(i,:,i); "
1591  "V$2(#2)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#2))(i,:,i); "
1592  "V$3(#3)+=comp(NonLin$2(#1," + aux_fems + ").vBase(#3))(i,:,i)");
1593  assem.push_mi(mim);
1594  assem.push_mf(mf_u1);
1595  assem.push_mf(mf_u2);
1596  assem.push_mf(mf_lambda);
1597  if (pmf_coeff)
1598  assem.push_mf(*pmf_coeff);
1599  assem.push_nonlinear_term(&nterm1);
1600  assem.push_nonlinear_term(&nterm2);
1601  assem.push_vec(Ru1);
1602  assem.push_vec(Ru2);
1603  assem.push_vec(Rl);
1604  assem.assembly(rg);
1605 
1606  gmm::scale(Ru2, scalar_type(-1));
1607  }
1608 
1609  struct integral_contact_nonmatching_meshes_brick : public virtual_brick {
1610 
1611  size_type rg1, rg2; // ids of mesh regions on mf_u1 and mf_u2 that are
1612  // expected to come in contact.
1613  mutable getfem::pfem pfem_proj; // cached fem for the projection between nonmatching meshes
1614  bool contact_only;
1615  int option;
1616 
1617  // option = 1 : Alart-Curnier
1618  // option = 2 : symmetric Alart-Curnier (almost symmetric with friction)
1619  // option = 3 : New method
1620  // option = 4 : New method based on De-Saxce.
1621 
1622  virtual void asm_real_tangent_terms(const model &md, size_type /* ib */,
1623  const model::varnamelist &vl,
1624  const model::varnamelist &dl,
1625  const model::mimlist &mims,
1626  model::real_matlist &matl,
1627  model::real_veclist &vecl,
1628  model::real_veclist &,
1629  size_type region,
1630  build_version version) const {
1631  // Integration method
1632  GMM_ASSERT1(mims.size() == 1,
1633  "Integral contact between nonmatching meshes bricks need a single mesh_im");
1634  const mesh_im &mim = *mims[0];
1635 
1636  // Variables : u1, u2, lambda
1637  // the variable lambda should be scalar in the frictionless
1638  // case and vector valued in the case with friction.
1639  GMM_ASSERT1(vl.size() == 3,
1640  "Integral contact between nonmatching meshes bricks need three variables");
1641  const model_real_plain_vector &u1 = md.real_variable(vl[0]);
1642  const model_real_plain_vector &u2 = md.real_variable(vl[1]);
1643  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(vl[0]);
1644  const mesh_fem &mf_u2 = md.mesh_fem_of_variable(vl[1]);
1645  const model_real_plain_vector &lambda = md.real_variable(vl[2]);
1646  const mesh_fem &mf_lambda = md.mesh_fem_of_variable(vl[2]);
1647  GMM_ASSERT1(mf_lambda.get_qdim() == (contact_only ? 1 : mf_u1.get_qdim()),
1648  "The contact stress variable has not the right dimension");
1649 
1650  // Data : r, [friction_coeffs,] [alpha,] [WT1, WT2]
1651  // alpha, WT1, WT2 are optional and equal to 1, 0 and 0 by default respectively.
1652  if (contact_only) {
1653  GMM_ASSERT1(dl.size() == 1,
1654  "Wrong number of data for integral contact between nonmatching meshes "
1655  << "brick, the number of data should be equal to 1 .");
1656  }
1657  else {
1658  GMM_ASSERT1(dl.size() >= 2 && dl.size() <= 5,
1659  "Wrong number of data for integral contact between nonmatching meshes "
1660  << "brick, it should be between 2 and 5 instead of " << dl.size() << " .");
1661  }
1662  const model_real_plain_vector &vr = md.real_variable(dl[0]);
1663  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Parameter r should be a scalar");
1664 
1665  const model_real_plain_vector *f_coeffs = 0, *WT1 = 0, *WT2 = 0;
1666  const mesh_fem *pmf_coeff = 0;
1667  scalar_type alpha = 1;
1668  if (!contact_only) {
1669  f_coeffs = &(md.real_variable(dl[1]));
1670  pmf_coeff = md.pmesh_fem_of_variable(dl[1]);
1671 
1672  size_type sl = gmm::vect_size(*f_coeffs);
1673  if (pmf_coeff) { sl *= pmf_coeff->get_qdim(); sl /= pmf_coeff->nb_dof(); }
1674  GMM_ASSERT1(sl == 1 || sl == 2 || sl ==3,
1675  "the data corresponding to the friction coefficient "
1676  "has not the right format");
1677 
1678  if (dl.size() >= 3) {
1679  alpha = md.real_variable(dl[2])[0];
1680  GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[2])) == 1,
1681  "Parameter alpha should be a scalar");
1682  }
1683 
1684  if (dl.size() >= 4) {
1685  if (dl[3].compare(vl[0]) != 0)
1686  WT1 = &(md.real_variable(dl[3]));
1687  else if (md.n_iter_of_variable(vl[0]) > 1)
1688  WT1 = &(md.real_variable(vl[0],1));
1689  }
1690 
1691  if (dl.size() >= 5) {
1692  if (dl[4].compare(vl[1]) != 0)
1693  WT2 = &(md.real_variable(dl[4]));
1694  else if (md.n_iter_of_variable(vl[1]) > 1)
1695  WT2 = &(md.real_variable(vl[1],1));
1696  }
1697  }
1698 
1699  // Matrix terms (T_u1l, T_lu1, T_u2l, T_lu2, T_ll, T_u1u1, T_u2u2, T_u1u2)
1700  GMM_ASSERT1(matl.size() == size_type(3 + // U1L, U2L, LL
1701  2 * !is_symmetric() + // LU1, LU2
1702  3 * (option == 2) + // U1U1, U2U2, U1U2
1703  1 * (option == 2 && !is_symmetric())), // U2U1
1704  "Wrong number of terms for "
1705  "integral contact between nonmatching meshes brick");
1706 
1707  mesh_region rg(region);
1708  mf_u1.linked_mesh().intersect_with_mpi_region(rg);
1709 
1710  size_type N = mf_u1.linked_mesh().dim();
1711 
1712  // projection of the second mesh_fem onto the mesh of the first mesh_fem
1713  if (!pfem_proj)
1714  pfem_proj = new_projected_fem(mf_u2, mim, rg2, rg1);
1715 
1716  getfem::mesh_fem mf_u2_proj(mim.linked_mesh(), dim_type(N));
1717  mf_u2_proj.set_finite_element(mim.linked_mesh().convex_index(), pfem_proj);
1718 
1719  size_type nbdof1 = mf_u1.nb_dof();
1720  size_type nbdof_lambda = mf_lambda.nb_dof();
1721  size_type nbdof2 = mf_u2.nb_dof();
1722  size_type nbsub = mf_u2_proj.nb_basic_dof();
1723 
1724  std::vector<size_type> ind;
1725  mf_u2_proj.get_global_dof_index(ind);
1726  gmm::unsorted_sub_index SUBI(ind);
1727 
1728  gmm::csc_matrix<scalar_type> Esub(nbsub, nbdof2);
1729  if (mf_u2.is_reduced())
1730  gmm::copy(gmm::sub_matrix(mf_u2.extension_matrix(),
1731  SUBI, gmm::sub_interval(0, nbdof2)),
1732  Esub);
1733 
1734  model_real_plain_vector u2_proj(nbsub);
1735  if (mf_u2.is_reduced())
1736  gmm::mult(Esub, u2, u2_proj);
1737  else
1738  gmm::copy(gmm::sub_vector(u2, SUBI), u2_proj);
1739 
1740  model_real_plain_vector WT2_proj(0);
1741  if (WT2) {
1742  gmm::resize(WT2_proj, nbsub);
1743  if (mf_u2.is_reduced())
1744  gmm::mult(Esub, *WT2, WT2_proj);
1745  else
1746  gmm::copy(gmm::sub_vector(*WT2, SUBI), WT2_proj);
1747  }
1748 
1749  size_type U1L = 0;
1750  size_type LU1 = is_symmetric() ? size_type(-1) : 1;
1751  size_type U2L = is_symmetric() ? 1 : 2;
1752  size_type LU2 = is_symmetric() ? size_type(-1) : 3;
1753  size_type LL = is_symmetric() ? 2 : 4;
1754  size_type U1U1 = (option != 2) ? size_type(-1) : (is_symmetric() ? 3 : 5);
1755  size_type U2U2 = (option != 2) ? size_type(-1) : (is_symmetric() ? 4 : 6);
1756  size_type U1U2 = (option != 2) ? size_type(-1) : (is_symmetric() ? 5 : 7);
1757  size_type U2U1 = (option != 2 || is_symmetric()) ? size_type(-1) : 8;
1758 
1759  if (version & model::BUILD_MATRIX) {
1760  GMM_TRACE2("Integral contact between nonmatching meshes "
1761  "tangent term");
1762  for (size_type i = 0; i < matl.size(); i++) gmm::clear(matl[i]);
1763 
1764  model_real_sparse_matrix dummy_mat(0, 0);
1765  model_real_sparse_matrix &Klu1 = (LU1 == size_type(-1)) ? dummy_mat : matl[LU1];
1766  model_real_sparse_matrix &Ku1u1 = (U1U1 == size_type(-1)) ? dummy_mat : matl[U1U1];
1767 
1768  model_real_sparse_matrix Ku2l(nbsub, nbdof_lambda);
1769  model_real_sparse_matrix Klu2(nbdof_lambda, nbsub);
1770  model_real_sparse_matrix Ku2u2(nbsub, nbsub);
1771  model_real_sparse_matrix Ku1u2(nbdof1, nbsub);
1772  model_real_sparse_matrix Ku2u1(nbsub, nbdof1);
1773 
1774  if (contact_only)
1775  asm_Alart_Curnier_contact_nonmatching_meshes_tangent_matrix
1776  (matl[U1L], Klu1, Ku2l, Klu2, matl[LL], Ku1u1, Ku2u2, Ku1u2,
1777  mim, mf_u1, u1, mf_u2_proj, u2_proj, mf_lambda, lambda,
1778  vr[0], rg, option);
1779  else
1780  asm_Alart_Curnier_contact_nonmatching_meshes_tangent_matrix
1781  (matl[U1L], Klu1, Ku2l, Klu2, matl[LL], Ku1u1, Ku2u2, Ku1u2, Ku2u1,
1782  mim, mf_u1, u1, mf_u2_proj, u2_proj, mf_lambda, lambda,
1783  pmf_coeff, f_coeffs, vr[0], alpha, WT1, &WT2_proj, rg, option);
1784 
1785  if (mf_u2.is_reduced()) {
1786  gmm::mult(gmm::transposed(Esub), Ku2l, matl[U2L]);
1787  if (LU2 != size_type(-1)) gmm::mult(Klu2, Esub, matl[LU2]);
1788  if (U2U2 != size_type(-1)) {
1789  model_real_sparse_matrix tmp(nbsub, nbdof2);
1790  gmm::mult(Ku2u2, Esub, tmp);
1791  gmm::mult(gmm::transposed(Esub), tmp, matl[U2U2]);
1792  }
1793  if (U1U2 != size_type(-1)) gmm::mult(Ku1u2, Esub, matl[U1U2]);
1794  if (U2U1 != size_type(-1)) gmm::mult(gmm::transposed(Esub), Ku2u1, matl[U2U1]);
1795  }
1796  else {
1797  gmm::copy(Ku2l, gmm::sub_matrix(matl[U2L], SUBI, gmm::sub_interval(0, nbdof_lambda)));
1798  if (LU2 != size_type(-1))
1799  gmm::copy(Klu2, gmm::sub_matrix(matl[LU2], gmm::sub_interval(0, nbdof_lambda), SUBI));
1800  if (U2U2 != size_type(-1))
1801  gmm::copy(Ku2u2, gmm::sub_matrix(matl[U2U2], SUBI));
1802  if (U1U2 != size_type(-1))
1803  gmm::copy(Ku1u2, gmm::sub_matrix(matl[U1U2], gmm::sub_interval(0, nbdof1), SUBI));
1804  if (U2U1 != size_type(-1))
1805  gmm::copy(Ku2u1, gmm::sub_matrix(matl[U2U1], SUBI, gmm::sub_interval(0, nbdof1)));
1806  }
1807  }
1808 
1809  if (version & model::BUILD_RHS) {
1810  for (size_type i = 0; i < matl.size(); i++) gmm::clear(vecl[i]);
1811 
1812  model_real_plain_vector Ru2(nbsub);
1813 
1814  if (contact_only)
1815  asm_Alart_Curnier_contact_nonmatching_meshes_rhs
1816  (vecl[U1L], Ru2, vecl[LL], // u1, u2, lambda
1817  mim, mf_u1, u1, mf_u2_proj, u2_proj, mf_lambda, lambda,
1818  vr[0], rg, option);
1819  else
1820  asm_Alart_Curnier_contact_nonmatching_meshes_rhs
1821  (vecl[U1L], Ru2, vecl[LL], // u1, u2, lambda
1822  mim, mf_u1, u1, mf_u2_proj, u2_proj, mf_lambda, lambda,
1823  pmf_coeff, f_coeffs, vr[0], alpha, WT1, &WT2_proj, rg, option);
1824 
1825  if (mf_u2.is_reduced())
1826  gmm::mult(gmm::transposed(Esub), Ru2, vecl[U2L]);
1827  else
1828  gmm::copy(Ru2, gmm::sub_vector(vecl[U2L], SUBI));
1829  }
1830  }
1831 
1832  integral_contact_nonmatching_meshes_brick(size_type rg1_, size_type rg2_,
1833  bool contact_only_, int option_)
1834  : rg1(rg1_), rg2(rg2_), pfem_proj(0),
1835  contact_only(contact_only_), option(option_)
1836  {
1837  set_flags(contact_only
1838  ? "Integral contact between nonmatching meshes brick"
1839  : "Integral contact and friction between nonmatching "
1840  "meshes brick",
1841  false /* is linear*/,
1842  (option==2) && contact_only /* is symmetric */,
1843  false /* is coercive */, true /* is real */,
1844  false /* is complex */);
1845  }
1846 
1847  ~integral_contact_nonmatching_meshes_brick()
1848  { if (pfem_proj) del_projected_fem(pfem_proj); }
1849 
1850  };
1851 
1852 
1853  //=========================================================================
1854  // Add a frictionless contact condition between two bodies discretized
1855  // with nonmatching meshes.
1856  //=========================================================================
1857 
1858  size_type add_integral_contact_between_nonmatching_meshes_brick
1859  (model &md, const mesh_im &mim, const std::string &varname_u1,
1860  const std::string &varname_u2, const std::string &multname_n,
1861  const std::string &dataname_r,
1862  size_type region1, size_type region2, int option) {
1863 
1864  pbrick pbr = std::make_shared<integral_contact_nonmatching_meshes_brick>
1865  (region1, region2, true /* contact_only */, option);
1866 
1867  model::termlist tl;
1868 
1869  switch (option) {
1870  case 1 : case 3 :
1871  tl.push_back(model::term_description(varname_u1, multname_n, false)); // U1L
1872  tl.push_back(model::term_description(multname_n, varname_u1, false)); // LU1
1873  tl.push_back(model::term_description(varname_u2, multname_n, false)); // U2L
1874  tl.push_back(model::term_description(multname_n, varname_u2, false)); // LU2
1875  tl.push_back(model::term_description(multname_n, multname_n, true)); // LL
1876  break;
1877  case 2 :
1878  tl.push_back(model::term_description(varname_u1, multname_n, true)); // U1L
1879  tl.push_back(model::term_description(varname_u2, multname_n, true)); // U2L
1880  tl.push_back(model::term_description(multname_n, multname_n, true)); // LL
1881  tl.push_back(model::term_description(varname_u1, varname_u1, true)); // U1U1
1882  tl.push_back(model::term_description(varname_u2, varname_u2, true)); // U2U2
1883  tl.push_back(model::term_description(varname_u1, varname_u2, true)); // U1U2
1884  break;
1885  default : GMM_ASSERT1(false,
1886  "Incorrect option for integral contact brick");
1887  }
1888 
1889  model::varnamelist dl(1, dataname_r);
1890 
1891  model::varnamelist vl(1, varname_u1);
1892  vl.push_back(varname_u2);
1893  vl.push_back(multname_n);
1894 
1895  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region1);
1896  }
1897 
1898 
1899  //=========================================================================
1900  // Add a contact condition with Coulomb friction between two bodies
1901  // discretized with nonmatching meshes.
1902  //=========================================================================
1903 
1904  size_type add_integral_contact_between_nonmatching_meshes_brick
1905  (model &md, const mesh_im &mim, const std::string &varname_u1,
1906  const std::string &varname_u2, const std::string &multname,
1907  const std::string &dataname_r, const std::string &dataname_friction_coeffs,
1908  size_type region1, size_type region2, int option,
1909  const std::string &dataname_alpha,
1910  const std::string &dataname_wt1, const std::string &dataname_wt2) {
1911 
1912  pbrick pbr = std::make_shared<integral_contact_nonmatching_meshes_brick>
1913  (region1, region2, false /* contact_only */, option);
1914 
1915  model::termlist tl;
1916 
1917  switch (option) {
1918  case 1 : case 3 : case 4 :
1919  tl.push_back(model::term_description(varname_u1, multname, false)); // 0: U1L
1920  tl.push_back(model::term_description(multname, varname_u1, false)); // 1: LU1
1921  tl.push_back(model::term_description(varname_u2, multname, false)); // 2: U2L
1922  tl.push_back(model::term_description(multname, varname_u2, false)); // 3: LU2
1923  tl.push_back(model::term_description(multname, multname, true)); // 4: LL
1924  break;
1925  case 2 :
1926  tl.push_back(model::term_description(varname_u1, multname, false)); // 0: U1L
1927  tl.push_back(model::term_description(multname, varname_u1, false)); // 1: LU1
1928  tl.push_back(model::term_description(varname_u2, multname, false)); // 2: U2L
1929  tl.push_back(model::term_description(multname, varname_u2, false)); // 3: LU2
1930  tl.push_back(model::term_description(multname, multname, true)); // 4: LL
1931  tl.push_back(model::term_description(varname_u1, varname_u1, true)); // 5: U1U1
1932  tl.push_back(model::term_description(varname_u2, varname_u2, true)); // 6: U2U2
1933  tl.push_back(model::term_description(varname_u1, varname_u2, true)); // 7: U1U2
1934  tl.push_back(model::term_description(varname_u2, varname_u1, true)); // 8: U2U1
1935  break;
1936  default : GMM_ASSERT1(false,
1937  "Incorrect option for integral contact brick");
1938  }
1939 
1940  model::varnamelist dl(1, dataname_r); // 0 -> r
1941  dl.push_back(dataname_friction_coeffs); // 1 -> f_coeff,[tau_adh,tresca_lim]
1942  if (dataname_alpha.size()) {
1943  dl.push_back(dataname_alpha); // 2 -> alpha
1944  if (dataname_wt1.size()) {
1945  dl.push_back(dataname_wt1); // 3 -> WT1
1946  if (dataname_wt2.size()) {
1947  dl.push_back(dataname_wt2); // 4 -> WT2
1948  // TODO: VT1, VT2, gamma
1949  }
1950  }
1951  }
1952 
1953  model::varnamelist vl(1, varname_u1);
1954  vl.push_back(varname_u2);
1955  vl.push_back(multname);
1956 
1957  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region1);
1958  }
1959 
1960 
1961  //=========================================================================
1962  //
1963  // Integral penalized contact (with friction) between non-matching meshes.
1964  //
1965  //=========================================================================
1966 
1967  template<typename MAT, typename VECT1>
1968  void asm_penalized_contact_nonmatching_meshes_tangent_matrix // frictionless
1969  (MAT &Ku1u1, MAT &Ku2u2, MAT &Ku1u2,
1970  const mesh_im &mim,
1971  const getfem::mesh_fem &mf_u1, const VECT1 &U1,
1972  const getfem::mesh_fem &mf_u2, const VECT1 &U2,
1973  const getfem::mesh_fem *pmf_lambda, const VECT1 *lambda,
1974  const scalar_type r, const mesh_region &rg, int option = 1) {
1975 
1976  contact_nonmatching_meshes_nonlinear_term
1977  nterm((option == 1) ? K_UU_V1 : K_UU_V2, r,
1978  mf_u1, U1, mf_u2, U2, pmf_lambda, lambda);
1979 
1980  const std::string aux_fems = pmf_lambda ? "#1,#2,#3" : "#1,#2";
1981 
1982  getfem::generic_assembly assem;
1983  assem.set("M$1(#1,#1)+=comp(NonLin(#1," + aux_fems + ").vBase(#1).vBase(#1))(i,j,:,i,:,j); "
1984  "M$2(#2,#2)+=comp(NonLin(#1," + aux_fems + ").vBase(#2).vBase(#2))(i,j,:,i,:,j); "
1985  "M$3(#1,#2)+=comp(NonLin(#1," + aux_fems + ").vBase(#1).vBase(#2))(i,j,:,i,:,j)");
1986  assem.push_mi(mim);
1987  assem.push_mf(mf_u1);
1988  assem.push_mf(mf_u2);
1989  if (pmf_lambda)
1990  assem.push_mf(*pmf_lambda);
1991  assem.push_nonlinear_term(&nterm);
1992  assem.push_mat(Ku1u1);
1993  assem.push_mat(Ku2u2);
1994  assem.push_mat(Ku1u2);
1995  assem.assembly(rg);
1996 
1997  gmm::scale(Ku1u2, scalar_type(-1));
1998  }
1999 
2000  template<typename VECT1>
2001  void asm_penalized_contact_nonmatching_meshes_rhs // frictionless
2002  (VECT1 &Ru1, VECT1 &Ru2,
2003  const mesh_im &mim,
2004  const getfem::mesh_fem &mf_u1, const VECT1 &U1,
2005  const getfem::mesh_fem &mf_u2, const VECT1 &U2,
2006  const getfem::mesh_fem *pmf_lambda, const VECT1 *lambda,
2007  scalar_type r, const mesh_region &rg, int option = 1) {
2008 
2009  contact_nonmatching_meshes_nonlinear_term
2010  nterm((option == 1) ? RHS_U_V5 : RHS_U_V2, r,
2011  mf_u1, U1, mf_u2, U2, pmf_lambda, lambda);
2012 
2013  const std::string aux_fems = pmf_lambda ? "#1,#2,#3": "#1,#2";
2014  getfem::generic_assembly assem;
2015  assem.set("V$1(#1)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1))(i,:,i); "
2016  "V$2(#2)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#2))(i,:,i)");
2017  assem.push_mi(mim);
2018  assem.push_mf(mf_u1);
2019  assem.push_mf(mf_u2);
2020  if (pmf_lambda)
2021  assem.push_mf(*pmf_lambda);
2022  assem.push_nonlinear_term(&nterm);
2023  assem.push_vec(Ru1);
2024  assem.push_vec(Ru2);
2025  assem.assembly(rg);
2026 
2027  gmm::scale(Ru2, scalar_type(-1));
2028  }
2029 
2030 
2031  template<typename MAT, typename VECT1>
2032  void asm_penalized_contact_nonmatching_meshes_tangent_matrix // with friction
2033  (MAT &Ku1u1, MAT &Ku2u2, MAT &Ku1u2, MAT &Ku2u1,
2034  const mesh_im &mim,
2035  const getfem::mesh_fem &mf_u1, const VECT1 &U1,
2036  const getfem::mesh_fem &mf_u2, const VECT1 &U2,
2037  const getfem::mesh_fem *pmf_lambda, const VECT1 *lambda,
2038  const getfem::mesh_fem *pmf_coeff, const VECT1 *f_coeffs, scalar_type r,
2039  scalar_type alpha, const VECT1 *WT1, const VECT1 *WT2,
2040  const mesh_region &rg, int option = 1) {
2041 
2042  size_type subterm = 0;
2043  switch (option) {
2044  case 1 : subterm = K_UU_FRICT_V4; break;
2045  case 2 : subterm = K_UU_FRICT_V3; break;
2046  case 3 : subterm = K_UU_FRICT_V5; break;
2047  }
2048 
2049  contact_nonmatching_meshes_nonlinear_term
2050  nterm(subterm, r, mf_u1, U1, mf_u2, U2, pmf_lambda, lambda,
2051  pmf_coeff, f_coeffs, alpha, WT1, WT2);
2052 
2053  const std::string aux_fems = pmf_coeff ? "#1,#2,#3,#4"
2054  : (pmf_lambda ? "#1,#2,#3": "#1,#2");
2055 
2056  getfem::generic_assembly assem;
2057  assem.set("M$1(#1,#1)+=comp(NonLin(#1," + aux_fems + ").vBase(#1).vBase(#1))(i,j,:,i,:,j); "
2058  "M$2(#2,#2)+=comp(NonLin(#1," + aux_fems + ").vBase(#2).vBase(#2))(i,j,:,i,:,j); "
2059  "M$3(#1,#2)+=comp(NonLin(#1," + aux_fems + ").vBase(#1).vBase(#2))(i,j,:,i,:,j); "
2060  "M$4(#2,#1)+=comp(NonLin(#1," + aux_fems + ").vBase(#2).vBase(#1))(i,j,:,i,:,j)");
2061  assem.push_mi(mim);
2062  assem.push_mf(mf_u1);
2063  assem.push_mf(mf_u2);
2064 
2065  if (pmf_lambda)
2066  assem.push_mf(*pmf_lambda);
2067  else if (pmf_coeff)
2068  assem.push_mf(*pmf_coeff); // dummy
2069 
2070  if (pmf_coeff)
2071  assem.push_mf(*pmf_coeff);
2072 
2073  assem.push_nonlinear_term(&nterm);
2074  assem.push_mat(Ku1u1);
2075  assem.push_mat(Ku2u2);
2076  assem.push_mat(Ku1u2);
2077  assem.push_mat(Ku2u1);
2078  assem.assembly(rg);
2079 
2080  gmm::scale(Ku1u2, scalar_type(-1));
2081  gmm::scale(Ku2u1, scalar_type(-1));
2082  }
2083 
2084 
2085  template<typename VECT1>
2086  void asm_penalized_contact_nonmatching_meshes_rhs // with friction
2087  (VECT1 &Ru1, VECT1 &Ru2,
2088  const mesh_im &mim,
2089  const getfem::mesh_fem &mf_u1, const VECT1 &U1,
2090  const getfem::mesh_fem &mf_u2, const VECT1 &U2,
2091  const getfem::mesh_fem *pmf_lambda, const VECT1 *lambda,
2092  const getfem::mesh_fem *pmf_coeff, const VECT1 *f_coeffs, scalar_type r,
2093  scalar_type alpha, const VECT1 *WT1, const VECT1 *WT2,
2094  const mesh_region &rg, int option = 1) {
2095 
2096  size_type subterm = 0;
2097  switch (option) {
2098  case 1 : subterm = RHS_U_FRICT_V7; break;
2099  case 2 : subterm = RHS_U_FRICT_V6; break;
2100  case 3 : subterm = RHS_U_FRICT_V8; break;
2101  }
2102 
2103  contact_nonmatching_meshes_nonlinear_term
2104  nterm(subterm, r, mf_u1, U1, mf_u2, U2, pmf_lambda, lambda,
2105  pmf_coeff, f_coeffs, alpha, WT1, WT2);
2106 
2107  const std::string aux_fems = pmf_coeff ? "#1,#2,#3,#4"
2108  : (pmf_lambda ? "#1,#2,#3": "#1,#2");
2109  getfem::generic_assembly assem;
2110  assem.set("V$1(#1)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#1))(i,:,i); "
2111  "V$2(#2)+=comp(NonLin$1(#1," + aux_fems + ").vBase(#2))(i,:,i)");
2112 
2113  assem.push_mi(mim);
2114  assem.push_mf(mf_u1);
2115  assem.push_mf(mf_u2);
2116 
2117  if (pmf_lambda)
2118  assem.push_mf(*pmf_lambda);
2119  else if (pmf_coeff)
2120  assem.push_mf(*pmf_coeff); // dummy
2121 
2122  if (pmf_coeff)
2123  assem.push_mf(*pmf_coeff);
2124 
2125  assem.push_nonlinear_term(&nterm);
2126  assem.push_vec(Ru1);
2127  assem.push_vec(Ru2);
2128  assem.assembly(rg);
2129 
2130  gmm::scale(Ru2, scalar_type(-1));
2131  }
2132 
2133 
2134  struct penalized_contact_nonmatching_meshes_brick : public virtual_brick {
2135 
2136  size_type rg1, rg2; // ids of mesh regions on mf_u1 and mf_u2 that are
2137  // expected to come in contact.
2138  mutable getfem::pfem pfem_proj; // cached fem for the projection between nonmatching meshes
2139  bool contact_only;
2140  int option;
2141 
2142  virtual void asm_real_tangent_terms(const model &md, size_type /* ib */,
2143  const model::varnamelist &vl,
2144  const model::varnamelist &dl,
2145  const model::mimlist &mims,
2146  model::real_matlist &matl,
2147  model::real_veclist &vecl,
2148  model::real_veclist &,
2149  size_type region,
2150  build_version version) const {
2151  // Integration method
2152  GMM_ASSERT1(mims.size() == 1,
2153  "Penalized contact between nonmatching meshes bricks need a single mesh_im");
2154  const mesh_im &mim = *mims[0];
2155 
2156  // Variables : u1, u2
2157  GMM_ASSERT1(vl.size() == 2,
2158  "Penalized contact between nonmatching meshes bricks need two variables");
2159  const model_real_plain_vector &u1 = md.real_variable(vl[0]);
2160  const model_real_plain_vector &u2 = md.real_variable(vl[1]);
2161  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(vl[0]);
2162  const mesh_fem &mf_u2 = md.mesh_fem_of_variable(vl[1]);
2163 
2164  size_type N = mf_u1.linked_mesh().dim();
2165 
2166  // Data : r, [lambda,] [friction_coeffs,] [alpha,] [WT1, WT2]
2167  size_type nb_data_1 = ((option == 1) ? 1 : 2) + (contact_only ? 0 : 1);
2168  size_type nb_data_2 = nb_data_1 + (contact_only ? 0 : 3);
2169  GMM_ASSERT1(dl.size() >= nb_data_1 && dl.size() <= nb_data_2,
2170  "Wrong number of data for penalized contact between nonmatching meshes "
2171  << "brick, " << dl.size() << " should be between "
2172  << nb_data_1 << " and " << nb_data_2 << ".");
2173 
2174  size_type nd = 0;
2175  const model_real_plain_vector &vr = md.real_variable(dl[nd]);
2176  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Parameter r should be a scalar");
2177 
2178  size_type sl;
2179  const model_real_plain_vector *lambda = 0;
2180  const mesh_fem *pmf_lambda = 0;
2181  if (option != 1) {
2182  nd++;
2183  lambda = &(md.real_variable(dl[nd]));
2184  pmf_lambda = md.pmesh_fem_of_variable(dl[nd]);
2185  sl = gmm::vect_size(*lambda) * pmf_lambda->get_qdim() / pmf_lambda->nb_dof();
2186  GMM_ASSERT1(sl == (contact_only ? 1 : N),
2187  "the data corresponding to the contact stress "
2188  "has not the right format");
2189  }
2190 
2191  const model_real_plain_vector *f_coeffs = 0;
2192  const mesh_fem *pmf_coeff = 0;
2193  scalar_type alpha = 1;
2194  const model_real_plain_vector *WT1 = 0;
2195  const model_real_plain_vector *WT2 = 0;
2196  if (!contact_only) {
2197  nd++;
2198  f_coeffs = &(md.real_variable(dl[nd]));
2199  pmf_coeff = md.pmesh_fem_of_variable(dl[nd]);
2200  sl = gmm::vect_size(*f_coeffs);
2201  if (pmf_coeff) { sl *= pmf_coeff->get_qdim(); sl /= pmf_coeff->nb_dof(); }
2202  GMM_ASSERT1(sl == 1 || sl == 2 || sl == 3,
2203  "the data corresponding to the friction coefficient "
2204  "has not the right format");
2205 
2206  if (dl.size() > nd+1) {
2207  nd++;
2208  alpha = md.real_variable(dl[nd])[0];
2209  GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[nd])) == 1,
2210  "Parameter alpha should be a scalar");
2211  }
2212 
2213  if (dl.size() > nd+1) {
2214  nd++;
2215  if (dl[nd].compare(vl[0]) != 0)
2216  WT1 = &(md.real_variable(dl[nd]));
2217  else if (md.n_iter_of_variable(vl[0]) > 1)
2218  WT1 = &(md.real_variable(vl[0],1));
2219  }
2220 
2221  if (dl.size() > nd+1) {
2222  nd++;
2223  if (dl[nd].compare(vl[1]) != 0)
2224  WT2 = &(md.real_variable(dl[nd]));
2225  else if (md.n_iter_of_variable(vl[1]) > 1)
2226  WT2 = &(md.real_variable(vl[1],1));
2227  }
2228  }
2229 
2230  GMM_ASSERT1(matl.size() == (contact_only ? 3 : 4),
2231  "Wrong number of terms for penalized contact "
2232  "between nonmatching meshes brick");
2233 
2234  mesh_region rg(region);
2235  mf_u1.linked_mesh().intersect_with_mpi_region(rg); // FIXME: mfu_2?
2236 
2237  // projection of the second mesh_fem onto the mesh of the first mesh_fem
2238  if (!pfem_proj)
2239  pfem_proj = new_projected_fem(mf_u2, mim, rg2, rg1);
2240 
2241  getfem::mesh_fem mf_u2_proj(mim.linked_mesh(), dim_type(N));
2242  mf_u2_proj.set_finite_element(mim.linked_mesh().convex_index(), pfem_proj);
2243 
2244  size_type nbdof1 = mf_u1.nb_dof();
2245  size_type nbdof2 = mf_u2.nb_dof();
2246  size_type nbsub = mf_u2_proj.nb_dof();
2247 
2248  std::vector<size_type> ind;
2249  mf_u2_proj.get_global_dof_index(ind);
2250  gmm::unsorted_sub_index SUBI(ind);
2251 
2252  gmm::csc_matrix<scalar_type> Esub(nbsub, nbdof2);
2253  if (mf_u2.is_reduced())
2254  gmm::copy(gmm::sub_matrix(mf_u2.extension_matrix(),
2255  SUBI, gmm::sub_interval(0, nbdof2)),
2256  Esub);
2257 
2258  model_real_plain_vector u2_proj(nbsub);
2259  if (mf_u2.is_reduced())
2260  gmm::mult(Esub, u2, u2_proj);
2261  else
2262  gmm::copy(gmm::sub_vector(u2, SUBI), u2_proj);
2263 
2264  model_real_plain_vector WT2_proj(0);
2265  if (WT2) {
2266  gmm::resize(WT2_proj, nbsub);
2267  if (mf_u2.is_reduced())
2268  gmm::mult(Esub, *WT2, WT2_proj);
2269  else
2270  gmm::copy(gmm::sub_vector(*WT2, SUBI), WT2_proj);
2271  }
2272 
2273  if (version & model::BUILD_MATRIX) {
2274  GMM_TRACE2("Penalized contact between nonmatching meshes tangent term");
2275  gmm::clear(matl[0]);
2276  gmm::clear(matl[1]);
2277  gmm::clear(matl[2]);
2278 
2279  model_real_sparse_matrix Ku2u2(nbsub,nbsub);
2280  model_real_sparse_matrix Ku1u2(nbdof1,nbsub);
2281 
2282  if (contact_only) {
2283  asm_penalized_contact_nonmatching_meshes_tangent_matrix
2284  (matl[0], Ku2u2, Ku1u2, mim, mf_u1, u1, mf_u2_proj, u2_proj,
2285  pmf_lambda, lambda, vr[0], rg, option);
2286  }
2287  else {
2288  gmm::clear(matl[3]);
2289  model_real_sparse_matrix Ku2u1(nbsub,nbdof1);
2290  asm_penalized_contact_nonmatching_meshes_tangent_matrix
2291  (matl[0], Ku2u2, Ku1u2, Ku2u1, mim, mf_u1, u1, mf_u2_proj, u2_proj,
2292  pmf_lambda, lambda, pmf_coeff, f_coeffs, vr[0], alpha, WT1, &WT2_proj, rg, option);
2293  if (mf_u2.is_reduced())
2294  gmm::mult(gmm::transposed(Esub), Ku2u1, matl[3]);
2295  else
2296  gmm::copy(Ku2u1, gmm::sub_matrix(matl[3], SUBI, gmm::sub_interval(0, nbdof1)));
2297  }
2298 
2299  if (mf_u2.is_reduced()) {
2300  model_real_sparse_matrix tmp(nbsub, nbdof2);
2301  gmm::mult(Ku2u2, Esub, tmp);
2302  gmm::mult(gmm::transposed(Esub), tmp, matl[1]);
2303  gmm::mult(Ku1u2, Esub, matl[2]);
2304  }
2305  else {
2306  gmm::copy(Ku2u2, gmm::sub_matrix(matl[1], SUBI));
2307  gmm::copy(Ku1u2, gmm::sub_matrix(matl[2], gmm::sub_interval(0, nbdof1), SUBI));
2308  }
2309  }
2310 
2311  if (version & model::BUILD_RHS) {
2312  gmm::clear(vecl[0]);
2313  gmm::clear(vecl[1]);
2314 
2315  model_real_plain_vector Ru2(nbsub);
2316 
2317  if (contact_only)
2318  asm_penalized_contact_nonmatching_meshes_rhs
2319  (vecl[0], Ru2, mim, mf_u1, u1, mf_u2_proj, u2_proj, pmf_lambda, lambda,
2320  vr[0], rg, option);
2321  else
2322  asm_penalized_contact_nonmatching_meshes_rhs
2323  (vecl[0], Ru2, mim, mf_u1, u1, mf_u2_proj, u2_proj, pmf_lambda, lambda,
2324  pmf_coeff, f_coeffs, vr[0], alpha, WT1, &WT2_proj, rg, option);
2325 
2326  if (mf_u2.is_reduced())
2327  gmm::mult(gmm::transposed(Esub), Ru2, vecl[1]);
2328  else
2329  gmm::copy(Ru2, gmm::sub_vector(vecl[1], SUBI));
2330  }
2331  }
2332 
2333  penalized_contact_nonmatching_meshes_brick(size_type rg1_, size_type rg2_,
2334  bool contact_only_, int option_)
2335  : rg1(rg1_), rg2(rg2_), pfem_proj(0),
2336  contact_only(contact_only_), option(option_)
2337  {
2338  set_flags(contact_only
2339  ? "Integral penalized contact between nonmatching meshes brick"
2340  : "Integral penalized contact and friction between nonmatching "
2341  "meshes brick",
2342  false /* is linear*/, contact_only /* is symmetric */,
2343  true /* is coercive */, true /* is real */,
2344  false /* is complex */);
2345  }
2346 
2347  ~penalized_contact_nonmatching_meshes_brick()
2348  { if (pfem_proj) del_projected_fem(pfem_proj); }
2349 
2350  };
2351 
2352 
2353  //=========================================================================
2354  // Add a frictionless contact condition between two bodies discretized
2355  // with nonmatching meshes.
2356  //=========================================================================
2357 
2358  size_type add_penalized_contact_between_nonmatching_meshes_brick
2359  (model &md, const mesh_im &mim, const std::string &varname_u1,
2360  const std::string &varname_u2, const std::string &dataname_r,
2361  size_type region1, size_type region2,
2362  int option, const std::string &dataname_n) {
2363 
2364  pbrick pbr = std::make_shared<penalized_contact_nonmatching_meshes_brick>
2365  (region1, region2, true /* contact_only */, option);
2366  model::termlist tl;
2367  tl.push_back(model::term_description(varname_u1, varname_u1, true));
2368  tl.push_back(model::term_description(varname_u2, varname_u2, true));
2369  tl.push_back(model::term_description(varname_u1, varname_u2, true));
2370 
2371  model::varnamelist dl(1, dataname_r);
2372  switch (option) {
2373  case 1: break;
2374  case 2: dl.push_back(dataname_n); break;
2375  default: GMM_ASSERT1(false, "Penalized contact brick : invalid option");
2376  }
2377 
2378  model::varnamelist vl(1, varname_u1);
2379  vl.push_back(varname_u2);
2380 
2381  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region1);
2382  }
2383 
2384 
2385  //=========================================================================
2386  // Add a contact condition with friction between two bodies discretized
2387  // with nonmatching meshes.
2388  //=========================================================================
2389 
2390  size_type add_penalized_contact_between_nonmatching_meshes_brick
2391  (model &md, const mesh_im &mim, const std::string &varname_u1,
2392  const std::string &varname_u2, const std::string &dataname_r,
2393  const std::string &dataname_friction_coeffs,
2394  size_type region1, size_type region2, int option,
2395  const std::string &dataname_lambda, const std::string &dataname_alpha,
2396  const std::string &dataname_wt1, const std::string &dataname_wt2) {
2397 
2398  pbrick pbr = std::make_shared<penalized_contact_nonmatching_meshes_brick>
2399  (region1, region2, false /* contact_only */, option);
2400  model::termlist tl;
2401  tl.push_back(model::term_description(varname_u1, varname_u1, true)); // 0: U1U1
2402  tl.push_back(model::term_description(varname_u2, varname_u2, true)); // 1: U2U2
2403  tl.push_back(model::term_description(varname_u1, varname_u2, true)); // 2: U1U2
2404  tl.push_back(model::term_description(varname_u2, varname_u1, true)); // 3: U2U1
2405 
2406  model::varnamelist dl(1, dataname_r);
2407  switch (option) {
2408  case 1: break;
2409  case 2: case 3: dl.push_back(dataname_lambda); break;
2410  default: GMM_ASSERT1(false, "Penalized contact brick : invalid option");
2411  }
2412  dl.push_back(dataname_friction_coeffs);
2413  if (dataname_alpha.size() > 0) {
2414  dl.push_back(dataname_alpha);
2415  if (dataname_wt1.size() > 0) {
2416  dl.push_back(dataname_wt1);
2417  if (dataname_wt2.size() > 0)
2418  dl.push_back(dataname_wt2);
2419  }
2420  }
2421 
2422  model::varnamelist vl(1, varname_u1);
2423  vl.push_back(varname_u2);
2424 
2425  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region1);
2426  }
2427 
2428 
2429  // Computation of contact area and contact forces
2430  void compute_integral_contact_area_and_force
2431  (model &md, size_type indbrick, scalar_type &area,
2432  model_real_plain_vector &F) {
2433 
2434  pbrick pbr = md.brick_pointer(indbrick);
2435  const model::mimlist &ml = md.mimlist_of_brick(indbrick);
2436  const model::varnamelist &vl = md.varnamelist_of_brick(indbrick);
2437  const model::varnamelist &dl = md.datanamelist_of_brick(indbrick);
2438  size_type reg = md.region_of_brick(indbrick);
2439 
2440  GMM_ASSERT1(ml.size() == 1, "Wrong size");
2441 
2442  if (pbr->brick_name() == "Integral contact with rigid obstacle brick" ||
2443  pbr->brick_name() == "Integral contact and friction with rigid obstacle brick") {
2444  integral_contact_rigid_obstacle_brick *p
2445  = dynamic_cast<integral_contact_rigid_obstacle_brick *>
2446  (const_cast<virtual_brick *>(pbr.get()));
2447  GMM_ASSERT1(p, "Wrong type of brick");
2448 
2449  GMM_ASSERT1(vl.size() >= 2, "Wrong size");
2450  const model_real_plain_vector &u = md.real_variable(vl[0]);
2451  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
2452  const model_real_plain_vector &lambda = md.real_variable(vl[1]);
2453  const mesh_fem &mf_lambda = md.mesh_fem_of_variable(vl[1]);
2454  GMM_ASSERT1(dl.size() >= 1, "Wrong size");
2455  const model_real_plain_vector &obs = md.real_variable(dl[0]);
2456  const mesh_fem &mf_obs = md.mesh_fem_of_variable(dl[0]);
2457 
2458  //FIXME: use an adapted integration method
2459  area = asm_level_set_contact_area(*ml[0], mf_u, u, mf_obs, obs, reg, -1e-3,
2460  &mf_lambda, &lambda, 1e-1);
2461 
2462  gmm::resize(F, mf_u.nb_dof());
2463  asm_level_set_normal_source_term
2464  (F, *ml[0], mf_u, mf_obs, obs, mf_lambda, lambda, reg);
2465  }
2466  else if (pbr->brick_name() == "Integral penalized contact with rigid obstacle brick" ||
2467  pbr->brick_name() == "Integral penalized contact and friction with rigid "
2468  "obstacle brick") {
2469  penalized_contact_rigid_obstacle_brick *p
2470  = dynamic_cast<penalized_contact_rigid_obstacle_brick *>
2471  (const_cast<virtual_brick *>(pbr.get()));
2472  GMM_ASSERT1(p, "Wrong type of brick");
2473  GMM_ASSERT1(false, "Not implemented yet");
2474  }
2475  else if (pbr->brick_name() == "Integral contact between nonmatching meshes brick" ||
2476  pbr->brick_name() == "Integral contact and friction between nonmatching "
2477  "meshes brick") {
2478  integral_contact_nonmatching_meshes_brick *p
2479  = dynamic_cast<integral_contact_nonmatching_meshes_brick *>
2480  (const_cast<virtual_brick *>(pbr.get()));
2481  GMM_ASSERT1(p, "Wrong type of brick");
2482 
2483  GMM_ASSERT1(vl.size() == 3, "Wrong size");
2484  const model_real_plain_vector &u1 = md.real_variable(vl[0]);
2485  const model_real_plain_vector &u2 = md.real_variable(vl[1]);
2486  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(vl[0]);
2487  const mesh_fem &mf_u2 = md.mesh_fem_of_variable(vl[1]);
2488  const model_real_plain_vector &lambda = md.real_variable(vl[2]);
2489  const mesh_fem &mf_lambda = md.mesh_fem_of_variable(vl[2]);
2490 
2491  getfem::pfem pfem_proj = new_projected_fem(mf_u2, *ml[0], p->rg2, p->rg1);
2492  getfem::mesh_fem mf_u2_proj(mf_u1.linked_mesh(), mf_u1.linked_mesh().dim());
2493  mf_u2_proj.set_finite_element(mf_u1.linked_mesh().convex_index(), pfem_proj);
2494 
2495  std::vector<size_type> ind;
2496  mf_u2_proj.get_global_dof_index(ind);
2497  gmm::unsorted_sub_index SUBI(ind);
2498 
2499  size_type nbdof2 = mf_u2.nb_dof();
2500  size_type nbsub = mf_u2_proj.nb_basic_dof();
2501  model_real_plain_vector u2_proj(nbsub);
2502 
2503  if (mf_u2.is_reduced()) {
2504  gmm::csc_matrix<scalar_type> Esub(nbsub, nbdof2);
2505  gmm::copy(gmm::sub_matrix(mf_u2.extension_matrix(),
2506  SUBI, gmm::sub_interval(0, nbdof2)),
2507  Esub);
2508  gmm::mult(Esub, u2, u2_proj);
2509  }
2510  else
2511  gmm::copy(gmm::sub_vector(u2, SUBI), u2_proj);
2512 
2513  //FIXME: use an adapted integration method
2514  area = asm_nonmatching_meshes_contact_area
2515  (*ml[0], mf_u1, u1, mf_u2_proj, u2_proj, reg, -1e-3,
2516  &mf_lambda, &lambda, 1e-1);
2517 
2518  gmm::resize(F, mf_u1.nb_dof());
2519  asm_nonmatching_meshes_normal_source_term
2520  (F, *ml[0], mf_u1, mf_u2_proj, mf_lambda, lambda, reg);
2521 
2522  del_projected_fem(pfem_proj);
2523  }
2524  else if (pbr->brick_name() == "Integral penalized contact between nonmatching meshes brick" ||
2525  pbr->brick_name() == "Integral penalized contact and friction between nonmatching "
2526  "meshes brick") {
2527  penalized_contact_nonmatching_meshes_brick *p
2528  = dynamic_cast<penalized_contact_nonmatching_meshes_brick *>
2529  (const_cast<virtual_brick *>(pbr.get()));
2530  GMM_ASSERT1(p, "Wrong type of brick");
2531  GMM_ASSERT1(false, "Not implemented yet");
2532  }
2533 
2534  }
2535 
2536  //=========================================================================
2537  //
2538  // Contact condition with a rigid obstacle : generic Nitsche's method
2539  //
2540  //=========================================================================
2541 
2542  size_type add_Nitsche_contact_with_rigid_obstacle_brick
2543  (model &md, const mesh_im &mim, const std::string &varname_u,
2544  const std::string &Neumannterm,
2545  const std::string &obs, const std::string &dataname_gamma0,
2546  scalar_type theta_,
2547  std::string dataname_friction_coeff,
2548  const std::string &dataname_alpha,
2549  const std::string &dataname_wt,
2550  size_type region) {
2551 
2552  std::string theta = std::to_string(theta_);
2553  ga_workspace workspace(md, ga_workspace::inherit::ALL); // reenables vars
2554  size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
2555  GMM_ASSERT1(order == 0, "Wrong expression of the Neumann term");
2556  // model::varnamelist vl, vl_test1, vl_test2, dl;
2557  // bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 1);
2558 
2559  std::string gamma = "(("+dataname_gamma0+")/element_size)";
2560  std::string thetagamma = "("+theta+"/"+gamma+")";
2561  std::string contact_normal = "(-Normalized(Grad_"+obs+"))";
2562  std::string gap = "("+obs+"-"+varname_u+"."+contact_normal+")";
2563  std::string Vs = "("+varname_u +
2564  (dataname_wt.size() ? "-("+dataname_wt+"))" : ")");
2565  if (dataname_alpha.size()) Vs = "(("+dataname_alpha+")*"+Vs;
2566  if (dataname_friction_coeff.size() == 0)
2567  dataname_friction_coeff = std::string("0");
2568 
2569  std::string expr = "Coulomb_friction_coupled_projection("+Neumannterm
2570  +", "+contact_normal+", "+Vs+", "+gap+", "+dataname_friction_coeff
2571  +","+gamma+").(";
2572 
2573  if (theta_ != scalar_type(0)) {
2574  std::string derivative_Neumann=workspace.extract_order1_term(varname_u);
2575  expr = "-"+thetagamma+"*("+Neumannterm+").("+derivative_Neumann
2576  +") + "+expr+thetagamma+"*("+derivative_Neumann+")";
2577  }
2578  expr += "-Test_"+varname_u+")";
2579 
2580  return add_nonlinear_term
2581  (md, mim, expr, region, false, false,
2582  "Nitsche contact with rigid obstacle");
2583  }
2584 
2585 #ifdef EXPERIMENTAL_PURPOSE_ONLY
2586 
2587  //=========================================================================
2588  //
2589  // Contact condition with a rigid obstacle : generic Nitsche's method
2590  // Experimental for midpoint scheme
2591  //
2592  //=========================================================================
2593 
2594 
2595  size_type add_Nitsche_contact_with_rigid_obstacle_brick_modified_midpoint
2596  (model &md, const mesh_im &mim, const std::string &varname_u,
2597  const std::string &Neumannterm,
2598  const std::string &Neumannterm_wt,
2599  const std::string &obs, const std::string &dataname_gamma0,
2600  scalar_type theta_,
2601  std::string dataname_friction_coeff,
2602  const std::string &/* dataname_alpha*/,
2603  const std::string &dataname_wt,
2604  size_type region) {
2605 
2606  std::string theta = std::to_string(theta_);
2607  ga_workspace workspace(md, ga_workspace::inherit::ALL);
2608  size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
2609  GMM_ASSERT1(order == 0, "Wrong expression of the Neumann term");
2610 
2611  std::string gamma = "((1/("+dataname_gamma0+"))*element_size)";
2612  std::string thetagamma = "("+theta+"*"+gamma+")";
2613  std::string contact_normal = "(-Normalized(Grad_"+obs+"))";
2614  std::string gap = "("+obs+"-"+varname_u+"."+contact_normal+")";
2615  std::string gap_wt = "("+obs+"-"+dataname_wt+"."+contact_normal+")";
2616  std::string Vs = "("+varname_u +
2617  (dataname_wt.size() ? "-("+dataname_wt+"))" : ")");
2618  std::string Vs_wt = Vs; // To be adapted ...
2619  if (dataname_friction_coeff.size() == 0)
2620  dataname_friction_coeff = std::string("0");
2621 
2622  std::string Pg_wt = "(-"+gamma+"*("+Neumannterm_wt+")."+contact_normal+
2623  " - "+gap_wt+")";
2624 
2625  std::string expr
2626  ="((Heaviside("+Pg_wt+")*Coulomb_friction_coupled_projection("
2627  +Neumannterm+", "+contact_normal+", "+Vs+", "+gap+", "
2628  +dataname_friction_coeff+",1/"+gamma+"))+"
2629  +"((1-Heaviside("+Pg_wt+"))*(Coulomb_friction_coupled_projection("
2630  +Neumannterm_wt+", "+contact_normal+", "+Vs_wt+", "+gap_wt+", "
2631  +dataname_friction_coeff+",1/"+gamma+")*0.5+"
2632  "Coulomb_friction_coupled_projection(2*("
2633  +Neumannterm+")-("+Neumannterm_wt+"), "+contact_normal
2634  +", 2*"+Vs+"-"+Vs_wt+", 2*"+gap+"-"+gap_wt+", "
2635  +dataname_friction_coeff+",1/"+gamma+")*0.5))).(";
2636 
2637  if (theta_ != scalar_type(0)) {
2638  std::string derivative_Neumann=workspace.extract_order1_term(varname_u);
2639  expr = "-"+thetagamma+"*("+Neumannterm+").("+derivative_Neumann
2640  +") + "+expr+thetagamma+"*("+derivative_Neumann+")";
2641  }
2642  expr += "-Test_"+varname_u+")";
2643 
2644  return add_nonlinear_term
2645  (md, mim, expr, region, false, false,
2646  "Nitsche contact with rigid obstacle");
2647  }
2648 
2649 
2650 #endif
2651 
2652 
2653 } /* end of namespace getfem. */
bgeot_rtree.h
region-tree for window/point search on a set of rectangles.
getfem::generic_assembly
Generic assembly of vectors, matrices.
Definition: getfem_assembling_tensors.h:606
getfem::generic_assembly::push_vec
void push_vec(VEC &v)
Add a new output vector.
Definition: getfem_assembling_tensors.h:662
getfem::generic_assembly::push_mat
void push_mat(const MAT &m)
Add a new output matrix (fake const version..)
Definition: getfem_assembling_tensors.h:670
getfem::generic_assembly::assembly
void assembly(const mesh_region &region=mesh_region::all_convexes())
do the assembly on the specified region (boundary or set of convexes)
Definition: getfem_assembling_tensors.cc:1820
getfem::generic_assembly::push_nonlinear_term
void push_nonlinear_term(pnonlinear_elem_term net)
Add a new non-linear term.
Definition: getfem_assembling_tensors.h:654
getfem::generic_assembly::push_mi
void push_mi(const mesh_im &im_)
Add a new mesh_im.
Definition: getfem_assembling_tensors.h:651
getfem::generic_assembly::push_mf
void push_mf(const mesh_fem &mf_)
Add a new mesh_fem.
Definition: getfem_assembling_tensors.h:649
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::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
getfem::model::add_brick
size_type add_brick(pbrick pbr, const varnamelist &varnames, const varnamelist &datanames, const termlist &terms, const mimlist &mims, size_type region)
Add a brick to the model.
Definition: getfem_models.cc:1031
getfem_contact_and_friction_common.h
Comomon tools for unilateral contact and Coulomb friction bricks.
getfem_contact_and_friction_integral.h
Unilateral contact and Coulomb friction condition brick.
getfem_projected_fem.h
FEM which projects a mesh_fem on a different mesh.
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::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:263
gmm_condition_number.h
computation of the condition number of dense matrices.
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:243
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:46
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::add_integral_contact_between_nonmatching_meshes_brick
size_type add_integral_contact_between_nonmatching_meshes_brick(model &md, const mesh_im &mim, const std::string &varname_u1, const std::string &varname_u2, const std::string &multname_n, const std::string &dataname_r, size_type region1, size_type region2, int option=1)
Add a frictionless contact condition between nonmatching meshes to the model.
Definition: getfem_contact_and_friction_integral.cc:1859
getfem::add_nonlinear_term
size_type APIDECL add_nonlinear_term(model &md, const mesh_im &mim, const std::string &expr, size_type region=size_type(-1), bool is_sym=false, bool is_coercive=false, const std::string &brickname="")
Add a nonlinear term given by the weak form language expression expr which will be assembled in regio...
Definition: getfem_models.cc:3675
getfem::add_integral_contact_with_rigid_obstacle_brick
size_type add_integral_contact_with_rigid_obstacle_brick(model &md, const mesh_im &mim, const std::string &varname_u, const std::string &multname_n, const std::string &dataname_obs, const std::string &dataname_r, size_type region, int option=1)
Add a frictionless contact condition with a rigid obstacle to the model, which is defined in an integ...
Definition: getfem_contact_and_friction_integral.cc:921
getfem::add_penalized_contact_with_rigid_obstacle_brick
size_type add_penalized_contact_with_rigid_obstacle_brick(model &md, const mesh_im &mim, const std::string &varname_u, const std::string &dataname_obs, const std::string &dataname_r, size_type region, int option=1, const std::string &dataname_lambda_n="")
Add a penalized contact frictionless condition with a rigid obstacle to the model.
Definition: getfem_contact_and_friction_integral.cc:1296
getfem::asm_level_set_normal_source_term
void asm_level_set_normal_source_term(VEC &R, const mesh_im &mim, const getfem::mesh_fem &mf_u, const getfem::mesh_fem &mf_obs, const VEC &obs, const getfem::mesh_fem &mf_lambda, const VEC &lambda, const mesh_region &rg)
Specific assembly procedure for the use of an Uzawa algorithm to solve contact problems.
Definition: getfem_contact_and_friction_integral.h:626
getfem::del_projected_fem
void del_projected_fem(pfem pf)
release a projected fem
Definition: getfem_projected_fem.h:178
getfem::add_penalized_contact_between_nonmatching_meshes_brick
size_type add_penalized_contact_between_nonmatching_meshes_brick(model &md, const mesh_im &mim, const std::string &varname_u1, const std::string &varname_u2, const std::string &dataname_r, size_type region1, size_type region2, int option=1, const std::string &dataname_lambda_n="")
Add a penalized contact frictionless condition between nonmatching meshes to the model.
Definition: getfem_contact_and_friction_integral.cc:2359
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:662
getfem::pbrick
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
Definition: getfem_models.h:48
getfem::new_projected_fem
pfem new_projected_fem(const mesh_fem &mf_source, const mesh_im &mim_target, size_type rg_source_=size_type(-1), size_type rg_target_=size_type(-1), dal::bit_vector blocked_dofs=dal::bit_vector(), bool store_val=true)
create a new projected FEM.
Definition: getfem_projected_fem.cc:880
getfem::add_Nitsche_contact_with_rigid_obstacle_brick
size_type add_Nitsche_contact_with_rigid_obstacle_brick(model &md, const mesh_im &mim, const std::string &varname_u, const std::string &Neumannterm, const std::string &expr_obs, const std::string &dataname_gamma0, scalar_type theta_, std::string dataexpr_friction_coeff, const std::string &dataname_alpha, const std::string &dataname_wt, size_type region)
Adds a contact condition with or without Coulomb friction on the variable varname_u and the mesh boun...
Definition: getfem_contact_and_friction_integral.cc:2543