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