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