GetFEM  5.5
getfem_contact_and_friction_integral.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2011-2026 Yves Renard, Konstantinos Poulios.
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program. If not, see https://www.gnu.org/licenses/.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 
31 /** @file getfem_contact_and_friction_integral.h
32  @author Yves Renard <Yves.Renard@insa-lyon.fr>
33  @author Konstantinos Poulios <logari81@googlemail.com>
34  @date November, 2011.
35  @brief Unilateral contact and Coulomb friction condition brick.
36  */
37 #ifndef GETFEM_CONTACT_AND_FRICTION_INTEGRAL_H__
38 #define GETFEM_CONTACT_AND_FRICTION_INTEGRAL_H__
39 
40 #include "getfem_models.h"
42 
43 namespace getfem {
44 
45 
46  /** Add a frictionless contact condition with a rigid obstacle
47  to the model, which is defined in an integral way. It is the direct
48  approximation of an augmented Lagrangian formulation (see Getfem user
49  documentation) defined at the continuous level. The advantage should be
50  a better scalability: the number of Newton iterations should be more or
51  less independent of the mesh size.
52  The condition is applied on the variable `varname_u`
53  on the boundary corresponding to `region`. The rigid obstacle should
54  be described with the data `dataname_obstacle` being a signed distance to
55  the obstacle (interpolated on a finite element method).
56  `multname_n` should be a fem variable representing the contact stress.
57  An inf-sup condition between `multname_n` and `varname_u` is required.
58  The augmentation parameter `dataname_r` should be chosen in a
59  range of acceptable values.
60  Possible values for `option` is 1 for the non-symmetric Alart-Curnier
61  augmented Lagrangian method, 2 for the symmetric one, 3 for the
62  non-symmetric Alart-Curnier method with an additional augmentation.
63  The default value is 1.
64  */
66  (model &md, const mesh_im &mim, const std::string &varname_u,
67  const std::string &multname_n, const std::string &dataname_obs,
68  const std::string &dataname_r, size_type region, int option = 1);
69 
70  /** Add a contact with friction condition with a rigid obstacle
71  to the model, which is defined in an integral way. It is the direct
72  approximation of an augmented Lagrangian formulation (see Getfem user
73  documentation) defined at the continuous level. The advantage should be
74  a better scalability: the number of Newton iterations should be more or
75  less independent of the mesh size.
76  The condition is applied on the variable `varname_u`
77  on the boundary corresponding to `region`. The rigid obstacle should
78  be described with the data `dataname_obstacle` being a signed distance
79  to the obstacle (interpolated on a finite element method).
80  `multname` should be a fem variable representing the contact and
81  friction stress.
82  An inf-sup condition between `multname` and `varname_u` is required.
83  The augmentation parameter `dataname_r` should be chosen in a
84  range of acceptable values.
85  The parameter `dataname_friction_coeffs` contains the Coulomb friction
86  coefficient and optionally an adhesional shear stress threshold and the
87  tresca limit shear stress. For constant coefficients its size is from
88  1 to 3. For coefficients described on a finite element method, this
89  vector contains a number of single values, value pairs or triplets
90  equal to the number of the corresponding mesh_fem's basic dofs.
91  Possible values for `option` is 1 for the non-symmetric Alart-Curnier
92  augmented Lagrangian method, 2 for the symmetric one, 3 for the
93  non-symmetric Alart-Curnier method with an additional augmentation
94  and 4 for a new unsymmetric method. The default value is 1.
95  (Option 4 ignores any adhesional stress and tresca limit coefficients.)
96  `dataname_alpha` and `dataname_wt` are optional parameters to solve
97  evolutionary friction problems. `dataname_gamma` and `dataname_vt`
98  represent optional data for adding a parameter-dependent sliding
99  velocity to the friction condition.
100  */
102  (model &md, const mesh_im &mim, const std::string &varname_u,
103  const std::string &multname, const std::string &dataname_obs,
104  const std::string &dataname_r, const std::string &dataname_friction_coeffs,
105  size_type region, int option = 1, const std::string &dataname_alpha = "",
106  const std::string &dataname_wt = "",
107  const std::string &dataname_gamma = "",
108  const std::string &dataname_vt = "");
109 
110 
111  /** Add a penalized contact frictionless condition with a rigid obstacle
112  to the model.
113  The condition is applied on the variable `varname_u`
114  on the boundary corresponding to `region`. The rigid obstacle should
115  be described with the data `dataname_obstacle` being a signed distance to
116  the obstacle (interpolated on a finite element method).
117  The penalization parameter `dataname_r` should be chosen
118  large enough to prescribe an approximate non-penetration condition
119  but not too large not to deteriorate too much the conditionning of
120  the tangent system. `dataname_n` is an optional parameter used if option
121  is 2. In that case, the penalization term is shifted by lambda_n (this
122  allows the use of an Uzawa algorithm on the corresponding augmented
123  Lagrangian formulation)
124  */
126  (model &md, const mesh_im &mim, const std::string &varname_u,
127  const std::string &dataname_obs, const std::string &dataname_r,
128  size_type region, int option = 1, const std::string &dataname_lambda_n = "");
129 
130  /** Add a penalized contact condition with Coulomb friction with a
131  rigid obstacle to the model.
132  The condition is applied on the variable `varname_u`
133  on the boundary corresponding to `region`. The rigid obstacle should
134  be described with the data `dataname_obstacle` being a signed distance to
135  the obstacle (interpolated on a finite element method).
136  The parameter `dataname_friction_coeffs` contains the Coulomb friction
137  coefficient and optionally an adhesional shear stress threshold and the
138  tresca limit shear stress. For constant coefficients its size is from
139  1 to 3. For coefficients described on a finite element method, this
140  vector contains a number of single values, value pairs or triplets
141  equal to the number of the corresponding mesh_fem's basic dofs.
142  The penalization parameter `dataname_r` should be chosen
143  large enough to prescribe approximate non-penetration and friction
144  conditions but not too large not to deteriorate too much the
145  conditionning of the tangent system.
146  `dataname_lambda` is an optional parameter used if option
147  is 2. In that case, the penalization term is shifted by lambda (this
148  allows the use of an Uzawa algorithm on the corresponding augmented
149  Lagrangian formulation)
150  `dataname_alpha` and `dataname_wt` are optional parameters to solve
151  evolutionary friction problems.
152  */
154  (model &md, const mesh_im &mim, const std::string &varname_u,
155  const std::string &dataname_obs, const std::string &dataname_r,
156  const std::string &dataname_friction_coeffs,
157  size_type region, int option = 1, const std::string &dataname_lambda = "",
158  const std::string &dataname_alpha = "",
159  const std::string &dataname_wt = "");
160 
161 
162  /** Add a frictionless contact condition between nonmatching meshes
163  to the model. This brick adds a contact which is defined
164  in an integral way. It is the direct approximation of an augmented
165  Lagrangian formulation (see Getfem user documentation) defined at the
166  continuous level. The advantage should be a better scalability:
167  the number of Newton iterations should be more or less independent
168  of the mesh size.
169  The condition is applied on the variables `varname_u1` and `varname_u2`
170  on the boundaries corresponding to `region1` and `region2`.
171  `multname_n` should be a fem variable representing the contact stress.
172  An inf-sup condition between `multname_n` and `varname_u1` and
173  `varname_u2` is required.
174  The augmentation parameter `dataname_r` should be chosen in a
175  range of acceptable values.
176  Possible values for `option` is 1 for the non-symmetric Alart-Curnier
177  augmented Lagrangian method, 2 for the symmetric one, 3 for the
178  non-symmetric Alart-Curnier method with an additional augmentation.
179  The default value is 1.
180  */
182  (model &md, const mesh_im &mim, const std::string &varname_u1,
183  const std::string &varname_u2, const std::string &multname_n,
184  const std::string &dataname_r,
185  size_type region1, size_type region2, int option = 1);
186 
187  /** Add a contact with friction condition between nonmatching meshes
188  to the model. This brick adds a contact which is defined
189  in an integral way. It is the direct approximation of an augmented
190  Lagrangian formulation (see Getfem user documentation) defined at the
191  continuous level. The advantage should be a better scalability:
192  the number of Newton iterations should be more or less independent
193  of the mesh size.
194  The condition is applied on the variables `varname_u1` and `varname_u2`
195  on the boundaries corresponding to `region1` and `region2`.
196  `multname` should be a fem variable representing the contact and
197  friction stress.
198  An inf-sup condition between `multname` and `varname_u1` and
199  `varname_u2` is required.
200  The augmentation parameter `dataname_r` should be chosen in a
201  range of acceptable values.
202  The parameter `dataname_friction_coeffs` contains the Coulomb friction
203  coefficient and optionally an adhesional shear stress threshold and the
204  tresca limit shear stress. For constant coefficients its size is from
205  1 to 3. For coefficients described on a finite element method on the
206  same mesh as `varname_u1`, this vector contains a number of single
207  values, value pairs or triplets equal to the number of the
208  corresponding mesh_fem's basic dofs.
209  Possible values for `option` is 1 for the non-symmetric Alart-Curnier
210  augmented Lagrangian method, 2 for the symmetric one, 3 for the
211  non-symmetric Alart-Curnier method with an additional augmentation
212  and 4 for a new unsymmetric method. The default value is 1.
213  (Option 4 ignores any adhesional stress and tresca limit coefficients.)
214  `dataname_alpha`, `dataname_wt1` and `dataname_wt2` are optional
215  parameters to solve evolutionary friction problems.
216  */
218  (model &md, const mesh_im &mim, const std::string &varname_u1,
219  const std::string &varname_u2, const std::string &multname,
220  const std::string &dataname_r, const std::string &dataname_friction_coeffs,
221  size_type region1, size_type region2, int option = 1,
222  const std::string &dataname_alpha = "",
223  const std::string &dataname_wt1 = "",
224  const std::string &dataname_wt2 = "");
225 
226 
227  /** Add a penalized contact frictionless condition between nonmatching
228  meshes to the model.
229  The condition is applied on the variables `varname_u1` and `varname_u2`
230  on the boundaries corresponding to `region1` and `region2`.
231  The penalization parameter `dataname_r` should be chosen
232  large enough to prescribe an approximate non-penetration condition
233  but not too large not to deteriorate too much the conditionning of
234  the tangent system. `dataname_n` is an optional parameter used if option
235  is 2. In that case, the penalization term is shifted by lambda_n (this
236  allows the use of an Uzawa algorithm on the corresponding augmented
237  Lagrangian formulation)
238  */
240  (model &md, const mesh_im &mim, const std::string &varname_u1,
241  const std::string &varname_u2, const std::string &dataname_r,
242  size_type region1, size_type region2,
243  int option = 1, const std::string &dataname_lambda_n = "");
244 
245 
246  /** Add a penalized contact condition with Coulomb friction between
247  nonmatching meshes to the model.
248  The condition is applied on the variables `varname_u1` and `varname_u2`
249  on the boundaries corresponding to `region1` and `region2`.
250  The penalization parameter `dataname_r` should be chosen
251  large enough to prescribe approximate non-penetration and friction
252  conditions but not too large not to deteriorate too much the
253  conditionning of the tangent system.
254  The parameter `dataname_friction_coeffs` contains the Coulomb friction
255  coefficient and optionally an adhesional shear stress threshold and the
256  tresca limit shear stress. For constant coefficients its size is from
257  1 to 3. For coefficients described on a finite element method on the
258  same mesh as `varname_u1`, this vector contains a number of single
259  values, value pairs or triplets equal to the number of the
260  corresponding mesh_fem's basic dofs.
261  `dataname_lambda` is an optional parameter used if option
262  is 2. In that case, the penalization term is shifted by lambda (this
263  allows the use of an Uzawa algorithm on the corresponding augmented
264  Lagrangian formulation)
265  `dataname_alpha`, `dataname_wt1` and `dataname_wt2` are optional parameters
266  to solve evolutionary friction problems.
267  */
269  (model &md, const mesh_im &mim, const std::string &varname_u1,
270  const std::string &varname_u2, const std::string &dataname_r,
271  const std::string &dataname_friction_coeffs,
272  size_type region1, size_type region2, int option = 1,
273  const std::string &dataname_lambda = "",
274  const std::string &dataname_alpha = "",
275  const std::string &dataname_wt1 = "",
276  const std::string &dataname_wt2 = "");
277 
278 
279  enum contact_nonlinear_term_version { RHS_L_V1,
280  RHS_L_V2,
281  K_LL_V1,
282  K_LL_V2,
283  UZAWA_PROJ,
284  CONTACT_FLAG,
285  CONTACT_PRESSURE,
286 
287  RHS_U_V1,
288  RHS_U_V2,
289  RHS_U_V4,
290  RHS_U_V5,
291  RHS_U_FRICT_V1,
292  RHS_U_FRICT_V4,
293  RHS_U_FRICT_V6,
294  RHS_U_FRICT_V7,
295  RHS_U_FRICT_V8,
296  RHS_U_FRICT_V5,
297  RHS_L_FRICT_V1,
298  RHS_L_FRICT_V2,
299  RHS_L_FRICT_V4,
300  K_UL_V1,
301  K_UL_V2,
302  K_UL_V3,
303  UZAWA_PROJ_FRICT,
304  UZAWA_PROJ_FRICT_SAXCE,
305 
306  K_UU_V1,
307  K_UU_V2,
308  K_UL_FRICT_V1, // negative EYE
309  K_UL_FRICT_V2,
310  K_UL_FRICT_V3,
311  K_UL_FRICT_V4,
312  K_UL_FRICT_V5,
313  K_UL_FRICT_V7,
314  K_UL_FRICT_V8,
315  K_LL_FRICT_V1,
316  K_LL_FRICT_V2,
317  K_LL_FRICT_V3,
318  K_LL_FRICT_V4,
319  K_UU_FRICT_V1,
320  K_UU_FRICT_V2,
321  K_UU_FRICT_V3,
322  K_UU_FRICT_V4,
323  K_UU_FRICT_V5
324  };
325 
326  class contact_nonlinear_term : public nonlinear_elem_term {
327 
328  protected:
329  // the following variables are used in the compute method and their values
330  // have to be calculated during the preceding calls to the prepare method
331  // at the current interpolation context
332  base_small_vector lnt, lt; // multiplier lambda and its tangential component lambda_t
333  scalar_type ln; // normal component lambda_n of the multiplier
334  base_small_vector zt; // tangential relative displacement
335  scalar_type un; // normal relative displacement (positive when the first
336  // elastic body surface moves outwards)
337  base_small_vector no; // surface normal, pointing outwards with respect
338  // to the (first) elastic body
339  scalar_type g; // gap value
340  scalar_type f_coeff; // coefficient of friction value
341  scalar_type tau_adh; // adhesional shear resistance of the interface
342  scalar_type tresca_lim; // tresca shear limit for the interface
343 
344  // these variables are used as temporary storage and they will usually contain
345  // garbage from previous calculations
346  base_small_vector aux1, auxN, V; // helper vectors of size 1, N and N respectively
347  base_matrix GP; // helper matrix of size NxN
348 
349  void adjust_tensor_size(void);
350 
351  public:
352  dim_type N;
353  size_type option;
354  scalar_type r;
355  bool contact_only;
356  scalar_type alpha;
357 
358  bgeot::multi_index sizes_;
359 
360  contact_nonlinear_term(dim_type N_, size_type option_, scalar_type r_,
361  bool contact_only_ = true,
362  scalar_type alpha_ = scalar_type(1)) :
363  tau_adh(0), tresca_lim(gmm::default_max(scalar_type())),
364  N(N_), option(option_), r(r_), contact_only(contact_only_), alpha(alpha_) {
365 
366  adjust_tensor_size();
367  }
368 
369  const bgeot::multi_index &sizes(size_type) const { return sizes_; }
370 
371  virtual void friction_law(scalar_type p, scalar_type &tau);
372  virtual void friction_law(scalar_type p, scalar_type &tau, scalar_type &tau_grad);
373 
374  virtual void compute(fem_interpolation_context&, bgeot::base_tensor &t);
375  virtual void prepare(fem_interpolation_context& /*ctx*/, size_type /*nb*/)
376  { GMM_ASSERT1(false, "the prepare method has to be reimplemented in "
377  "a derived class"); }
378 
379  };
380 
381 
382  class contact_rigid_obstacle_nonlinear_term : public contact_nonlinear_term {
383 
384  // temporary variables to be used inside the prepare method
385  base_small_vector vt; // of size N
386  base_vector coeff; // of variable size
387  base_matrix grad; // of size 1xN
388 
389  public:
390  // class specific objects to take into account inside the prepare method
391  const mesh_fem &mf_u; // mandatory
392  const mesh_fem &mf_obs; // mandatory
393  const mesh_fem *pmf_lambda; // optional for terms involving lagrange multipliers
394  const mesh_fem *pmf_coeff; // optional for terms involving fem described coefficient of friction
395  base_vector U, obs, lambda, friction_coeff, tau_adhesion, tresca_limit, WT, VT;
396  scalar_type gamma;
397 
398  template <typename VECT1>
399  contact_rigid_obstacle_nonlinear_term
400  (size_type option_, scalar_type r_,
401  const mesh_fem &mf_u_, const VECT1 &U_,
402  const mesh_fem &mf_obs_, const VECT1 &obs_,
403  const mesh_fem *pmf_lambda_ = 0, const VECT1 *lambda_ = 0,
404  const mesh_fem *pmf_coeff_ = 0, const VECT1 *f_coeffs_ = 0,
405  scalar_type alpha_ = scalar_type(1), const VECT1 *WT_ = 0,
406  scalar_type gamma_ = scalar_type(1), const VECT1 *VT_ = 0
407  )
408  : contact_nonlinear_term(mf_u_.linked_mesh().dim(), option_, r_,
409  (f_coeffs_ == 0), alpha_
410  ),
411  mf_u(mf_u_), mf_obs(mf_obs_),
412  pmf_lambda(pmf_lambda_), pmf_coeff(pmf_coeff_),
413  U(mf_u.nb_basic_dof()), obs(mf_obs.nb_basic_dof()),
414  lambda(0), friction_coeff(0), tau_adhesion(0), tresca_limit(0),
415  WT(0), VT(0), gamma(gamma_)
416  {
417 
418  mf_u.extend_vector(U_, U);
419  mf_obs.extend_vector(obs_, obs);
420 
421  if (pmf_lambda) {
422  lambda.resize(pmf_lambda->nb_basic_dof());
423  pmf_lambda->extend_vector(*lambda_, lambda);
424  }
425 
426  if (!contact_only) {
427  if (!pmf_coeff) {
428  f_coeff = (*f_coeffs_)[0];
429  if (gmm::vect_size(*f_coeffs_) > 1) tau_adh = (*f_coeffs_)[1];
430  if (gmm::vect_size(*f_coeffs_) > 2) tresca_lim = (*f_coeffs_)[2];
431  }
432  else {
433  size_type ncoeffs = gmm::vect_size(*f_coeffs_)/pmf_coeff->nb_dof();
434  GMM_ASSERT1(ncoeffs >= 1 && ncoeffs <= 3, "Wrong vector dimension for friction coefficients");
435  gmm::resize(friction_coeff, pmf_coeff->nb_basic_dof());
436  pmf_coeff->extend_vector(gmm::sub_vector(*f_coeffs_, gmm::sub_slice(0,pmf_coeff->nb_dof(),ncoeffs)),
437  friction_coeff);
438  if (ncoeffs > 1) {
439  gmm::resize(tau_adhesion, pmf_coeff->nb_basic_dof());
440  pmf_coeff->extend_vector(gmm::sub_vector(*f_coeffs_, gmm::sub_slice(1,pmf_coeff->nb_dof(),ncoeffs)),
441  tau_adhesion);
442  }
443  if (ncoeffs > 2) {
444  gmm::resize(tresca_limit, pmf_coeff->nb_basic_dof());
445  pmf_coeff->extend_vector(gmm::sub_vector(*f_coeffs_, gmm::sub_slice(2,pmf_coeff->nb_dof(),ncoeffs)),
446  tresca_limit);
447  }
448  }
449 
450  if (WT_ && gmm::vect_size(*WT_)) {
451  WT.resize(mf_u.nb_basic_dof());
452  mf_u.extend_vector(*WT_, WT);
453  }
454 
455  if (VT_ && gmm::vect_size(*VT_)) {
456  VT.resize(mf_u.nb_basic_dof());
457  mf_u.extend_vector(*VT_, VT);
458  }
459  }
460 
461  vt.resize(N);
462  gmm::resize(grad, 1, N);
463 
464  GMM_ASSERT1(mf_u.get_qdim() == N, "wrong qdim for the mesh_fem");
465  }
466 
467  // this methode prepares all necessary data for the compute method
468  // of the base class
469  virtual void prepare(fem_interpolation_context& ctx, size_type nb);
470 
471  };
472 
473 
474  class contact_nonmatching_meshes_nonlinear_term : public contact_nonlinear_term {
475 
476  // temporary variables to be used inside the prepare method
477  base_vector coeff; // of variable size
478  base_matrix grad; // of size 1xN
479 
480  public:
481  const mesh_fem &mf_u1; // displacements on the non-mortar side
482  const mesh_fem &mf_u2; // displacements of the mortar side projected on the non-mortar side
483  const mesh_fem *pmf_lambda; // Lagrange multipliers defined on the non-mortar side
484  const mesh_fem *pmf_coeff; // coefficient of friction defined on the non-mortar side
485  base_vector U1, U2, lambda, friction_coeff, tau_adhesion, tresca_limit, WT1, WT2;
486 
487  template <typename VECT1>
488  contact_nonmatching_meshes_nonlinear_term
489  (size_type option_, scalar_type r_,
490  const mesh_fem &mf_u1_, const VECT1 &U1_,
491  const mesh_fem &mf_u2_, const VECT1 &U2_,
492  const mesh_fem *pmf_lambda_ = 0, const VECT1 *lambda_ = 0,
493  const mesh_fem *pmf_coeff_ = 0, const VECT1 *f_coeffs_ = 0,
494  scalar_type alpha_ = scalar_type(1),
495  const VECT1 *WT1_ = 0, const VECT1 *WT2_ = 0
496  )
497  : contact_nonlinear_term(mf_u1_.linked_mesh().dim(), option_, r_,
498  (f_coeffs_ == 0), alpha_
499  ),
500  mf_u1(mf_u1_), mf_u2(mf_u2_),
501  pmf_lambda(pmf_lambda_), pmf_coeff(pmf_coeff_),
502  U1(mf_u1.nb_basic_dof()), U2(mf_u2.nb_basic_dof()),
503  lambda(0), friction_coeff(0), tau_adhesion(0), tresca_limit(0),
504  WT1(0), WT2(0)
505  {
506 
507  GMM_ASSERT1(mf_u2.linked_mesh().dim() == N,
508  "incompatible mesh dimensions for the given mesh_fem's");
509 
510  mf_u1.extend_vector(U1_, U1);
511  mf_u2.extend_vector(U2_, U2);
512 
513  if (pmf_lambda) {
514  lambda.resize(pmf_lambda->nb_basic_dof());
515  pmf_lambda->extend_vector(*lambda_, lambda);
516  }
517 
518  if (!contact_only) {
519  if (!pmf_coeff) {
520  f_coeff = (*f_coeffs_)[0];
521  if (gmm::vect_size(*f_coeffs_) > 1) tau_adh = (*f_coeffs_)[1];
522  if (gmm::vect_size(*f_coeffs_) > 2) tresca_lim = (*f_coeffs_)[2];
523  }
524  else {
525  size_type ncoeffs = gmm::vect_size(*f_coeffs_)/pmf_coeff->nb_dof();
526  GMM_ASSERT1(ncoeffs >= 1 && ncoeffs <= 3, "Wrong vector dimension for friction coefficients");
527  gmm::resize(friction_coeff, pmf_coeff->nb_basic_dof());
528  pmf_coeff->extend_vector(gmm::sub_vector(*f_coeffs_, gmm::sub_slice(0,pmf_coeff->nb_dof(),ncoeffs)),
529  friction_coeff);
530  if (ncoeffs > 1) {
531  gmm::resize(tau_adhesion, pmf_coeff->nb_basic_dof());
532  pmf_coeff->extend_vector(gmm::sub_vector(*f_coeffs_, gmm::sub_slice(1,pmf_coeff->nb_dof(),ncoeffs)),
533  tau_adhesion);
534  }
535  if (ncoeffs > 2) {
536  gmm::resize(tresca_limit, pmf_coeff->nb_basic_dof());
537  pmf_coeff->extend_vector(gmm::sub_vector(*f_coeffs_, gmm::sub_slice(2,pmf_coeff->nb_dof(),ncoeffs)),
538  tresca_limit);
539  }
540  }
541  if (WT1_ && WT2_ && gmm::vect_size(*WT1_) && gmm::vect_size(*WT2_)) {
542  WT1.resize(mf_u1.nb_basic_dof());
543  mf_u1.extend_vector(*WT1_, WT1);
544  WT2.resize(mf_u2.nb_basic_dof());
545  mf_u2.extend_vector(*WT2_, WT2);
546  }
547  }
548 
549  gmm::resize(grad, 1, N);
550 
551  GMM_ASSERT1(mf_u1.get_qdim() == N, "wrong qdim for the 1st mesh_fem");
552  GMM_ASSERT1(mf_u2.get_qdim() == N, "wrong qdim for the 2nd mesh_fem");
553  }
554 
555  // this method prepares all necessary data for the compute method
556  // of the base class
557  virtual void prepare(fem_interpolation_context& ctx, size_type nb);
558 
559  };
560 
561 
562  /** Specific assembly procedure for the use of an Uzawa algorithm to solve
563  contact with rigid obstacle problems with friction.
564  */
565  template<typename VECT1>
567  (VECT1 &R, const mesh_im &mim,
568  const getfem::mesh_fem &mf_u, const VECT1 &U,
569  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
570  const getfem::mesh_fem &mf_lambda, const VECT1 &lambda,
571  const getfem::mesh_fem *pmf_coeff, const VECT1 &f_coeff, const VECT1 *WT,
572  scalar_type r, scalar_type alpha, const mesh_region &rg, int option = 1) {
573 
574  contact_rigid_obstacle_nonlinear_term
575  nterm1((option == 1) ? UZAWA_PROJ_FRICT : UZAWA_PROJ_FRICT_SAXCE, r,
576  mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
577  pmf_coeff, &f_coeff, alpha, WT);
578 
580  if (pmf_coeff) // variable coefficient of friction
581  assem.set("V(#3)+=comp(NonLin$1(#1,#1,#2,#3,#4).vBase(#3))(i,:,i); ");
582  else // constant coefficient of friction
583  assem.set("V(#3)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#3))(i,:,i); ");
584  assem.push_mi(mim);
585  assem.push_mf(mf_u);
586  assem.push_mf(mf_obs);
587  assem.push_mf(mf_lambda);
588  if (pmf_coeff)
589  assem.push_mf(*pmf_coeff);
590  assem.push_nonlinear_term(&nterm1);
591  assem.push_vec(R);
592  assem.assembly(rg);
593  }
594 
595  /** Specific assembly procedure for the use of an Uzawa algorithm to solve
596  contact problems.
597  */
598  template<typename VECT1>
600  (VECT1 &R, const mesh_im &mim,
601  const getfem::mesh_fem &mf_u, const VECT1 &U,
602  const getfem::mesh_fem &mf_obs, const VECT1 &obs,
603  const getfem::mesh_fem &mf_lambda, const VECT1 &lambda,
604  scalar_type r, const mesh_region &rg) {
605 
606  contact_rigid_obstacle_nonlinear_term
607  nterm1(UZAWA_PROJ, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda);
608 
610  assem.set("V(#3)+=comp(NonLin$1(#1,#1,#2,#3).Base(#3))(i,:); ");
611  assem.push_mi(mim);
612  assem.push_mf(mf_u);
613  assem.push_mf(mf_obs);
614  assem.push_mf(mf_lambda);
615  assem.push_nonlinear_term(&nterm1);
616  assem.push_vec(R);
617  assem.assembly(rg);
618  }
619 
620 
621  /** Specific assembly procedure for the use of an Uzawa algorithm to solve
622  contact problems.
623  */
624  template<typename VEC>
626  (VEC &R, const mesh_im &mim,
627  const getfem::mesh_fem &mf_u,
628  const getfem::mesh_fem &mf_obs, const VEC &obs,
629  const getfem::mesh_fem &mf_lambda, const VEC &lambda,
630  const mesh_region &rg) {
631 
632  bool contact_only = (mf_lambda.get_qdim() == 1);
633 
634  VEC U;
635  gmm::resize(U, mf_u.nb_dof());
636  scalar_type dummy_r(0.);
637  VEC dummy_f_coeff; gmm::resize(dummy_f_coeff,1);
638  contact_rigid_obstacle_nonlinear_term
639  nterm(RHS_U_V1, dummy_r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
640  0, contact_only ? 0 : &dummy_f_coeff);
641 
643  assem.set("V(#1)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#1))(i,:,i); ");
644  assem.push_mi(mim);
645  assem.push_mf(mf_u);
646  assem.push_mf(mf_obs);
647  assem.push_mf(mf_lambda);
648  assem.push_nonlinear_term(&nterm);
649  assem.push_vec(R);
650  assem.assembly(rg);
651  }
652 
653  template<typename VEC>
654  scalar_type asm_level_set_contact_area
655  (const mesh_im &mim,
656  const getfem::mesh_fem &mf_u, const VEC &U,
657  const getfem::mesh_fem &mf_obs, const VEC &obs,
658  const mesh_region &rg, scalar_type threshold_factor=0.0,
659  const getfem::mesh_fem *mf_lambda=0, const VEC *lambda=0,
660  scalar_type threshold_pressure_factor=0.0) {
661 
662  if (!rg.get_parent_mesh())
663  rg.from_mesh(mim.linked_mesh());
664  getfem::mesh_fem mf_ca(mf_u.linked_mesh());
665  mf_ca.set_classical_finite_element(rg.index(),1);
666 
667  VEC mesh_size(mf_ca.nb_dof());
668  VEC mesh_size2(mf_ca.nb_dof());
669  { // assemble an estimator of the mesh size
670  getfem::generic_assembly assem_mesh_size;
671  assem_mesh_size.set("V(#1)+=comp(Base(#1))");
672  assem_mesh_size.push_mi(mim);
673  assem_mesh_size.push_mf(mf_ca);
674  assem_mesh_size.push_vec(mesh_size2);
675  assem_mesh_size.assembly(rg);
676  for (dal::bv_visitor_c dof(mf_ca.basic_dof_on_region(rg));
677  !dof.finished(); ++dof)
678  mesh_size[dof] = sqrt(mesh_size2[dof]);
679  }
680 
681  VEC threshold(mf_ca.nb_dof());
682  if (mf_lambda && lambda) {
683  VEC pressure(mf_ca.nb_dof());
684  VEC dummy_f_coeff(1);
685  bool contact_only = (mf_lambda->get_qdim() == 1);
686  contact_rigid_obstacle_nonlinear_term
687  nterm_pressure(CONTACT_PRESSURE, 0., mf_u, U, mf_obs, obs,
688  mf_lambda, lambda, 0, contact_only ? 0 : &dummy_f_coeff);
689 
690  getfem::generic_assembly assem_pressure;
691  assem_pressure.set("V(#4)+=comp(NonLin(#1,#1,#2,#3).Base(#4))(i,:)");
692  assem_pressure.push_mi(mim);
693  assem_pressure.push_mf(mf_u);
694  assem_pressure.push_mf(mf_obs);
695  assem_pressure.push_mf(*mf_lambda);
696  assem_pressure.push_mf(mf_ca);
697  assem_pressure.push_nonlinear_term(&nterm_pressure);
698  assem_pressure.push_vec(pressure);
699  assem_pressure.assembly(rg);
700  for (dal::bv_visitor_c dof(mf_ca.basic_dof_on_region(rg));
701  !dof.finished(); ++dof)
702  pressure[dof] /= mesh_size2[dof];
703 
704  // in areas where pressure is clearly non-zero set a low threshold
705  // in order to avoid false negative contact detection
706  scalar_type threshold_pressure(threshold_pressure_factor *
707  gmm::vect_norminf(pressure));
708  gmm::copy(gmm::scaled(mesh_size, scalar_type(-1)), threshold);
709  for (getfem::mr_visitor v(rg); !v.finished(); ++v) {
710  size_type nbdof = mf_ca.nb_basic_dof_of_face_of_element(v.cv(), v.f());
711  mesh_fem::ind_dof_face_ct::const_iterator
712  itdof = mf_ca.ind_basic_dof_of_face_of_element(v.cv(), v.f()).begin();
713  bool all_positive = true;
714  for (size_type k=0; k < nbdof; ++k, ++itdof)
715  if (pressure[*itdof] < threshold_pressure) { all_positive = false; break; }
716  if (!all_positive) {
717  itdof = mf_ca.ind_basic_dof_of_face_of_element(v.cv(), v.f()).begin();
718  for (size_type k=0; k < nbdof; ++k, ++itdof)
719  threshold[*itdof] = threshold_factor * mesh_size[*itdof];
720  }
721  }
722  }
723  else
724  gmm::copy(gmm::scaled(mesh_size, threshold_factor), threshold);
725 
726  // compute the total contact area
727  // remark: the CONTACT_FLAG option misuses lambda as a threshold field
728  contact_rigid_obstacle_nonlinear_term
729  nterm(CONTACT_FLAG, 0., mf_u, U, mf_obs, obs, &mf_ca, &threshold);
730 
732  assem.set("V()+=comp(NonLin(#1,#1,#2,#3))(i)");
733  assem.push_mi(mim);
734  assem.push_mf(mf_u);
735  assem.push_mf(mf_obs);
736  assem.push_mf(mf_ca);
737  assem.push_nonlinear_term(&nterm);
738  std::vector<scalar_type> v(1);
739  assem.push_vec(v);
740  assem.assembly(rg);
741  return v[0];
742  }
743 
744 
745  template<typename VEC>
746  void asm_nonmatching_meshes_normal_source_term
747  (VEC &R, const mesh_im &mim,
748  const getfem::mesh_fem &mf_u1,
749  const getfem::mesh_fem &mf_u2_proj,
750  const getfem::mesh_fem &mf_lambda, const VEC &lambda,
751  const mesh_region &rg) {
752 
753  bool contact_only = (mf_lambda.get_qdim() == 1);
754 
755  VEC U1, U2_proj;
756  gmm::resize(U1, mf_u1.nb_dof());
757  gmm::resize(U2_proj, mf_u2_proj.nb_dof());
758  scalar_type dummy_r(0);
759  VEC dummy_f_coeff; gmm::resize(dummy_f_coeff,1);
760  contact_nonmatching_meshes_nonlinear_term
761  nterm(RHS_U_V1, dummy_r, mf_u1, U1, mf_u2_proj, U2_proj, &mf_lambda, &lambda,
762  0, contact_only ? 0 : &dummy_f_coeff);
763 
765  assem.set("V(#1)+=comp(NonLin(#1,#1,#2,#3).vBase(#1))(i,:,i)");
766  assem.push_mi(mim);
767  assem.push_mf(mf_u1);
768  assem.push_mf(mf_u2_proj);
769  assem.push_mf(mf_lambda);
770  assem.push_nonlinear_term(&nterm);
771  assem.push_vec(R);
772  assem.assembly(rg);
773  }
774 
775  template<typename VEC>
776  scalar_type asm_nonmatching_meshes_contact_area
777  (const mesh_im &mim,
778  const getfem::mesh_fem &mf_u1, const VEC &U1,
779  const getfem::mesh_fem &mf_u2_proj, const VEC &U2_proj,
780  const mesh_region &rg, scalar_type threshold_factor=0.0,
781  const getfem::mesh_fem *mf_lambda=0, const VEC *lambda=0,
782  scalar_type threshold_pressure_factor=0.0) {
783 
784  if (!rg.get_parent_mesh())
785  rg.from_mesh(mim.linked_mesh());
786  getfem::mesh_fem mf_ca(mf_u1.linked_mesh());
787  mf_ca.set_classical_finite_element(rg.index(),1);
788 
789  VEC mesh_size(mf_ca.nb_dof());
790  VEC mesh_size2(mf_ca.nb_dof());
791  { // assemble an estimator of the mesh size
792  getfem::generic_assembly assem_mesh_size;
793  assem_mesh_size.set("V(#1)+=comp(Base(#1))");
794  assem_mesh_size.push_mi(mim);
795  assem_mesh_size.push_mf(mf_ca);
796  assem_mesh_size.push_vec(mesh_size2);
797  assem_mesh_size.assembly(rg);
798  for (dal::bv_visitor_c dof(mf_ca.basic_dof_on_region(rg));
799  !dof.finished(); ++dof)
800  mesh_size[dof] = sqrt(mesh_size2[dof]);
801  }
802 
803  VEC threshold(mf_ca.nb_dof());
804  if (mf_lambda && lambda) {
805  VEC pressure(mf_ca.nb_dof());
806  VEC dummy_f_coeff(1);
807  bool contact_only = (mf_lambda->get_qdim() == 1);
808  contact_nonmatching_meshes_nonlinear_term
809  nterm_pressure(CONTACT_PRESSURE, 0., mf_u1, U1, mf_u2_proj, U2_proj,
810  mf_lambda, lambda, 0, contact_only ? 0 : &dummy_f_coeff);
811 
812  getfem::generic_assembly assem_pressure;
813  assem_pressure.set("V(#4)+=comp(NonLin(#1,#1,#2,#3).Base(#4))(i,:)");
814  assem_pressure.push_mi(mim);
815  assem_pressure.push_mf(mf_u1);
816  assem_pressure.push_mf(mf_u2_proj);
817  assem_pressure.push_mf(*mf_lambda);
818  assem_pressure.push_mf(mf_ca);
819  assem_pressure.push_nonlinear_term(&nterm_pressure);
820  assem_pressure.push_vec(pressure);
821  assem_pressure.assembly(rg);
822  for (dal::bv_visitor_c dof(mf_ca.basic_dof_on_region(rg));
823  !dof.finished(); ++dof)
824  pressure[dof] /= mesh_size2[dof];
825 
826  // in areas where pressure is clearly non-zero set a low threshold
827  // in order to avoid false negative contact detection
828  scalar_type threshold_pressure(threshold_pressure_factor *
829  gmm::vect_norminf(pressure));
830  gmm::copy(gmm::scaled(mesh_size, scalar_type(-1)), threshold);
831  for (getfem::mr_visitor v(rg); !v.finished(); ++v) {
832  size_type nbdof = mf_ca.nb_basic_dof_of_face_of_element(v.cv(), v.f());
833  mesh_fem::ind_dof_face_ct::const_iterator
834  itdof = mf_ca.ind_basic_dof_of_face_of_element(v.cv(), v.f()).begin();
835  bool all_positive = true;
836  for (size_type k=0; k < nbdof; ++k, ++itdof)
837  if (pressure[*itdof] < threshold_pressure) { all_positive = false; break; }
838  if (!all_positive) {
839  itdof = mf_ca.ind_basic_dof_of_face_of_element(v.cv(), v.f()).begin();
840  for (size_type k=0; k < nbdof; ++k, ++itdof)
841  threshold[*itdof] = threshold_factor * mesh_size[*itdof];
842  }
843  }
844  }
845  else
846  gmm::copy(gmm::scaled(mesh_size, threshold_factor), threshold);
847 
848 
849  // compute the total contact area
850  // remark: the CONTACT_FLAG option misuses lambda as a threshold field
851  contact_nonmatching_meshes_nonlinear_term
852  nterm(CONTACT_FLAG, 0., mf_u1, U1, mf_u2_proj, U2_proj, &mf_ca, &threshold);
853 
855  assem.set("V()+=comp(NonLin(#1,#1,#2,#3))(i)");
856  assem.push_mi(mim);
857  assem.push_mf(mf_u1);
858  assem.push_mf(mf_u2_proj);
859  assem.push_mf(mf_ca);
860  assem.push_nonlinear_term(&nterm);
861  std::vector<scalar_type> v(1);
862  assem.push_vec(v);
863  assem.assembly(rg);
864  return v[0];
865  }
866 
867 
868  void compute_integral_contact_area_and_force
869  (model &md, size_type indbrick,
870  scalar_type &area, model_real_plain_vector &Forces);
871 
872 
873 
874  /** Adds a contact condition with or without Coulomb friction on the variable
875  `varname_u` and the mesh boundary `region`. `Neumannterm`
876  is the expression of the Neumann term (obtained by the Green formula)
877  described as an expression of the high-level
878  generic assembly language. This term can be obtained with
879  md.Neumann_term(varname, region) once all volumic bricks have
880  been added to the model. The contact condition
881  is prescribed with Nitsche's method. The rigid obstacle should
882  be described with the data `dataname_obstacle` being a signed distance to
883  the obstacle (interpolated on a finite element method).
884  `gamma0name` is the Nitsche's method parameter.
885  `theta` is a scalar value which can be
886  positive or negative. `theta = 1` corresponds to the standard symmetric
887  method which is conditionnaly coercive for `gamma0` small.
888  `theta = -1` corresponds to the skew-symmetric method which is
889  inconditionnaly coercive. `theta = 0` is the simplest method
890  for which the second derivative of the Neumann term is not necessary.
891  The optional parameter `dataexpr_friction_coeff` is the friction
892  coefficient which could be any expression of the assembly language.
893  Returns the brick index in the model.
894  */
896  (model &md, const mesh_im &mim, const std::string &varname_u,
897  const std::string &Neumannterm,
898  const std::string &expr_obs, const std::string &dataname_gamma0,
899  scalar_type theta_,
900  std::string dataexpr_friction_coeff,
901  const std::string &dataname_alpha,
902  const std::string &dataname_wt,
903  size_type region);
904 
905 #ifdef EXPERIMENTAL_PURPOSE_ONLY
906 
907  /** Adds a contact condition with or without Coulomb friction on the variable
908  `varname_u` and the mesh boundary `region`. The contact condition
909  is prescribed with Nitsche's method. The rigid obstacle should
910  be described with the data `dataname_obstacle` being a signed distance to
911  the obstacle (interpolated on a finite element method).
912  `gamma0name` is the Nitsche's method parameter.
913  `theta` is a scalar value which can be
914  positive or negative. `theta = 1` corresponds to the standard symmetric
915  method which is conditionnaly coercive for `gamma0` small.
916  `theta = -1` corresponds to the skew-symmetric method which is
917  inconditionnaly coercive. `theta = 0` is the simplest method
918  for which the second derivative of the Neumann term is not necessary.
919  The optional parameter `dataname_friction_coeff` is the friction
920  coefficient which could be constant or defined on a finite element
921  method.
922  Returns the brick index in the model.
923  */
924  size_type add_Nitsche_contact_with_rigid_obstacle_brick_modified_midpoint
925  (model &md, const mesh_im &mim, const std::string &varname_u,
926  const std::string &Neumannterm, const std::string &Neumannterm_wt,
927  const std::string &obs, const std::string &dataname_gamma0,
928  scalar_type theta_,
929  std::string dataname_friction_coeff,
930  const std::string &dataname_alpha,
931  const std::string &dataname_wt,
932  size_type region);
933 
934 #endif
935 
936 
937 
938  /** Adds a contact condition with or without Coulomb friction between
939  two bodies in a fictitious domain. The contact condition is applied on
940  the variable `varname_u1` corresponds with the first
941  and slave body with Nitsche's method and on the variable `varname_u2`
942  corresponds with the second and master body with Nitsche's method.
943  The contact condition is evaluated on the fictitious slave bondary.
944  The first body should be described by the level-set `dataname_d1`
945  and the second body should be described by the level-set
946 `dataname_d2`. `gamma0name` is the Nitsche's method parameter.
947  `theta` is a scalar value which can be positive or negative.
948  `theta = 1` corresponds to the standard symmetric method which
949  is conditionnaly coercive for `gamma0` small.
950  `theta = -1` corresponds to the skew-symmetric method which is
951  inconditionnaly coercive. `theta = 0` is the simplest method for
952  which the second derivative of the Neumann term is not necessary.
953 The optional parameter `dataname_friction_coeff` is the friction
954 coefficient which could be constant or defined on a finite element method.
955 CAUTION: This brick has to be added in the model
956  after all the bricks corresponding to partial differential
957  terms having a Neumann term. Moreover, This brick can only
958  be applied to bricks declaring their Neumann terms. Returns the brick index in the model.
959 
960  */
962  (model &md, const mesh_im &mim, const std::string &varname_u1,
963  const std::string &varname_u2, const std::string &dataname_d1,
964  const std::string &dataname_d2, const std::string &dataname_gamma0,
965  scalar_type theta,
966  const std::string &dataname_friction_coeff,
967  const std::string &dataname_alpha,
968  const std::string &dataname_wt1, const std::string &dataname_wt2);
969 
970 
971 } /* end of namespace getfem. */
972 
973 
974 #endif /* GETFEM_CONTACT_AND_FRICTION_INTEGRAL_H__ */
Generic assembly of vectors, matrices.
void push_vec(VEC &v)
Add a new output vector.
void assembly(const mesh_region &region=mesh_region::all_convexes())
do the assembly on the specified region (boundary or set of convexes)
void push_nonlinear_term(pnonlinear_elem_term net)
Add a new non-linear term.
void push_mi(const mesh_im &im_)
Add a new mesh_im.
void push_mf(const mesh_fem &mf_)
Add a new mesh_fem.
Describe a finite element method linked to a mesh.
virtual dim_type get_qdim() const
Return the Q dimension.
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Describe an integration method linked to a mesh.
"iterator" class for regions.
structure used to hold a set of convexes and/or convex faces.
`‘Model’' variables store the variables, the data and the description of a model.
Generic assembly implementation.
Model representation in Getfem.
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:209
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
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
GEneric Tool for Finite Element Methods.
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.
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...
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.
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.
size_type add_Nitsche_fictitious_domain_contact_brick(model &md, const mesh_im &mim, const std::string &varname_u1, const std::string &varname_u2, const std::string &dataname_d1, const std::string &dataname_d2, const std::string &dataname_gamma0, scalar_type theta, const std::string &dataname_friction_coeff, const std::string &dataname_alpha, const std::string &dataname_wt1, const std::string &dataname_wt2)
Adds a contact condition with or without Coulomb friction between two bodies in a fictitious domain.
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.
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...
void asm_integral_contact_Uzawa_proj(VECT1 &R, const mesh_im &mim, const getfem::mesh_fem &mf_u, const VECT1 &U, const getfem::mesh_fem &mf_obs, const VECT1 &obs, const getfem::mesh_fem &mf_lambda, const VECT1 &lambda, const getfem::mesh_fem *pmf_coeff, const VECT1 &f_coeff, const VECT1 *WT, scalar_type r, scalar_type alpha, const mesh_region &rg, int option=1)
Specific assembly procedure for the use of an Uzawa algorithm to solve contact with rigid obstacle pr...