GetFEM  5.4.2
getfem_contact_and_friction_large_sliding.cc
1 /*===========================================================================
2 
3  Copyright (C) 2013-2020 Yves Renard, Konstantinos Poulios.
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 #include "getfem/bgeot_rtree.h"
23 #include "getfem/getfem_contact_and_friction_integral.h"
24 #include "getfem/getfem_contact_and_friction_common.h"
25 #include "getfem/getfem_contact_and_friction_large_sliding.h"
26 #include "getfem/getfem_assembling.h"
27 #include "gmm/gmm_condition_number.h"
28 
29 namespace getfem {
30 
31 
32  //=========================================================================
33  // Augmented friction law
34  //=========================================================================
35 
36 
37 #define FRICTION_LAW 1
38 
39 
40 #if FRICTION_LAW == 1 // Complete law with friction
41 
42  template <typename VEC, typename VEC2, typename VECR>
43  void aug_friction(const VEC &lambda, scalar_type g, const VEC &Vs,
44  const VEC &n, scalar_type r, const VEC2 &f, VECR &F) {
45  scalar_type nn = gmm::vect_norm2(n);
46  scalar_type lambdan = gmm::vect_sp(lambda, n)/nn;
47  scalar_type lambdan_aug = gmm::neg(lambdan + r * g);
48  size_type i = gmm::vect_size(f);
49  scalar_type tau = ((i >= 3) ? f[2] : scalar_type(0)) + f[0]*lambdan_aug;
50  if (i >= 2) tau = std::min(tau, f[1]);
51 
52  if (tau > scalar_type(0)) {
53  gmm::add(lambda, gmm::scaled(Vs, -r), F);
54  scalar_type mu = gmm::vect_sp(F, n)/nn;
55  gmm::add(gmm::scaled(n, -mu/nn), F);
56  scalar_type norm = gmm::vect_norm2(F);
57  if (norm > tau) gmm::scale(F, tau / norm);
58  } else { gmm::clear(F); }
59 
60  gmm::add(gmm::scaled(n, -lambdan_aug/nn), F);
61  }
62 
63  template <typename VEC, typename VEC2, typename VECR, typename MAT>
64  void aug_friction_grad(const VEC &lambda, scalar_type g, const VEC &Vs,
65  const VEC &n, scalar_type r, const VEC2 &f, VECR &F,
66  MAT &dlambda, VECR &dg, MAT &dn, MAT &dVs) {
67  size_type N = gmm::vect_size(lambda);
68  scalar_type nn = gmm::vect_norm2(n);
69  scalar_type lambdan = gmm::vect_sp(lambda, n)/nn;
70  scalar_type lambdan_aug = gmm::neg(lambdan + r * g);
71  size_type i = gmm::vect_size(f);
72  scalar_type tau = ((i >= 3) ? f[2] : scalar_type(0)) + f[0]*lambdan_aug;
73  if (i >= 2) tau = std::min(tau, f[1]);
74  scalar_type norm(0);
75 
76  if (tau > scalar_type(0)) {
77  gmm::add(lambda, gmm::scaled(Vs, -r), F);
78  scalar_type mu = gmm::vect_sp(F, n)/nn;
79  gmm::add(gmm::scaled(n, -mu/nn), F);
80  norm = gmm::vect_norm2(F);
81  gmm::copy(gmm::identity_matrix(), dn);
82  gmm::scale(dn, -mu/nn);
83  gmm::rank_one_update(dn, gmm::scaled(n, mu/(nn*nn*nn)), n);
84  gmm::rank_one_update(dn, gmm::scaled(n, scalar_type(-1)/(nn*nn)), F);
85  gmm::copy(gmm::identity_matrix(), dVs);
86  gmm::rank_one_update(dVs, n, gmm::scaled(n, scalar_type(-1)/(nn*nn)));
87 
88  if (norm > tau) {
89  gmm::rank_one_update(dVs, F,
90  gmm::scaled(F, scalar_type(-1)/(norm*norm)));
91  gmm::scale(dVs, tau / norm);
92  gmm::copy(gmm::scaled(F, scalar_type(1)/norm), dg);
93  gmm::rank_one_update(dn, gmm::scaled(F, mu/(norm*norm*nn)), F);
94  gmm::scale(dn, tau / norm);
95  gmm::scale(F, tau / norm);
96  } else gmm::clear(dg);
97 
98  } else { gmm::clear(dg); gmm::clear(dVs); gmm::clear(F); gmm::clear(dn); }
99  // At this stage, F = P_{B_T}, dVs = d_v P_{B_T}, dn = d_n P_{B_T}
100  // and dg = d_tau P_{B_T}.
101 
102  gmm::copy(dVs, dlambda);
103  if (norm > tau && ((i <= 1) || tau < f[1]) && ((i <= 2) || tau > f[2])) {
104  gmm::rank_one_update(dn, dg, gmm::scaled(lambda, -f[0]/nn));
105  gmm::rank_one_update(dn, dg, gmm::scaled(n, f[0]*lambdan/(nn*nn)));
106  gmm::rank_one_update(dlambda, dg, gmm::scaled(n, -f[0]/nn));
107  gmm::scale(dg, -f[0]*r);
108  } else gmm::clear(dg);
109  if (lambdan_aug > scalar_type(0)) {
110  gmm::add(gmm::scaled(n, r/nn), dg);
111  gmm::rank_one_update(dlambda, n, gmm::scaled(n, scalar_type(1)/(nn*nn)));
112  gmm::rank_one_update(dn, gmm::scaled(n, scalar_type(1)/(nn*nn)), lambda);
113  gmm::rank_one_update(dn,
114  gmm::scaled(n,(lambdan_aug-lambdan)/(nn*nn*nn)), n);
115  for (size_type j = 0; j < N; ++j) dn(j,j) -= lambdan_aug/nn;
116  }
117  gmm::add(gmm::scaled(n, -lambdan_aug/nn), F);
118 
119  gmm::scale(dVs, -r);
120  }
121 
122 #elif FRICTION_LAW == 2 // Contact only
123 
124  template <typename VEC, typename VEC2, typename VECR>
125  void aug_friction(const VEC &lambda, scalar_type g, const VEC &,
126  const VEC &n, scalar_type r, const VEC2 &, VECR &F) {
127  scalar_type nn = gmm::vect_norm2(n);
128  scalar_type lambdan = gmm::vect_sp(lambda, n)/nn;
129  scalar_type lambdan_aug = gmm::neg(lambdan + r * g);
130  gmm::copy(gmm::scaled(n, -lambdan_aug/nn), F);
131  }
132 
133  template <typename VEC, typename VEC2, typename VECR, typename MAT>
134  void aug_friction_grad(const VEC &lambda, scalar_type g, const VEC &,
135  const VEC &n, scalar_type r, const VEC2 &, VECR &F,
136  MAT &dlambda, VECR &dg, MAT &dn, MAT &dVs) {
137  size_type N = gmm::vect_size(lambda);
138  scalar_type nn = gmm::vect_norm2(n);
139  scalar_type lambdan = gmm::vect_sp(lambda, n)/nn;
140  scalar_type lambdan_aug = gmm::neg(lambdan + r * g);
141 
142  gmm::clear(dg); gmm::clear(dVs); gmm::clear(F);
143  gmm::clear(dn); gmm::clear(dlambda);
144  // At this stage, F = P_{B_T}, dVs = d_v P_{B_T}, dn = d_n P_{B_T}
145  // and dg = d_tau P_{B_T}.
146 
147  if (lambdan_aug > scalar_type(0)) {
148  gmm::add(gmm::scaled(n, r/nn), dg);
149  gmm::rank_one_update(dlambda, n, gmm::scaled(n, scalar_type(1)/(nn*nn)));
150  gmm::rank_one_update(dn, gmm::scaled(n, scalar_type(1)/(nn*nn)), lambda);
151  gmm::rank_one_update(dn,
152  gmm::scaled(n,(lambdan_aug-lambdan)/(nn*nn*nn)), n);
153  for (size_type j = 0; j < N; ++j) dn(j,j) -= lambdan_aug/nn;
154  }
155  gmm::add(gmm::scaled(n, -lambdan_aug/nn), F);
156 
157  gmm::scale(dVs, -r);
158  }
159 
160 
161 
162 #elif FRICTION_LAW == 3 // Dummy law for test
163 
164  template <typename VEC, typename VEC2, typename VECR>
165  void aug_friction(const VEC &lambda, scalar_type g, const VEC &Vs,
166  const VEC &n, scalar_type r, const VEC2 &f, VECR &F) {
167  gmm::copy(gmm::scaled(lambda, g*r*f[0]), F); // dummy
168  gmm::copy(gmm::scaled(Vs, g*r*f[0]), F); // dummy
169 
170  gmm::copy(n, F);
171  }
172 
173  template <typename VEC, typename VEC2, typename VECR, typename MAT>
174  void aug_friction_grad(const VEC &lambda, scalar_type g, const VEC &Vs,
175  const VEC &n, scalar_type r, const VEC2 &f, VECR &F,
176  MAT &dlambda, VECR &dg, MAT &dn, MAT &dVs) {
177  gmm::copy(gmm::scaled(lambda, g*r*f[0]), F); // dummy
178  gmm::copy(gmm::scaled(Vs, g*r*f[0]), F); // dummy
179 
180  gmm::copy(n, F);
181  gmm::clear(dlambda);
182  gmm::clear(dg);
183  gmm::clear(dVs);
184  gmm::copy(gmm::identity_matrix(), dn);
185  }
186 
187 #elif FRICTION_LAW == 4 // Dummy law for test
188 
189  template <typename VEC, typename VEC2, typename VECR>
190  void aug_friction(const VEC &lambda, scalar_type g, const VEC &Vs,
191  const VEC &n, scalar_type r, const VEC2 &f, VECR &F) {
192  gmm::copy(gmm::scaled(lambda, g*r*f[0]*n[0]*Vs[0]), F); // dummy
193  gmm::copy(lambda, F);
194  }
195 
196  template <typename VEC, typename VEC2, typename VECR, typename MAT>
197  void aug_friction_grad(const VEC &lambda, scalar_type g, const VEC &Vs,
198  const VEC &n, scalar_type r, const VEC2 &f, VECR &F,
199  MAT &dlambda, VECR &dg, MAT &dn, MAT &dVs) {
200  gmm::copy(gmm::scaled(lambda, g*r*f[0]*n[0]*Vs[0]), F); // dummy
201  gmm::clear(dn);
202  gmm::clear(dg);
203  gmm::clear(dVs);
204  gmm::copy(lambda, F);
205  gmm::copy(gmm::identity_matrix(), dlambda);
206  }
207 
208 #elif FRICTION_LAW == 5 // Dummy law for test
209 
210  template <typename VEC, typename VEC2, typename VECR>
211  void aug_friction(const VEC &lambda, scalar_type g, const VEC &Vs,
212  const VEC &n, scalar_type r, const VEC2 &f, VECR &F) {
213  gmm::copy(gmm::scaled(lambda, g*r*f[0]*n[0]*Vs[0]), F); // dummy
214  gmm::clear(F); F[0] = g;
215  }
216 
217  template <typename VEC, typename VEC2, typename VECR, typename MAT>
218  void aug_friction_grad(const VEC &lambda, scalar_type g, const VEC &Vs,
219  const VEC &n, scalar_type r, const VEC2 &f, VECR &F,
220  MAT &dlambda, VECR &dg, MAT &dn, MAT &dVs) {
221  gmm::copy(gmm::scaled(lambda, g*r*f[0]*n[0]*Vs[0]), F); // dummy
222  gmm::clear(dlambda);
223  gmm::clear(dn);
224  gmm::clear(dg);
225  gmm::clear(dVs);
226  gmm::clear(F); F[0] = g;
227  dg[0] = 1.;
228  }
229 
230 #endif
231 
232 
233  //=========================================================================
234  //
235  // Large sliding brick. Work in progress
236  //
237  //=========================================================================
238 
239  // For the moment, with raytrace detection and integral unsymmetric
240  // Alart-Curnier augmented Lagrangian
241 
242 
243  struct integral_large_sliding_contact_brick : public virtual_brick {
244 
245  multi_contact_frame &mcf;
246  bool with_friction;
247 
248 
249  virtual void asm_real_tangent_terms(const model &md, size_type /* ib */,
250  const model::varnamelist &vl,
251  const model::varnamelist &dl,
252  const model::mimlist &mims,
253  model::real_matlist &matl,
254  model::real_veclist &vecl,
255  model::real_veclist &,
256  size_type region,
257  build_version version) const;
258 
259  integral_large_sliding_contact_brick(multi_contact_frame &mcff,
260  bool with_fric)
261  : mcf(mcff), with_friction(with_fric) {
262  set_flags("Integral large sliding contact brick",
263  false /* is linear*/, false /* is symmetric */,
264  false /* is coercive */, true /* is real */,
265  false /* is complex */);
266  }
267 
268  };
269 
270 
271  struct gauss_point_precomp {
272  size_type N;
273  fem_precomp_pool fppool;
274  const multi_contact_frame &mcf;
275  const model &md;
276  const multi_contact_frame::contact_pair *cp;
277 
278  const base_node &x(void) const { return cp->slave_point; }
279  const base_node &nx(void) const { return cp->slave_n; }
280  const base_node &y(void) const { return cp->master_point; }
281  const base_node &y_ref(void) const { return cp->master_point_ref; }
282  const base_node &ny(void) const { return cp->master_n; }
283  scalar_type g(void) const { return cp->signed_dist; }
284 
285  base_matrix I_nxnx_;
286  bool I_nxnx_computed;
287  const base_matrix &I_nxnx(void) {
288  if (!I_nxnx_computed) {
289  gmm::copy(gmm::identity_matrix(), I_nxnx_);
290  gmm::rank_one_update(I_nxnx_, nx(), gmm::scaled(nx(),scalar_type(-1)));
291  I_nxnx_computed = true;
292  }
293  return I_nxnx_;
294  }
295 
296  base_matrix I_nyny_;
297  bool I_nyny_computed;
298  const base_matrix &I_nyny(void) {
299  if (!I_nyny_computed) {
300  gmm::copy(gmm::identity_matrix(), I_nyny_);
301  gmm::rank_one_update(I_nyny_, ny(), gmm::scaled(ny(),scalar_type(-1)));
302  I_nyny_computed = true;
303  }
304  return I_nyny_;
305  }
306 
307  base_matrix I_nxny_;
308  bool I_nxny_computed;
309  const base_matrix &I_nxny(void) {
310  if (!I_nxny_computed) {
311  gmm::copy(gmm::identity_matrix(), I_nxny_);
312  gmm::rank_one_update(I_nxny_, nx(),
313  gmm::scaled(ny(),scalar_type(-1)/nxny));
314  I_nxny_computed = true;
315  }
316  return I_nxny_;
317  }
318 
319  scalar_type nxny;
320  scalar_type nxdotny(void) const { return nxny; }
321 
322  bool isrigid_;
323  bool isrigid(void) { return isrigid_; }
324 
325  base_tensor base_ux, base_uy, base_lx, base_ly;
326  base_matrix vbase_ux_, vbase_uy_, vbase_lx_, vbase_ly_;
327  bool vbase_ux_init, vbase_uy_init, vbase_lx_init, vbase_ly_init;
328  base_tensor grad_base_ux, grad_base_uy, vgrad_base_ux_, vgrad_base_uy_;
329  bool vgrad_base_ux_init, vgrad_base_uy_init;
330  bool have_lx, have_ly;
331 
332  fem_interpolation_context ctx_ux_, ctx_uy_, ctx_lx_, ctx_ly_;
333  bool ctx_ux_init, ctx_uy_init, ctx_lx_init, ctx_ly_init;
334  base_matrix Gx, Gy;
335  const mesh_fem *mf_ux_, *mf_uy_, *mf_lx_, *mf_ly_;
336  gmm::sub_interval I_ux_, I_uy_, I_lx_, I_ly_;
337  pfem pf_ux, pf_uy, pf_lx, pf_ly;
338  size_type ndof_ux_, qdim_ux, ndof_uy_, qdim_uy, ndof_lx_, qdim_lx;
339  size_type ndof_ly_, qdim_ly, cvx_, cvy_, ibx, iby;
340  short_type fx, fy;
341  bgeot::pgeometric_trans pgtx, pgty;
342  const mesh_im *mim;
343  pintegration_method pim;
344  scalar_type weight_;
345 
346  scalar_type weight(void) { return weight_; }
347 
348  const mesh &meshx(void) const { return mf_ux_->linked_mesh(); }
349  const mesh &meshy(void) const { return mf_uy_->linked_mesh(); }
350  const mesh_fem *mf_ux(void) const { return mf_ux_; }
351  const mesh_fem *mf_uy(void) const { return mf_uy_; }
352  const mesh_fem *mf_lx(void) const { return mf_lx_; }
353  const mesh_fem *mf_ly(void) const { return mf_ly_; }
354  size_type ndof_ux(void) const { return ndof_ux_; }
355  size_type ndof_uy(void) const { return ndof_uy_; }
356  size_type ndof_lx(void) const { return ndof_lx_; }
357  size_type ndof_ly(void) const { return ndof_ly_; }
358  size_type cvx(void) const { return cvx_; }
359  size_type cvy(void) const { return cvy_; }
360  const gmm::sub_interval I_ux(void) const { return I_ux_; }
361  const gmm::sub_interval I_uy(void) const { return I_uy_; }
362  const gmm::sub_interval I_lx(void) const { return I_lx_; }
363  const gmm::sub_interval I_ly(void) const { return I_ly_; }
364 
365 
366  fem_interpolation_context &ctx_ux(void) {
367  if (!ctx_ux_init) {
368  bgeot::vectors_to_base_matrix(Gx, meshx().points_of_convex(cvx_));
369  pfem_precomp pfp_ux
370  = fppool(pf_ux, pim->approx_method()->pintegration_points());
371  ctx_ux_ = fem_interpolation_context(pgtx, pfp_ux, cp->slave_ind_pt,
372  Gx, cvx_, fx);
373  ctx_ux_init = true;
374  }
375  return ctx_ux_;
376  }
377 
378  fem_interpolation_context &ctx_lx(void) {
379  GMM_ASSERT1(have_lx, "No multiplier defined on the slave surface");
380  if (!ctx_lx_init) {
381  pfem_precomp pfp_lx
382  = fppool(pf_lx, pim->approx_method()->pintegration_points());
383  ctx_lx_ = fem_interpolation_context(pgtx, pfp_lx, cp->slave_ind_pt,
384  ctx_ux().G(), cvx_, fx);
385  ctx_lx_init = true;
386  }
387  return ctx_lx_;
388  }
389 
390  fem_interpolation_context &ctx_uy(void) {
391  GMM_ASSERT1(!isrigid(), "Rigid obstacle master node: no fem defined");
392  if (!ctx_uy_init) {
393  bgeot::vectors_to_base_matrix(Gy, meshy().points_of_convex(cvy_));
394  ctx_uy_ = fem_interpolation_context(pgty, pf_uy, y_ref(), Gy, cvy_, fy);
395  ctx_uy_init = true;
396  }
397  return ctx_uy_;
398  }
399 
400  fem_interpolation_context &ctx_ly(void) {
401  GMM_ASSERT1(have_ly, "No multiplier defined on the master surface");
402  if (!ctx_ly_init) {
403  ctx_ly_ = fem_interpolation_context(pgty, pf_ly, y_ref(),
404  ctx_uy().G(), cvy_, fy);
405  ctx_ly_init = true;
406  }
407  return ctx_ly_;
408  }
409 
410  const base_matrix &vbase_ux(void) {
411  if (!vbase_ux_init) {
412  ctx_ux().base_value(base_ux);
413  vectorize_base_tensor(base_ux, vbase_ux_, ndof_ux_, qdim_ux, N);
414  vbase_ux_init = true;
415  }
416  return vbase_ux_;
417  }
418 
419  const base_matrix &vbase_uy(void) {
420  if (!vbase_uy_init) {
421  ctx_uy().base_value(base_uy);
422  vectorize_base_tensor(base_uy, vbase_uy_, ndof_uy_, qdim_uy, N);
423  vbase_uy_init = true;
424  }
425  return vbase_uy_;
426  }
427 
428  const base_matrix &vbase_lx(void) {
429  if (!vbase_lx_init) {
430  ctx_lx().base_value(base_lx);
431  vectorize_base_tensor(base_lx, vbase_lx_, ndof_lx_, qdim_lx, N);
432  vbase_lx_init = true;
433  }
434  return vbase_lx_;
435  }
436 
437  const base_matrix &vbase_ly(void) {
438  if (!vbase_ly_init) {
439  ctx_ly().base_value(base_ly);
440  vectorize_base_tensor(base_ly, vbase_ly_, ndof_ly_, qdim_ly, N);
441  vbase_ly_init = true;
442  }
443  return vbase_ly_;
444  }
445 
446  const base_tensor &vgrad_base_ux(void) {
447  if (!vgrad_base_ux_init) {
448  ctx_ux().grad_base_value(grad_base_ux);
449  vectorize_grad_base_tensor(grad_base_ux, vgrad_base_ux_, ndof_ux_,
450  qdim_ux, N);
451  vgrad_base_ux_init = true;
452  }
453  return vgrad_base_ux_;
454  }
455 
456  const base_tensor &vgrad_base_uy(void) {
457  if (!vgrad_base_uy_init) {
458  ctx_uy().grad_base_value(grad_base_uy);
459  vectorize_grad_base_tensor(grad_base_uy, vgrad_base_uy_, ndof_uy_,
460  qdim_uy, N);
461  vgrad_base_uy_init = true;
462  }
463  return vgrad_base_uy_;
464  }
465 
466  base_small_vector lambda_x_, lambda_y_;
467  bool lambda_x_init, lambda_y_init;
468  base_vector coeff;
469 
470  const base_small_vector &lx(void) {
471  if (!lambda_x_init) {
472  pfem pf = ctx_lx().pf();
473  slice_vector_on_basic_dof_of_element(*mf_lx_,mcf.mult_of_boundary(ibx),
474  cvx_, coeff);
475  pf->interpolation(ctx_lx(), coeff, lambda_x_, dim_type(N));
476  lambda_x_init = true;
477  }
478  return lambda_x_;
479  }
480 
481  const base_small_vector &ly(void) {
482  if (!lambda_y_init) {
483  pfem pf = ctx_ly().pf();
484  slice_vector_on_basic_dof_of_element(*mf_ly_,mcf.mult_of_boundary(iby),
485  cvy_, coeff);
486  pf->interpolation(ctx_ly(), coeff, lambda_y_, dim_type(N));
487  lambda_y_init = true;
488  }
489  return lambda_y_;
490  }
491 
492  base_matrix grad_phix_, grad_phix_inv_, grad_phiy_, grad_phiy_inv_;
493  bool grad_phix_init, grad_phix_inv_init;
494  bool grad_phiy_init, grad_phiy_inv_init;
495 
496  const base_matrix &grad_phix(void) {
497  if (!grad_phix_init) {
498  pfem pf = ctx_ux().pf();
499  slice_vector_on_basic_dof_of_element(*mf_ux_, mcf.disp_of_boundary(ibx),
500  cvx_, coeff);
501  pf->interpolation_grad(ctx_ux(), coeff, grad_phix_, dim_type(N));
502  gmm::add(gmm::identity_matrix(), grad_phix_);
503  grad_phix_init = true;
504  }
505  return grad_phix_;
506  }
507 
508  const base_matrix &grad_phix_inv(void) {
509  if (!grad_phix_inv_init) {
510  gmm::copy(grad_phix(), grad_phix_inv_);
511  /* scalar_type J = */ gmm::lu_inverse(grad_phix_inv_);
512  // if (J <= scalar_type(0)) GMM_WARNING1("Inverted element !" << J);
513  grad_phix_inv_init = true;
514  }
515  return grad_phix_inv_;
516  }
517 
518  const base_matrix &grad_phiy(void) {
519  if (!grad_phiy_init) {
520  pfem pf = ctx_uy().pf();
521  slice_vector_on_basic_dof_of_element(*mf_uy_, mcf.disp_of_boundary(iby),
522  cvy_, coeff);
523  pf->interpolation_grad(ctx_uy(), coeff, grad_phiy_, dim_type(N));
524  gmm::add(gmm::identity_matrix(), grad_phiy_);
525  grad_phiy_init = true;
526  }
527  return grad_phiy_;
528  }
529 
530  const base_matrix &grad_phiy_inv(void) {
531  if (!grad_phiy_inv_init) {
532  gmm::copy(grad_phiy(), grad_phiy_inv_);
533  /* scalar_type J = */ gmm::lu_inverse(grad_phiy_inv_);
534  // if (J <= scalar_type(0)) GMM_WARNING1("Inverted element !" << J);
535  grad_phiy_inv_init = true;
536  }
537  return grad_phiy_inv_;
538  }
539 
540  scalar_type alpha;
541  base_small_vector x0_, y0_, nx0_, Vs_;
542  bool x0_init, y0_init, nx0_init, Vs_init;
543  base_matrix grad_phiy0_;
544  bool grad_phiy0_init;
545 
546  const base_small_vector &x0(void) {
547  if (!x0_init) {
548  const model_real_plain_vector &all_x0 = mcf.disp0_of_boundary(ibx);
549  if (all_x0.size()) {
550  pfem pf = ctx_ux().pf();
551  slice_vector_on_basic_dof_of_element(*mf_ux_, all_x0, cvx_, coeff);
552  pf->interpolation(ctx_ux(), coeff, x0_, dim_type(N));
553  } else gmm::clear(x0_);
554  gmm::add(ctx_ux().xreal(), x0_);
555  x0_init = true;
556  }
557  return x0_;
558  }
559 
560  const base_small_vector &y0(void) {
561  if (!y0_init) {
562  if (!isrigid()) {
563  const model_real_plain_vector &all_y0 = mcf.disp0_of_boundary(iby);
564  if (all_y0.size()) {
565  pfem pf = ctx_uy().pf();
566  slice_vector_on_basic_dof_of_element(*mf_uy_, all_y0, cvy_, coeff);
567  pf->interpolation(ctx_uy(), coeff, y0_, dim_type(N));
568  } else gmm::clear(y0_);
569  gmm::add(ctx_uy().xreal(), y0_);
570  } else gmm::copy(y(), y0_);
571  y0_init = true;
572  }
573  return y0_;
574  }
575 
576  const base_small_vector &nx0(void) {
577  if (!nx0_init) {
578  const model_real_plain_vector &all_x0 = mcf.disp0_of_boundary(ibx);
579  if (all_x0.size()) {
580  pfem pf = ctx_ux().pf();
581  slice_vector_on_basic_dof_of_element(*mf_ux_, all_x0, cvx_, coeff);
582  base_small_vector nx00_(N);
583  base_matrix grad_phix0_(N,N);
584  compute_normal(ctx_ux(), fx, false, coeff, nx00_, nx0_, grad_phix0_);
585  nx0_ /= gmm::vect_norm2(nx0_);
586  } else gmm::clear(nx0_);
587  nx0_init = true;
588  }
589  return nx0_;
590  }
591 
592  const base_small_vector &Vs(void) { // relative velocity
593  if (!Vs_init) {
594  if (alpha != scalar_type(0)) {
595 #ifdef CONSIDER_FRAME_INDIFFERENCE
596  if (!isrigid()) {
597  gmm::add(y0(), gmm::scaled(x0(), scalar_type(-1)), Vs_);
598  gmm::add(gmm::scaled(nx0(), -g()), Vs_);
599  } else
600 #endif
601  {
602  gmm::add(x(), gmm::scaled(y(), scalar_type(-1)), Vs_);
603  gmm::add(gmm::scaled(x0(), scalar_type(-1)), Vs_);
604  gmm::add(y0(), Vs_);
605  }
606  gmm::scale(Vs_, alpha);
607  } else gmm::clear(Vs_);
608  Vs_init = true;
609  }
610  return Vs_;
611  }
612 
613  const base_matrix &grad_phiy0(void) { // grad_phiy of previous time step
614  // To be verified ...
615  if (!grad_phiy0_init) {
616  const model_real_plain_vector &all_y0 = mcf.disp0_of_boundary(iby);
617  if (!isrigid() && all_y0.size()) {
618  pfem pf = ctx_uy().pf();
619  slice_vector_on_basic_dof_of_element(*mf_uy_, all_y0, cvy_, coeff);
620  pf->interpolation_grad(ctx_uy(), coeff, grad_phiy0_, dim_type(N));
621  gmm::add(gmm::identity_matrix(), grad_phiy0_);
622  } else gmm::copy(gmm::identity_matrix(), grad_phiy0_);
623  grad_phiy0_init = true;
624  }
625  return grad_phiy0_;
626  }
627 
628  base_small_vector un;
629 
630  void set_pair(const multi_contact_frame::contact_pair &cp_) {
631  cp = &cp_;
632  I_nxnx_computed = I_nyny_computed = I_nxny_computed = false;
633  ctx_ux_init = ctx_uy_init = ctx_lx_init = ctx_ly_init = false;
634  vbase_ux_init = vbase_uy_init = vbase_lx_init = vbase_ly_init = false;
635  vgrad_base_ux_init = vgrad_base_uy_init = false;
636  lambda_x_init = lambda_y_init = false;
637  have_lx = have_ly = false;
638  grad_phix_init = grad_phiy_init = false;
639  grad_phix_inv_init = grad_phiy_inv_init = false;
640  x0_init = y0_init = Vs_init = grad_phiy0_init = false;
641  nxny = gmm::vect_sp(nx(), ny());
642  isrigid_ = (cp->irigid_obstacle != size_type(-1));
643 
644  cvx_ = cp->slave_ind_element;
645  ibx = cp->slave_ind_boundary;
646  mf_ux_ = &(mcf.mfdisp_of_boundary(ibx));
647  pf_ux = mf_ux_->fem_of_element(cvx_);
648  qdim_ux = pf_ux->target_dim();
649  ndof_ux_ = pf_ux->nb_dof(cvx_) * N / qdim_ux;
650  fx = cp->slave_ind_face;
651  pgtx = meshx().trans_of_convex(cvx_);
652  mim = &(mcf.mim_of_boundary(ibx));
653  pim = mim->int_method_of_element(cvx_);
654  weight_ = pim->approx_method()->coeff(cp->slave_ind_pt) * ctx_ux().J();
655  gmm::mult(ctx_ux().B(), pgtx->normals()[fx], un);
656  weight_ *= gmm::vect_norm2(un);
657  const std::string &name_ux = mcf.varname_of_boundary(ibx);
658  I_ux_ = md.interval_of_variable(name_ux);
659 
660  const std::string &name_lx = mcf.multname_of_boundary(ibx);
661  have_lx = (name_lx.size() > 0);
662  if (have_lx) {
663  mf_lx_ = &(mcf.mfmult_of_boundary(ibx));
664  I_lx_ = md.interval_of_variable(name_lx);
665  pf_lx = mf_lx_->fem_of_element(cvx_);
666  qdim_lx = pf_lx->target_dim();
667  ndof_lx_ = pf_lx->nb_dof(cvx_) * N / qdim_lx;
668  }
669 
670  if (!isrigid_) {
671  cvy_ = cp->master_ind_element;
672  iby = cp->master_ind_boundary;
673  fy = cp->master_ind_face;
674  mf_uy_ = &(mcf.mfdisp_of_boundary(iby));
675  pf_uy = mf_uy_->fem_of_element(cvy_);
676  qdim_uy = pf_uy->target_dim();
677  ndof_uy_ = pf_uy->nb_dof(cvy_) * N / qdim_uy;
678  pgty = meshy().trans_of_convex(cvy_);
679 
680  const std::string &name_uy = mcf.varname_of_boundary(iby);
681  I_uy_ = md.interval_of_variable(name_uy);
682  const std::string &name_ly = mcf.multname_of_boundary(iby);
683  have_ly = (name_ly.size() > 0);
684  if (have_ly) {
685  mf_ly_ = &(mcf.mfmult_of_boundary(iby));
686  I_ly_ = md.interval_of_variable(name_ly);
687  pf_ly = mf_ly_->fem_of_element(cvy_);
688  qdim_ly = pf_ly->target_dim();
689  ndof_ly_ = pf_ly->nb_dof(cvy_) * N / qdim_ly;
690  }
691  }
692  }
693 
694  gauss_point_precomp(size_type N_, const model &md_,
695  const multi_contact_frame &mcf_, scalar_type alpha_) :
696  N(N_), mcf(mcf_), md(md_),
697  I_nxnx_(N,N), I_nyny_(N,N), I_nxny_(N,N),
698  lambda_x_(N), lambda_y_(N),
699  grad_phix_(N, N), grad_phix_inv_(N, N),
700  grad_phiy_(N, N), grad_phiy_inv_(N, N), alpha(alpha_),
701  x0_(N), y0_(N), nx0_(N), Vs_(N), grad_phiy0_(N, N), un(N) {}
702 
703  };
704 
705 /* static void do_test_F(size_type N) { */
706 
707 /* base_matrix dlambdaF(N, N), dnF(N, N), dVsF(N, N); */
708 /* base_small_vector F(N), dgF(N); */
709 
710 /* scalar_type EPS = 5E-9; */
711 /* for (size_type k = 0; k < 100; ++k) { */
712 /* base_small_vector lambda_r(N), Vs_r(N), nx_r(N), f_coeff_r(3); */
713 /* base_small_vector F2(N), F3(N); */
714 /* scalar_type g_r = gmm::random(1.), r_r = gmm::random(); */
715 /* gmm::fill_random(lambda_r); */
716 /* gmm::fill_random(Vs_r); */
717 /* gmm::fill_random(nx_r); */
718 /* gmm::scale(nx_r, 1./gmm::vect_norm2(nx_r)); */
719 /* f_coeff_r[0] = gmm::random(); */
720 /* f_coeff_r[1] = gmm::random(); */
721 /* f_coeff_r[2] = gmm::random(); */
722 
723 /* aug_friction(lambda_r, g_r, Vs_r, nx_r, r_r, f_coeff_r, F); */
724 /* aug_friction_grad(lambda_r, g_r, Vs_r, nx_r, r_r, f_coeff_r, F2, */
725 /* dlambdaF, dgF, dnF, dVsF); */
726 /* GMM_ASSERT1(gmm::vect_dist2(F2, F) < 1E-7, "bad F"); */
727 
728 /* base_small_vector dlambda(N); */
729 /* gmm::fill_random(dlambda); */
730 
731 
732 /* gmm::add(gmm::scaled(dlambda, EPS), nx_r); */
733 /* aug_friction(lambda_r, g_r, Vs_r, nx_r, r_r, f_coeff_r, F2); */
734 
735 /* gmm::mult(dnF, gmm::scaled(dlambda, EPS), F, F3); */
736 /* if (gmm::vect_dist2(F2, F3)/EPS > 1E-4) { */
737 /* cout << "lambda_r = " << lambda_r << " Vs_r = " << Vs_r */
738 /* << " nx_r = " << nx_r << endl << "g_r = " << g_r */
739 /* << " r_r = " << r_r << " f = " << f_coeff_r << endl; */
740 /* cout << "diff = " << gmm::vect_dist2(F2, F3)/EPS << endl; */
741 /* GMM_ASSERT1(false, "bad n derivative"); */
742 /* } */
743 
744 /* gmm::add(gmm::scaled(dlambda, -EPS), nx_r); */
745 
746 
747 /* gmm::add(gmm::scaled(dlambda, EPS), lambda_r); */
748 /* aug_friction(lambda_r, g_r, Vs_r, nx_r, r_r, f_coeff_r, F2); */
749 /* gmm::mult(dlambdaF, gmm::scaled(dlambda, EPS), F, F3); */
750 /* if (gmm::vect_dist2(F2, F3)/EPS > 1E-6) { */
751 /* cout << "diff = " << gmm::vect_dist2(F2, F3)/EPS << endl; */
752 /* GMM_ASSERT1(false, "bad lambda derivative"); */
753 /* } */
754 /* gmm::add(gmm::scaled(dlambda, -EPS), lambda_r); */
755 
756 
757 /* gmm::add(gmm::scaled(dlambda, EPS), Vs_r); */
758 /* aug_friction(lambda_r, g_r, Vs_r, nx_r, r_r, f_coeff_r, F2); */
759 /* gmm::mult(dVsF, gmm::scaled(dlambda, EPS), F, F3); */
760 /* if (gmm::vect_dist2(F2, F3)/EPS > 1E-6) { */
761 /* cout << "diff = " << gmm::vect_dist2(F2, F3)/EPS << endl; */
762 /* GMM_ASSERT1(false, "bad Vs derivative"); */
763 /* } */
764 /* gmm::add(gmm::scaled(dlambda, -EPS), Vs_r); */
765 
766 
767 /* g_r += EPS; */
768 /* aug_friction(lambda_r, g_r, Vs_r, nx_r, r_r, f_coeff_r, F2); */
769 /* gmm::add(gmm::scaled(dgF, EPS), F, F3); */
770 /* if (gmm::vect_dist2(F2, F3)/EPS > 1E-6) { */
771 /* cout << "diff = " << gmm::vect_dist2(F2, F3)/EPS << endl; */
772 /* GMM_ASSERT1(false, "bad g derivative"); */
773 /* } */
774 /* g_r -= EPS; */
775 /* } */
776 /* } */
777 
778 
779  void integral_large_sliding_contact_brick::asm_real_tangent_terms
780  (const model &md, size_type /* ib */, const model::varnamelist &vl,
781  const model::varnamelist &dl, const model::mimlist &/* mims */,
782  model::real_matlist &matl, model::real_veclist &vecl,
783  model::real_veclist &, size_type /* region */,
784  build_version version) const {
785 
786  // Data : r, friction_coeff.
787  GMM_ASSERT1((with_friction && dl.size() >= 2 && dl.size() <= 3)
788  || (!with_friction && dl.size() >= 1 && dl.size() <= 2),
789  "Wrong number of data for integral large sliding contact brick");
790 
791  GMM_ASSERT1(vl.size() == mcf.nb_variables() + mcf.nb_multipliers(),
792  "For the moment, it is not allowed to add boundaries to "
793  "the multi contact frame object after a model brick has "
794  "been added.");
795 
796  const model_real_plain_vector &vr = md.real_variable(dl[0]);
797  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Large sliding contact "
798  "brick: parameter r should be a scalar");
799  scalar_type r = vr[0];
800 
801  model_real_plain_vector f_coeff;
802  if (with_friction) {
803  f_coeff = md.real_variable(dl[1]);
804  GMM_ASSERT1(gmm::vect_size(f_coeff) <= 3,
805  "Large sliding contact "
806  "brick: the friction law has less than 3 parameters");
807  }
808  if (gmm::vect_size(f_coeff) == 0) // default: no friction
809  { f_coeff.resize(1); f_coeff[0] = scalar_type(0); }
810 
811  scalar_type alpha(0);
812  size_type ind = with_friction ? 2:1;
813  if (dl.size() >= ind+1) {
814  GMM_ASSERT1(md.real_variable(dl[ind]).size() == 1,
815  "Large sliding contact "
816  "brick: parameter alpha should be a scalar");
817  alpha = md.real_variable(dl[ind])[0];
818  }
819 
820  GMM_ASSERT1(matl.size() == 1,
821  "Large sliding contact brick should have only one term");
822  model_real_sparse_matrix &M = matl[0]; gmm::clear(M);
823  model_real_plain_vector &V = vecl[0]; gmm::clear(V);
824 
825  mcf.set_raytrace(true);
826  mcf.set_nodes_mode(0);
827  mcf.compute_contact_pairs();
828 
829  size_type N = mcf.dim();
830  base_matrix Melem;
831  base_matrix dlambdaF(N, N), dnF(N, N), dVsF(N, N);
832  base_small_vector F(N), dgF(N);
833 
834  scalar_type FMULT = 1.;
835 
836  // Stabilization for non-contact zones
837  for (size_type i = 0; i < mcf.nb_boundaries(); ++i)
838  if (mcf.is_self_contact() || mcf.is_slave_boundary(i)) {
839  size_type region = mcf.region_of_boundary(i);
840  const std::string &name_lx = mcf.multname_of_boundary(i);
841  GMM_ASSERT1(name_lx.size() > 0, "This brick need "
842  "multipliers defined on the multi_contact_frame object");
843  const mesh_fem &mflambda = mcf.mfmult_of_boundary(i);
844  const mesh_im &mim = mcf.mim_of_boundary(i);
845  const gmm::sub_interval &I = md.interval_of_variable(name_lx);
846 
847  if (version & model::BUILD_MATRIX) { // LXLX term
848  model_real_sparse_matrix M1(mflambda.nb_dof(), mflambda.nb_dof());
849  asm_mass_matrix(M1, mim, mflambda, region);
850  gmm::add(gmm::scaled(M1, FMULT/r), gmm::sub_matrix(M, I, I));
851  }
852 
853  if (version & model::BUILD_RHS) { // LX term
854  model_real_plain_vector V1(mflambda.nb_dof());
855  asm_source_term
856  (V1, mim, mflambda, mflambda,
857  md.real_variable(mcf.multname_of_boundary(i)), region);
858  gmm::add(gmm::scaled(V1, -FMULT/r), gmm::sub_vector(V, I));
859  }
860  }
861 
862  gauss_point_precomp gpp(N, md, mcf, alpha);
863 
864  // do_test_F(2); do_test_F(3);
865 
866  base_matrix auxNN1(N, N), auxNN2(N, N);
867  base_small_vector auxN1(N), auxN2(N);
868 
869  // Iterations on the contact pairs
870  for (size_type icp = 0; icp < mcf.nb_contact_pairs(); ++icp) {
871  const multi_contact_frame::contact_pair &cp = mcf.get_contact_pair(icp);
872  gpp.set_pair(cp);
873  const base_small_vector &nx = gpp.nx(), &ny = gpp.ny();
874  const mesh_fem *mf_ux = gpp.mf_ux(), *mf_lx = gpp.mf_lx(), *mf_uy(0);
875  size_type ndof_ux = gpp.ndof_ux(), ndof_uy(0), ndof_lx = gpp.ndof_lx();
876  size_type cvx = gpp.cvx(), cvy(0);
877  const gmm::sub_interval &I_ux = gpp.I_ux(), &I_lx = gpp.I_lx();
878  gmm::sub_interval I_uy;
879  bool isrigid = gpp.isrigid();
880  if (!isrigid) {
881  ndof_uy = gpp.ndof_uy(); I_uy = gpp.I_uy();
882  mf_uy = gpp.mf_uy(); cvy = gpp.cvy();
883  }
884  scalar_type weight = gpp.weight(), g = gpp.g();
885  const base_small_vector &lambda = gpp.lx();
886 
887  base_vector auxUX(ndof_ux), auxUY(ndof_uy);
888 
889  if (version & model::BUILD_MATRIX) {
890 
891  base_matrix auxUYN(ndof_uy, N);
892  base_matrix auxLXN1(ndof_lx, N), auxLXN2(ndof_lx, N);
893 
894  aug_friction_grad(lambda, g, gpp.Vs(), nx, r, f_coeff, F, dlambdaF,
895  dgF, dnF, dVsF);
896 
897 
898  const base_tensor &vgrad_base_ux = gpp.vgrad_base_ux();
899  base_matrix graddeltaunx(ndof_ux, N);
900  for (size_type i = 0; i < ndof_ux; ++i)
901  for (size_type j = 0; j < N; ++j)
902  for (size_type k = 0; k < N; ++k)
903  graddeltaunx(i, j) += nx[k] * vgrad_base_ux(i, k, j);
904 
905 #define CONSIDER_TERM1
906 #define CONSIDER_TERM2
907 #define CONSIDER_TERM3
908 
909 
910 #ifdef CONSIDER_TERM1
911  // Term -\delta\lambda(X) . \delta v(X)
912  gmm::resize(Melem, ndof_ux, ndof_lx); gmm::clear(Melem);
913  gmm::mult(gpp.vbase_ux(), gmm::transposed(gpp.vbase_lx()), Melem);
914  gmm::scale(Melem, -weight);
915  mat_elem_assembly(M, I_ux, I_lx, Melem, *mf_ux, cvx, *mf_lx, cvx);
916 #endif
917 
918 #ifdef CONSIDER_TERM2
919 
920  if (!isrigid) {
921  // Term \delta\lambda(X) . \delta v(Y)
922  gmm::resize(Melem, ndof_uy, ndof_lx); gmm::clear(Melem);
923  gmm::mult(gpp.vbase_uy(), gmm::transposed(gpp.vbase_lx()), Melem);
924  gmm::scale(Melem, weight);
925  mat_elem_assembly(M, I_uy, I_lx, Melem, *mf_uy, cvy, *mf_lx, cvx);
926 
927  // Term \lambda(X) . (\nabla \delta v(Y) (\nabla phi)^(-1)\delta y
928  gmm::clear(auxUYN);
929  const base_tensor &vgrad_base_uy = gpp.vgrad_base_uy();
930  for (size_type i = 0; i < ndof_uy; ++i)
931  for (size_type j = 0; j < N; ++j)
932  for (size_type k = 0; k < N; ++k)
933  auxUYN(i, j) += lambda[k] * vgrad_base_uy(i, k, j);
934  base_matrix lgraddeltavgradphiyinv(ndof_uy, N);
935  gmm::mult(auxUYN, gpp.grad_phiy_inv(), lgraddeltavgradphiyinv);
936 
937  // first sub term
938  gmm::resize(Melem, ndof_uy, ndof_uy); gmm::clear(Melem);
939  gmm::mult(lgraddeltavgradphiyinv, gpp.I_nxny(), auxUYN);
940  gmm::mult(auxUYN, gmm::transposed(gpp.vbase_uy()), Melem);
941  // Caution: re-use of auxUYN in second sub term
942  gmm::scale(Melem, -weight);
943  mat_elem_assembly(M, I_uy, I_uy, Melem, *mf_uy, cvy, *mf_uy, cvy);
944 
945  // Second sub term
946  gmm::resize(Melem, ndof_uy, ndof_ux); gmm::clear(Melem);
947  // Caution: re-use of auxUYN
948  // gmm::mult(lgraddeltavgradphiyinv, gpp.I_nxny(), auxUYN);
949  gmm::mult(auxUYN, gmm::transposed(gpp.vbase_ux()), Melem);
950 
951  // Third sub term
952  base_matrix auxUYUX(ndof_uy, ndof_ux);
953  gmm::mult(gpp.I_nxny(), gmm::transposed(gpp.grad_phix_inv()), auxNN1);
954  gmm::mult(lgraddeltavgradphiyinv, auxNN1, auxUYN);
955  gmm::mult(auxUYN, gmm::transposed(graddeltaunx), auxUYUX);
956  gmm::scale(auxUYUX, -g);
957  gmm::add(auxUYUX, Melem);
958  gmm::scale(Melem, weight);
959  mat_elem_assembly(M, I_uy, I_ux, Melem, *mf_uy, cvy, *mf_ux, cvx);
960  }
961 
962 #endif
963 
964 
965 #ifdef CONSIDER_TERM3
966 
967  // LXLX term
968  // Term (1/r)(I-dlambdaF)\delta\lambda\delta\mu
969  // the I of (I-dlambdaF) is skipped because globally added before
970  gmm::resize(Melem, ndof_lx, ndof_lx); gmm::clear(Melem);
971  gmm::copy(gmm::scaled(dlambdaF, scalar_type(-1)/r), auxNN1);
972  gmm::mult(gpp.vbase_lx(), auxNN1, auxLXN1);
973  gmm::mult(auxLXN1, gmm::transposed(gpp.vbase_lx()), Melem);
974  gmm::scale(Melem, weight*FMULT);
975  mat_elem_assembly(M, I_lx, I_lx, Melem, *mf_lx, cvx, *mf_lx, cvx);
976 
977  // LXUX term
978  // Term -(1/r)dnF\delta nx\delta\mu
979  gmm::resize(Melem, ndof_lx, ndof_ux); gmm::clear(Melem);
980  gmm::mult(gpp.vbase_lx(), dnF, auxLXN1);
981  gmm::mult(auxLXN1, gpp.I_nxnx(), auxLXN2);
982  gmm::mult(auxLXN2, gmm::transposed(gpp.grad_phix_inv()), auxLXN1);
983  gmm::mult(auxLXN1, gmm::transposed(graddeltaunx), Melem);
984  gmm::scale(Melem, scalar_type(1)/r);
985  // assembly factorized with the next term
986 
987  // Term -(1/r)dgF\delta g\delta\mu
988  base_vector deltamudgF(ndof_lx);
989  gmm::mult(gpp.vbase_lx(),
990  gmm::scaled(dgF, scalar_type(1)/(r*gpp.nxdotny())),
991  deltamudgF);
992 
993  // first sub term
994  gmm::mult(gpp.vbase_ux(), ny, auxUX);
995 
996  // second sub term
997  gmm::mult(gpp.I_nxnx(), gmm::scaled(ny, -g), auxN1);
998  gmm::mult(gpp.grad_phix_inv(), auxN1, auxN2);
999  gmm::mult_add(graddeltaunx, auxN2, auxUX); // auxUX -> \delta u(X) - g Dn_x
1000  gmm::rank_one_update(Melem, deltamudgF, auxUX);
1001  gmm::scale(Melem, weight*FMULT);
1002  mat_elem_assembly(M, I_lx, I_ux, Melem, *mf_lx, cvx, *mf_ux, cvx);
1003 
1004  if (!isrigid) {
1005  // LXUY term
1006  // third sub term
1007  gmm::resize(Melem, ndof_lx, ndof_uy); gmm::clear(Melem);
1008  gmm::mult(gpp.vbase_uy(), ny, auxUY); // auxUY -> \delta u(Y)
1009  gmm::rank_one_update(Melem, deltamudgF, auxUY);
1010  gmm::scale(Melem, -weight*FMULT);
1011  mat_elem_assembly(M, I_lx, I_uy, Melem, *mf_lx, cvx, *mf_uy, cvy);
1012  }
1013 
1014  if (alpha != scalar_type(0)) {
1015  // Term -(1/r) d_Vs F \delta Vs\delta\mu
1016 
1017  if (!isrigid) {
1018 #ifdef CONSIDER_FRAME_INDIFFERENCE
1019  base_matrix gphiy0gphiyinv(N, N);
1020  gmm::mult(gpp.grad_phiy0(), gpp.grad_phiy_inv(), gphiy0gphiyinv);
1021  gmm::mult(gphiy0gphiyinv, gpp.I_nxny(), auxNN1);
1022  gmm::rank_one_update(auxNN1, gpp.nx0(),
1023  gmm::scaled(gpp.ny(),scalar_type(1)/gpp.nxdotny()));
1024  gmm::mult(dVsF, auxNN1, auxNN2);
1025  // Caution: auxNN2 re-used in the second sub term
1026 
1027  // LXUX term
1028  // first sub term
1029  gmm::resize(Melem, ndof_lx, ndof_ux); gmm::clear(Melem);
1030  gmm::mult(gpp.vbase_lx(), gmm::transposed(auxNN2), auxLXN1);
1031  // Caution: auxLXN1 re-used in the third sub term
1032  gmm::mult(auxLXN1, gmm::transposed(gpp.vbase_ux()), Melem);
1033 
1034  // second sub term
1035  base_matrix auxLXUX(ndof_lx, ndof_ux);
1036  gmm::mult(auxNN2, gpp.I_nxnx(), auxNN1);
1037  gmm::mult(auxNN1, gmm::transposed(gpp.grad_phix_inv()), auxNN2);
1038  gmm::mult(gpp.vbase_lx(), gmm::transposed(auxNN2), auxLXN2);
1039  gmm::mult(auxLXN2, gmm::transposed(graddeltaunx), auxLXUX);
1040  gmm::scale(auxLXUX, -g);
1041  gmm::add(auxLXUX, Melem);
1042  gmm::scale(Melem, -weight*alpha*FMULT/r);
1043  mat_elem_assembly(M, I_lx, I_ux, Melem, *mf_lx, cvx, *mf_ux, cvx);
1044 
1045  // LXUY term
1046  // third sub term
1047  gmm::resize(Melem, ndof_lx, ndof_uy); gmm::clear(Melem);
1048  // Caution: auxLXN1 re-used
1049  gmm::mult(auxLXN1, gmm::transposed(gpp.vbase_uy()), Melem);
1050  gmm::scale(Melem, weight*alpha*FMULT/r);
1051  mat_elem_assembly(M, I_lx, I_uy, Melem, *mf_lx, cvx, *mf_uy, cvy);
1052 #else
1053  base_matrix I_gphiy0gphiyinv(N, N);
1054  gmm::mult(gmm::scaled(gpp.grad_phiy0(), scalar_type(-1)),
1055  gpp.grad_phiy_inv(), I_gphiy0gphiyinv);
1056  gmm::add(gmm::identity_matrix(), I_gphiy0gphiyinv);
1057 
1058  // LXUX term
1059  // first sub term
1060  gmm::resize(Melem, ndof_lx, ndof_ux); gmm::clear(Melem);
1061  gmm::mult(I_gphiy0gphiyinv, gpp.I_nxny(), auxNN1);
1062  for (size_type j = 0; j < N; ++j) auxNN1(j,j) -= scalar_type(1);
1063  gmm::mult(dVsF, auxNN1, auxNN2);
1064  gmm::mult(gpp.vbase_lx(), gmm::transposed(auxNN2), auxLXN2);
1065  // Caution: auxLXN2 re-used in the third sub term
1066  gmm::mult(auxLXN2, gmm::transposed(gpp.vbase_ux()), Melem);
1067 
1068  // second sub term
1069  base_matrix auxLXUX(ndof_lx, ndof_ux);
1070  gmm::mult(dVsF, I_gphiy0gphiyinv, auxNN1);
1071  gmm::mult(auxNN1, gpp.I_nxny(), auxNN2);
1072  gmm::mult(auxNN2, gmm::transposed(gpp.grad_phix_inv()), auxNN1);
1073  gmm::mult(gpp.vbase_lx(), gmm::transposed(auxNN1), auxLXN1);
1074  gmm::mult(auxLXN1, gmm::transposed(graddeltaunx), auxLXUX);
1075  gmm::scale(auxLXUX, -g);
1076  gmm::add(auxLXUX, Melem);
1077  gmm::scale(Melem, weight*alpha*FMULT/r);
1078  mat_elem_assembly(M, I_lx, I_ux, Melem, *mf_lx, cvx, *mf_ux, cvx);
1079 
1080  // LXUY term
1081  // third sub term
1082  gmm::resize(Melem, ndof_lx, ndof_uy); gmm::clear(Melem);
1083  // Caution: auxLXN2 re-used
1084  gmm::mult(auxLXN2, gmm::transposed(gpp.vbase_uy()), Melem);
1085  gmm::scale(Melem, -weight*alpha*FMULT/r);
1086  mat_elem_assembly(M, I_lx, I_uy, Melem, *mf_lx, cvx, *mf_uy, cvy);
1087 #endif
1088  } else {
1089  // LXUX term
1090  gmm::mult(gpp.vbase_lx(), gmm::transposed(dVsF), auxLXN1);
1091  gmm::mult(auxLXN1, gmm::transposed(gpp.vbase_ux()), Melem);
1092  gmm::scale(Melem, -weight*alpha*FMULT/r);
1093  mat_elem_assembly(M, I_lx, I_ux, Melem, *mf_lx, cvx, *mf_ux, cvx);
1094  }
1095  }
1096 #endif
1097  }
1098 
1099  if (version & model::BUILD_RHS) {
1100 
1101  if (!(version & model::BUILD_MATRIX))
1102  aug_friction(lambda, g, gpp.Vs(), nx, r, f_coeff, F);
1103 
1104 #ifdef CONSIDER_TERM1
1105 
1106  // Term lambda.\delta v(X)
1107  gmm::mult(gpp.vbase_ux(), lambda, auxUX);
1108  gmm::scale(auxUX, weight);
1109  vec_elem_assembly(V, I_ux, auxUX, *mf_ux, cvx);
1110 #endif
1111 
1112 #ifdef CONSIDER_TERM2
1113 
1114  // Term -lambda.\delta v(Y)
1115  if (!isrigid) {
1116  gmm::mult(gpp.vbase_uy(), lambda, auxUY);
1117  gmm::scale(auxUY, -weight);
1118  vec_elem_assembly(V, I_uy, auxUY, *mf_uy, cvy);
1119  }
1120 #endif
1121 
1122 #ifdef CONSIDER_TERM3
1123 
1124  // Term -(1/r)(lambda - F).\delta \mu
1125  // (1/r)(lambda).\delta \mu is skipped because globally added before
1126  base_vector auxLX(ndof_lx);
1127  gmm::mult(gpp.vbase_lx(), gmm::scaled(F, weight*FMULT/r), auxLX);
1128  vec_elem_assembly(V, I_lx, auxLX, *mf_lx, cvx);
1129 #endif
1130  }
1131 
1132  }
1133  }
1134 
1135 
1136  size_type add_integral_large_sliding_contact_brick_raytrace
1137  (model &md, multi_contact_frame &mcf,
1138  const std::string &dataname_r, const std::string &dataname_friction_coeff,
1139  const std::string &dataname_alpha) {
1140 
1141  bool with_friction = (dataname_friction_coeff.size() > 0);
1142  pbrick pbri
1143  = std::make_shared<integral_large_sliding_contact_brick>(mcf,
1144  with_friction);
1145 
1146  model::termlist tl; // A unique global unsymmetric term
1147  tl.push_back(model::term_description(true, false));
1148 
1149  model::varnamelist dl(1, dataname_r);
1150  if (with_friction) dl.push_back(dataname_friction_coeff);
1151  if (dataname_alpha.size()) dl.push_back(dataname_alpha);
1152 
1153  model::varnamelist vl;
1154 
1155  bool selfcontact = mcf.is_self_contact();
1156 
1157  dal::bit_vector uvar, mvar;
1158  for (size_type i = 0; i < mcf.nb_boundaries(); ++i) {
1159  size_type ind_u = mcf.ind_varname_of_boundary(i);
1160  if (!(uvar.is_in(ind_u))) {
1161  vl.push_back(mcf.varname(ind_u));
1162  uvar.add(ind_u);
1163  }
1164  size_type ind_lambda = mcf.ind_multname_of_boundary(i);
1165 
1166  if (selfcontact || mcf.is_slave_boundary(i))
1167  GMM_ASSERT1(ind_lambda != size_type(-1), "Large sliding contact "
1168  "brick: a multiplier should be associated to each slave "
1169  "boundary in the multi_contact_frame object.");
1170  if (ind_lambda != size_type(-1) && !(mvar.is_in(ind_lambda))) {
1171  vl.push_back(mcf.multname(ind_lambda));
1172  mvar.add(ind_u);
1173  }
1174  }
1175 
1176  return md.add_brick(pbri, vl, dl, tl, model::mimlist(), size_type(-1));
1177  }
1178 
1179 
1180 
1181  //=========================================================================
1182  //
1183  // Large sliding brick with field extension principle.
1184  // Deprecated. To be adapated to the high-level generic assembly
1185  //
1186  //=========================================================================
1187 
1188 #if 0
1189 
1190 #include <getfem/getfem_arch_config.h>
1191 #if GETFEM_HAVE_MUPARSER_MUPARSER_H
1192 #include <muParser/muParser.h>
1193 #elif GETFEM_HAVE_MUPARSER_H
1194 #include <muParser.h>
1195 #endif
1196 
1197  //=========================================================================
1198  // 1)- Structure which stores the contact boundaries and rigid obstacles
1199  //=========================================================================
1200 
1201 
1202  struct contact_frame {
1203  bool frictionless;
1204  size_type N;
1205  scalar_type friction_coef;
1206  std::vector<const model_real_plain_vector *> Us;
1207  std::vector<model_real_plain_vector> ext_Us;
1208  std::vector<const model_real_plain_vector *> lambdas;
1209  std::vector<model_real_plain_vector> ext_lambdas;
1210  struct contact_boundary {
1211  size_type region; // Boundary number
1212  const getfem::mesh_fem *mfu; // F.e.m. for the displacement.
1213  size_type ind_U; // Index of displacement.
1214  const getfem::mesh_fem *mflambda; // F.e.m. for the multiplier.
1215  size_type ind_lambda; // Index of multiplier.
1216  };
1217  std::vector<contact_boundary> contact_boundaries;
1218 
1219  gmm::dense_matrix< model_real_sparse_matrix * > UU;
1220  gmm::dense_matrix< model_real_sparse_matrix * > UL;
1221  gmm::dense_matrix< model_real_sparse_matrix * > LU;
1222  gmm::dense_matrix< model_real_sparse_matrix * > LL;
1223 
1224  std::vector< model_real_plain_vector *> Urhs;
1225  std::vector< model_real_plain_vector *> Lrhs;
1226 
1227 
1228 
1229  std::vector<std::string> coordinates;
1230  base_node pt_eval;
1231 #if GETFEM_HAVE_MUPARSER_MUPARSER_H || GETFEM_HAVE_MUPARSER_H
1232  std::vector<mu::Parser> obstacles_parsers;
1233 #endif
1234  std::vector<std::string> obstacles;
1235  std::vector<std::string> obstacles_velocities;
1236 
1237  size_type add_U(const getfem::mesh_fem &mfu,
1238  const model_real_plain_vector &U) {
1239  size_type i = 0;
1240  for (; i < Us.size(); ++i) if (Us[i] == &U) return i;
1241  Us.push_back(&U);
1242  model_real_plain_vector ext_U(mfu.nb_basic_dof()); // means that the structure has to be build each time ... to be changed. ATTENTION : la meme variable ne doit pas etre etendue dans deux vecteurs differents.
1243  mfu.extend_vector(U, ext_U);
1244  ext_Us.push_back(ext_U);
1245  return i;
1246  }
1247 
1248  size_type add_lambda(const getfem::mesh_fem &mfl,
1249  const model_real_plain_vector &l) {
1250  size_type i = 0;
1251  for (; i < lambdas.size(); ++i) if (lambdas[i] == &l) return i;
1252  lambdas.push_back(&l);
1253  model_real_plain_vector ext_l(mfl.nb_basic_dof()); // means that the structure has to be build each time ... to be changed. ATTENTION : la meme variable ne doit pas etre etendue dans deux vecteurs differents.
1254  mfl.extend_vector(l, ext_l);
1255  ext_lambdas.push_back(ext_l);
1256  return i;
1257  }
1258 
1259  void extend_vectors(void) {
1260  for (size_type i = 0; i < contact_boundaries.size(); ++i) {
1261  size_type ind_U = contact_boundaries[i].ind_U;
1262  contact_boundaries[i].mfu->extend_vector(*(Us[ind_U]), ext_Us[ind_U]);
1263  size_type ind_lambda = contact_boundaries[i].ind_lambda;
1264  contact_boundaries[i].mflambda->extend_vector(*(lambdas[ind_lambda]),
1265  ext_lambdas[ind_lambda]);
1266  }
1267  }
1268 
1269 
1270  const getfem::mesh_fem &mfu_of_boundary(size_type n) const
1271  { return *(contact_boundaries[n].mfu); }
1272  const getfem::mesh_fem &mflambda_of_boundary(size_type n) const
1273  { return *(contact_boundaries[n].mflambda); }
1274  const model_real_plain_vector &disp_of_boundary(size_type n) const
1275  { return ext_Us[contact_boundaries[n].ind_U]; }
1276  const model_real_plain_vector &lambda_of_boundary(size_type n) const
1277  { return ext_lambdas[contact_boundaries[n].ind_lambda]; }
1278  size_type region_of_boundary(size_type n) const
1279  { return contact_boundaries[n].region; }
1280  model_real_sparse_matrix &UU_matrix(size_type n, size_type m) const
1281  { return *(UU(contact_boundaries[n].ind_U, contact_boundaries[m].ind_U)); }
1282  model_real_sparse_matrix &LU_matrix(size_type n, size_type m) const {
1283  return *(LU(contact_boundaries[n].ind_lambda,
1284  contact_boundaries[m].ind_U));
1285  }
1286  model_real_sparse_matrix &UL_matrix(size_type n, size_type m) const {
1287  return *(UL(contact_boundaries[n].ind_U,
1288  contact_boundaries[m].ind_lambda));
1289  }
1290  model_real_sparse_matrix &LL_matrix(size_type n, size_type m) const {
1291  return *(LL(contact_boundaries[n].ind_lambda,
1292  contact_boundaries[m].ind_lambda));
1293  }
1294  model_real_plain_vector &U_vector(size_type n) const
1295  { return *(Urhs[contact_boundaries[n].ind_U]); }
1296  model_real_plain_vector &L_vector(size_type n) const
1297  { return *(Lrhs[contact_boundaries[n].ind_lambda]); }
1298 
1299  contact_frame(size_type NN) : N(NN), coordinates(N), pt_eval(N) {
1300  if (N > 0) coordinates[0] = "x";
1301  if (N > 1) coordinates[1] = "y";
1302  if (N > 2) coordinates[2] = "z";
1303  if (N > 3) coordinates[3] = "w";
1304  GMM_ASSERT1(N <= 4, "Complete the definition for contact in "
1305  "dimension greater than 4");
1306  }
1307 
1308  size_type add_obstacle(const std::string &obs) {
1309  size_type ind = obstacles.size();
1310  obstacles.push_back(obs);
1311  obstacles_velocities.push_back("");
1312 #if GETFEM_HAVE_MUPARSER_MUPARSER_H || GETFEM_HAVE_MUPARSER_H
1313 
1314  mu::Parser mu;
1315  obstacles_parsers.push_back(mu);
1316  obstacles_parsers[ind].SetExpr(obstacles[ind]);
1317  for (size_type k = 0; k < N; ++k)
1318  obstacles_parsers[ind].DefineVar(coordinates[k], &pt_eval[k]);
1319 #else
1320  GMM_ASSERT1(false, "You have to link muparser with getfem to deal "
1321  "with rigid body obstacles");
1322 #endif
1323  return ind;
1324  }
1325 
1326  size_type add_boundary(const getfem::mesh_fem &mfu,
1327  const model_real_plain_vector &U,
1328  const getfem::mesh_fem &mfl,
1329  const model_real_plain_vector &l,
1330  size_type reg) {
1331  contact_boundary cb;
1332  cb.region = reg;
1333  cb.mfu = &mfu;
1334  cb.mflambda = &mfl;
1335  cb.ind_U = add_U(mfu, U);
1336  cb.ind_lambda = add_lambda(mfl, l);
1337  size_type ind = contact_boundaries.size();
1338  contact_boundaries.push_back(cb);
1339  gmm::resize(UU, ind+1, ind+1);
1340  gmm::resize(UL, ind+1, ind+1);
1341  gmm::resize(LU, ind+1, ind+1);
1342  gmm::resize(LL, ind+1, ind+1);
1343  gmm::resize(Urhs, ind+1);
1344  gmm::resize(Lrhs, ind+1);
1345  return ind;
1346  }
1347 
1348  };
1349 
1350 
1351  //=========================================================================
1352  // 2)- Structure which computes the contact pairs, rhs and tangent terms
1353  //=========================================================================
1354 
1355  struct contact_elements {
1356 
1357  contact_frame &cf; // contact frame description.
1358 
1359  // list des enrichissements pour ses points : y0, d0, element ...
1360  bgeot::rtree element_boxes; // influence regions of boundary elements
1361  // list des enrichissements of boundary elements
1362  std::vector<size_type> boundary_of_elements;
1363  std::vector<size_type> ind_of_elements;
1364  std::vector<size_type> face_of_elements;
1365  std::vector<base_node> unit_normal_of_elements;
1366 
1367  contact_elements(contact_frame &ccf) : cf(ccf) {}
1368  void init(void);
1369  bool add_point_contribution(size_type boundary_num,
1370  getfem::fem_interpolation_context &ctxu,
1371  getfem::fem_interpolation_context &ctxl,
1372  scalar_type weight, scalar_type f_coeff,
1373  scalar_type r, model::build_version version);
1374  };
1375 
1376 
1377  void contact_elements::init(void) {
1378  fem_precomp_pool fppool;
1379  // compute the influence regions of boundary elements. To be run
1380  // before the assembly of contact terms.
1381  element_boxes.clear();
1382  unit_normal_of_elements.resize(0);
1383  boundary_of_elements.resize(0);
1384  ind_of_elements.resize(0);
1385  face_of_elements.resize(0);
1386 
1387  size_type N = 0;
1388  base_matrix G;
1389  model_real_plain_vector coeff;
1390  cf.extend_vectors();
1391  for (size_type i = 0; i < cf.contact_boundaries.size(); ++i) {
1392  size_type bnum = cf.region_of_boundary(i);
1393  const mesh_fem &mfu = cf.mfu_of_boundary(i);
1394  const model_real_plain_vector &U = cf.disp_of_boundary(i);
1395  const mesh &m = mfu.linked_mesh();
1396  if (i == 0) N = m.dim();
1397  GMM_ASSERT1(m.dim() == N,
1398  "Meshes are of mixed dimensions, cannot deal with that");
1399  base_node val(N), bmin(N), bmax(N), n0(N), n(N), n_mean(N);
1400  base_matrix grad(N,N);
1401  mesh_region region = m.region(bnum);
1402  GMM_ASSERT1(mfu.get_qdim() == N,
1403  "Wrong mesh_fem qdim to compute contact pairs");
1404 
1405  dal::bit_vector points_already_interpolated;
1406  std::vector<base_node> transformed_points(m.nb_max_points());
1407  for (getfem::mr_visitor v(region,m); !v.finished(); ++v) {
1408  size_type cv = v.cv();
1409  bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
1410  pfem pf_s = mfu.fem_of_element(cv);
1411  size_type nbd_t = pgt->nb_points();
1412  slice_vector_on_basic_dof_of_element(mfu, U, cv, coeff);
1413  bgeot::vectors_to_base_matrix
1414  (G, mfu.linked_mesh().points_of_convex(cv));
1415 
1416  pfem_precomp pfp = fppool(pf_s, &(pgt->geometric_nodes()));
1417  fem_interpolation_context ctx(pgt, pfp, size_type(-1), G, cv);
1418 
1419  size_type nb_pt_on_face = 0;
1420  gmm::clear(n_mean);
1421  for (short_type ip = 0; ip < nbd_t; ++ip) {
1422  size_type ind = m.ind_points_of_convex(cv)[ip];
1423 
1424  // computation of transformed vertex
1425  if (!(points_already_interpolated.is_in(ind))) {
1426  ctx.set_ii(ip);
1427  pf_s->interpolation(ctx, coeff, val, dim_type(N));
1428  val += ctx.xreal();
1429  transformed_points[ind] = val;
1430  points_already_interpolated.add(ind);
1431  } else {
1432  val = transformed_points[ind];
1433  }
1434  // computation of unit normal vector if the vertex is on the face
1435  bool is_on_face = false;
1436  bgeot::pconvex_structure cvs = pgt->structure();
1437  for (size_type k = 0; k < cvs->nb_points_of_face(v.f()); ++k)
1438  if (cvs->ind_points_of_face(v.f())[k] == ip) is_on_face = true;
1439  if (is_on_face) {
1440  ctx.set_ii(ip);
1441  n0 = bgeot::compute_normal(ctx, v.f());
1442  pf_s->interpolation_grad(ctx, coeff, grad, dim_type(N));
1443  gmm::add(gmm::identity_matrix(), grad);
1444  scalar_type J = bgeot::lu_inverse(&(*(grad.begin())), N);
1445  if (J <= scalar_type(0)) GMM_WARNING1("Inverted element ! " << J);
1446  gmm::mult(gmm::transposed(grad), n0, n);
1447  n /= gmm::vect_norm2(n);
1448  n_mean += n;
1449  ++nb_pt_on_face;
1450  }
1451 
1452  if (ip == 0) // computation of bounding box
1453  bmin = bmax = val;
1454  else {
1455  for (size_type k = 0; k < N; ++k) {
1456  bmin[k] = std::min(bmin[k], val[k]);
1457  bmax[k] = std::max(bmax[k], val[k]);
1458  }
1459  }
1460  }
1461 
1462  GMM_ASSERT1(nb_pt_on_face,
1463  "This element has not vertex on considered face !");
1464 
1465  // Computation of influence box :
1466  // offset of the bounding box relatively to its "diameter"
1467  scalar_type h = bmax[0] - bmin[0];
1468  for (size_type k = 1; k < N; ++k)
1469  h = std::max(h, bmax[k] - bmin[k]);
1470  for (size_type k = 0; k < N; ++k)
1471  { bmin[k] -= h; bmax[k] += h; }
1472 
1473  // Store the influence box and additional information.
1474  element_boxes.add_box(bmin, bmax, unit_normal_of_elements.size());
1475  n_mean /= gmm::vect_norm2(n_mean);
1476  unit_normal_of_elements.push_back(n_mean);
1477  boundary_of_elements.push_back(i);
1478  ind_of_elements.push_back(cv);
1479  face_of_elements.push_back(v.f());
1480  }
1481  }
1482  element_boxes.build_tree();
1483  }
1484 
1485 
1486 
1487  bool contact_elements::add_point_contribution
1488  (size_type boundary_num, getfem::fem_interpolation_context &ctxu,
1489  getfem::fem_interpolation_context &ctxl, scalar_type weight,
1490  scalar_type /*f_coeff*/, scalar_type r, model::build_version version) {
1491  const mesh_fem &mfu = cf.mfu_of_boundary(boundary_num);
1492  const mesh_fem &mfl = cf.mflambda_of_boundary(boundary_num);
1493  const model_real_plain_vector &U = cf.disp_of_boundary(boundary_num);
1494  const model_real_plain_vector &L = cf.lambda_of_boundary(boundary_num);
1495  size_type N = mfu.get_qdim();
1496  base_node x0 = ctxu.xreal();
1497  bool noisy = false;
1498 
1499  // ----------------------------------------------------------
1500  // Computation of the point coordinates and the unit normal
1501  // vector in real configuration
1502  // ----------------------------------------------------------
1503 
1504  base_node n0 = bgeot::compute_normal(ctxu, ctxu.face_num());
1505  scalar_type face_factor = gmm::vect_norm2(n0);
1506  size_type cv = ctxu.convex_num();
1507  base_small_vector n(N), val(N), h(N);
1508  base_matrix gradinv(N,N), grad(N,N), gradtot(N,N), G;
1509  size_type cvnbdofu = mfu.nb_basic_dof_of_element(cv);
1510  size_type cvnbdofl = mfl.nb_basic_dof_of_element(cv);
1511  base_vector coeff(cvnbdofu);
1512  slice_vector_on_basic_dof_of_element(mfu, U, cv, coeff);
1513  ctxu.pf()->interpolation(ctxu, coeff, val, dim_type(N));
1514  base_node x = x0 + val;
1515 
1516  ctxu.pf()->interpolation_grad(ctxu, coeff, gradinv, dim_type(N));
1517  gmm::add(gmm::identity_matrix(), gradinv);
1518  scalar_type J = bgeot::lu_inverse(&(*(grad_inv.begin())), N); // remplacer par une resolution...
1519  if (J <= scalar_type(0)) {
1520  GMM_WARNING1("Inverted element !");
1521 
1522  GMM_ASSERT1(!(version & model::BUILD_MATRIX), "Impossible to build "
1523  "tangent matrix for large sliding contact");
1524  if (version & model::BUILD_RHS) {
1525  base_vector Velem(cvnbdofl);
1526  for (size_type i = 0; i < cvnbdofl; ++i) Velem[i] = 1E200;
1527  vec_elem_assembly(cf.L_vector(boundary_num), Velem, mfl, cv);
1528  return false;
1529  }
1530  }
1531 
1532  gmm::mult(gmm::transposed(gradinv), n0, n);
1533  n /= gmm::vect_norm2(n);
1534 
1535  // ----------------------------------------------------------
1536  // Selection of influence boxes
1537  // ----------------------------------------------------------
1538 
1539  bgeot::rtree::pbox_set bset;
1540  element_boxes.find_boxes_at_point(x, bset);
1541 
1542  if (noisy) cout << "Number of boxes found : " << bset.size() << endl;
1543 
1544  // ----------------------------------------------------------
1545  // Eliminates some influence boxes with the mean normal
1546  // criterion : should at least eliminate the original element.
1547  // ----------------------------------------------------------
1548 
1549  bgeot::rtree::pbox_set::iterator it = bset.begin(), itnext;
1550  for (; it != bset.end(); it = itnext) {
1551  itnext = it; ++itnext;
1552  if (gmm::vect_sp(unit_normal_of_elements[(*it)->id], n)
1553  >= -scalar_type(1)/scalar_type(20)) bset.erase(it);
1554  }
1555 
1556  if (noisy)
1557  cout << "Number of boxes satisfying the unit normal criterion : "
1558  << bset.size() << endl;
1559 
1560 
1561  // ----------------------------------------------------------
1562  // For each remaining influence box, compute y0, the corres-
1563  // ponding unit normal vector and eliminate wrong auto-contact
1564  // situations with a test on |x0-y0|
1565  // ----------------------------------------------------------
1566 
1567  it = bset.begin();
1568  std::vector<base_node> y0s;
1569  std::vector<base_small_vector> n0_y0s;
1570  std::vector<scalar_type> d0s;
1571  std::vector<scalar_type> d1s;
1572  std::vector<size_type> elt_nums;
1573  std::vector<fem_interpolation_context> ctx_y0s;
1574  for (; it != bset.end(); ++it) {
1575  size_type boundary_num_y0 = boundary_of_elements[(*it)->id];
1576  size_type cv_y0 = ind_of_elements[(*it)->id];
1577  short_type face_y0 = short_type(face_of_elements[(*it)->id]);
1578  const mesh_fem &mfu_y0 = cf.mfu_of_boundary(boundary_num_y0);
1579  pfem pf_s_y0 = mfu_y0.fem_of_element(cv_y0);
1580  const model_real_plain_vector &U_y0
1581  = cf.disp_of_boundary(boundary_num_y0);
1582  const mesh &m_y0 = mfu_y0.linked_mesh();
1583  bgeot::pgeometric_trans pgt_y0 = m_y0.trans_of_convex(cv_y0);
1584  bgeot::pconvex_structure cvs_y0 = pgt_y0->structure();
1585 
1586  // Find an interior point (in order to promote the more interior
1587  // y0 in case of locally non invertible transformation.
1588  size_type ind_dep_point = 0;
1589  for (; ind_dep_point < cvs_y0->nb_points(); ++ind_dep_point) {
1590  bool is_on_face = false;
1591  for (size_type k = 0;
1592  k < cvs_y0->nb_points_of_face(face_y0); ++k)
1593  if (cvs_y0->ind_points_of_face(face_y0)[k]
1594  == ind_dep_point) is_on_face = true;
1595  if (!is_on_face) break;
1596  }
1597  GMM_ASSERT1(ind_dep_point < cvs_y0->nb_points(),
1598  "No interior point found !");
1599 
1600  base_node y0_ref = pgt_y0->convex_ref()->points()[ind_dep_point];
1601 
1602  slice_vector_on_basic_dof_of_element(mfu_y0, U_y0, cv_y0, coeff);
1603  // if (pf_s_y0->need_G())
1604  bgeot::vectors_to_base_matrix(G, m_y0.points_of_convex(cv_y0));
1605 
1606  fem_interpolation_context ctx_y0(pgt_y0, pf_s_y0, y0_ref, G, cv_y0);
1607 
1608  size_type newton_iter = 0;
1609  for(;;) { // Newton algorithm to invert geometric transformation
1610 
1611  pf_s_y0->interpolation(ctx_y0, coeff, val, dim_type(N));
1612  val += ctx_y0.xreal() - x;
1613  scalar_type init_res = gmm::vect_norm2(val);
1614 
1615  if (init_res < 1E-12) break;
1616  if (newton_iter > 100) {
1617  GMM_WARNING1("Newton has failed to invert transformation"); // il faudrait faire qlq chose d'autre ... !
1618  GMM_ASSERT1(!(version & model::BUILD_MATRIX), "Impossible to build "
1619  "tangent matrix for large sliding contact");
1620  if (version & model::BUILD_RHS) {
1621  base_vector Velem(cvnbdofl);
1622  for (size_type i = 0; i < cvnbdofl; ++i) Velem[i] = 1E200;
1623  vec_elem_assembly(cf.L_vector(boundary_num), Velem, mfl, cv);
1624  return false;
1625  }
1626  }
1627 
1628  pf_s_y0->interpolation_grad(ctx_y0, coeff, grad, dim_type(N));
1629  gmm::add(gmm::identity_matrix(), grad);
1630  gmm::mult(grad, ctx_y0.K(), gradtot);
1631 
1632  std::vector<long> ipvt(N);
1633  size_type info = gmm::lu_factor(gradtot, ipvt);
1634  GMM_ASSERT1(!info, "Singular system, pivot = " << info); // il faudrait faire qlq chose d'autre ... perturber par exemple
1635  gmm::lu_solve(gradtot, ipvt, h, val);
1636 
1637  // line search
1638  bool ok = false;
1639  scalar_type alpha;
1640  for (alpha = 1; alpha >= 1E-5; alpha/=scalar_type(2)) {
1641 
1642  ctx_y0.set_xref(y0_ref - alpha*h);
1643  pf_s_y0->interpolation(ctx_y0, coeff, val, dim_type(N));
1644  val += ctx_y0.xreal() - x;
1645 
1646  if (gmm::vect_norm2(val) < init_res) { ok = true; break; }
1647  }
1648  if (!ok)
1649  GMM_WARNING1("Line search has failed to invert transformation");
1650  y0_ref -= alpha*h;
1651  ctx_y0.set_xref(y0_ref);
1652  newton_iter++;
1653  }
1654 
1655  base_node y0 = ctx_y0.xreal();
1656  base_node n0_y0 = bgeot::compute_normal(ctx_y0, face_y0);
1657  scalar_type d0_ref = pgt_y0->convex_ref()->is_in_face(face_y0, y0_ref);
1658  scalar_type d0 = d0_ref / gmm::vect_norm2(n0_y0);
1659 
1660 
1661  scalar_type d1 = d0_ref; // approximatively a distance to the element
1662  short_type ifd = short_type(-1);
1663 
1664  for (short_type k = 0; k < pgt_y0->structure()->nb_faces(); ++k) {
1665  scalar_type dd = pgt_y0->convex_ref()->is_in_face(k, y0_ref);
1666  if (dd > scalar_type(0) && dd > gmm::abs(d1)) { d1 = dd; ifd = k; }
1667  }
1668 
1669  if (ifd != short_type(-1)) {
1670  d1 /= gmm::vect_norm2(bgeot::compute_normal(ctx_y0, ifd));
1671  if (gmm::abs(d1) < gmm::abs(d0)) d1 = d0;
1672  } else d1 = d0;
1673 
1674 // size_type iptf = m_y0.ind_points_of_face_of_convex(cv_y0, face_y0)[0];
1675 // base_node ptf = x0 - m_y0.points()[iptf];
1676 // scalar_type d2 = gmm::vect_sp(ptf, n0_y0) / gmm::vect_norm2(n0_y0);
1677 
1678  if (noisy) cout << "gmm::vect_norm2(n0_y0) = " << gmm::vect_norm2(n0_y0) << endl;
1679  // Eliminates wrong auto-contact situations
1680  if (noisy) cout << "autocontact status : x0 = " << x0 << " y0 = " << y0 << " " << gmm::vect_dist2(y0, x0) << " : " << d0*0.75 << " : " << d1*0.75 << endl;
1681  if (noisy) cout << "n = " << n << " unit_normal_of_elements[(*it)->id] = " << unit_normal_of_elements[(*it)->id] << endl;
1682 
1683  if (d0 < scalar_type(0)
1684  && ((&U_y0 == &U
1685  && (gmm::vect_dist2(y0, x0) < gmm::abs(d1)*scalar_type(3)/scalar_type(4)))
1686  || gmm::abs(d1) > 0.05)) {
1687  if (noisy) cout << "Eliminated x0 = " << x0 << " y0 = " << y0
1688  << " d0 = " << d0 << endl;
1689  continue;
1690  }
1691 
1692 // if (d0 < scalar_type(0) && &(U_y0) == &U
1693 // && gmm::vect_dist2(y0, x0) < gmm::abs(d1) * scalar_type(2)
1694 // && d2 < -ctxu.J() / scalar_type(2)) {
1695 // if (noisy) cout << "Eliminated x0 = " << x0 << " y0 = " << y0
1696 // << " d0 = " << d0 << endl;
1697 // continue;
1698 // }
1699 
1700  y0s.push_back(ctx_y0.xreal()); // useful ?
1701  elt_nums.push_back((*it)->id);
1702  d0s.push_back(d0);
1703  d1s.push_back(d1);
1704  ctx_y0s.push_back(ctx_y0);
1705  n0_y0 /= gmm::vect_norm2(n0_y0);
1706  n0_y0s.push_back(n0_y0);
1707 
1708  if (noisy) cout << "dist0 = " << d0 << " dist0 * area = "
1709  << pgt_y0->convex_ref()->is_in(y0_ref) << endl;
1710  }
1711 
1712  // ----------------------------------------------------------
1713  // Compute the distance to rigid obstacles and selects the
1714  // nearest boundary/obstacle.
1715  // ----------------------------------------------------------
1716 
1717  dim_type state = 0;
1718  scalar_type d0 = 1E100, d1 = 1E100;
1719  base_small_vector grad_obs(N);
1720 
1721  size_type ibound = size_type(-1);
1722  for (size_type k = 0; k < y0s.size(); ++k)
1723  if (d1s[k] < d1) { d0 = d0s[k]; d1 = d1s[k]; ibound = k; state = 1; }
1724 
1725 
1726  size_type irigid_obstacle = size_type(-1);
1727 #if GETFEM_HAVE_MUPARSER_MUPARSER_H || GETFEM_HAVE_MUPARSER_H
1728  gmm::copy(x, cf.pt_eval);
1729  for (size_type i = 0; i < cf.obstacles.size(); ++i) {
1730  scalar_type d0_o = scalar_type(cf.obstacles_parsers[i].Eval());
1731  if (d0_o < d0) { d0 = d0_o; irigid_obstacle = i; state = 2; }
1732  }
1733  if (state == 2) {
1734  scalar_type EPS = face_factor * 1E-9;
1735  for (size_type k = 0; k < N; ++k) {
1736  cf.pt_eval[k] += EPS;
1737  grad_obs[k] =
1738  (scalar_type(cf.obstacles_parsers[irigid_obstacle].Eval())-d0)/EPS;
1739  cf.pt_eval[k] -= EPS;
1740  }
1741  }
1742 
1743 #else
1744  if (cf.obstacles.size() > 0)
1745  GMM_WARNING1("Rigid obstacles are ignored. Recompile with "
1746  "muParser to account for rigid obstacles");
1747 #endif
1748 
1749 
1750  // ----------------------------------------------------------
1751  // Print the found contact state ...
1752  // ----------------------------------------------------------
1753 
1754 
1755  if (noisy && state == 1) {
1756  cout << "Point : " << x0 << " of boundary " << boundary_num
1757  << " and element " << cv << " state = " << int(state);
1758  if (version & model::BUILD_RHS) cout << " RHS";
1759  if (version & model::BUILD_MATRIX) cout << " MATRIX";
1760  }
1761  if (state == 1) {
1762  size_type boundary_num_y0 = boundary_of_elements[elt_nums[ibound]];
1763  const mesh_fem &mfu_y0 = cf.mfu_of_boundary(boundary_num_y0);
1764  const mesh &m_y0 = mfu_y0.linked_mesh();
1765  size_type cv_y0 = ind_of_elements[elt_nums[ibound]];
1766 
1767  if (noisy) cout << " y0 = " << y0s[ibound] << " of element "
1768  << cv_y0 << " of boundary " << boundary_num_y0 << endl;
1769  for (size_type k = 0; k < m_y0.nb_points_of_convex(cv_y0); ++k)
1770  if (noisy) cout << "point " << k << " : "
1771  << m_y0.points()[m_y0.ind_points_of_convex(cv_y0)[k]] << endl;
1772  if (boundary_num_y0 == 0 && boundary_num == 0 && d0 < 0.0 && (version & model::BUILD_MATRIX)) GMM_ASSERT1(false, "oups");
1773  }
1774  if (noisy) cout << " d0 = " << d0 << endl;
1775 
1776  // ----------------------------------------------------------
1777  // Add the contributions to the tangent matrices and rhs
1778  // ----------------------------------------------------------
1779 
1780  GMM_ASSERT1(ctxu.pf()->target_dim() == 1 && ctxl.pf()->target_dim() == 1,
1781  "Large sliding contact assembly procedure has to be adapted "
1782  "to intrinsic vectorial elements. To be done.");
1783 
1784  // Eviter les calculs inutiles dans le cas state == 2 ... a voir a la fin
1785  // regarder aussi si on peut factoriser des mat_elem_assembly ...
1786 
1787  base_matrix Melem;
1788  base_vector Velem;
1789 
1790  base_tensor tl, tu;
1791  ctxl.base_value(tl);
1792  ctxu.base_value(tu);
1793 
1794  base_small_vector lambda(N);
1795  slice_vector_on_basic_dof_of_element(mfl, L, cv, coeff);
1796  ctxl.pf()->interpolation(ctxl, coeff, lambda, dim_type(N));
1797  GMM_ASSERT1(!(isnan(lambda[0])), "internal error");
1798 
1799  // Unstabilized frictionless case for the moment
1800 
1801  // auxiliary variables
1802  scalar_type aux1, aux2;
1803 
1804  if (state) {
1805 
1806  // zeta = lambda + d0 * r * n
1807  base_small_vector zeta(N);
1808  gmm::add(lambda, gmm::scaled(n, r*d0), zeta);
1809 
1810  base_tensor tgradu;
1811  ctxu.grad_base_value(tgradu);
1812 
1813  // variables for y0
1814  base_tensor tu_y0;
1815  size_type boundary_num_y0 = 0, cv_y0 = 0, cvnbdofu_y0 = 0;
1816  if (state == 1) {
1817  ctx_y0s[ibound].base_value(tu_y0);
1818  boundary_num_y0 = boundary_of_elements[elt_nums[ibound]];
1819  cv_y0 = ind_of_elements[elt_nums[ibound]];
1820  cvnbdofu_y0 = cf.mfu_of_boundary(boundary_num_y0).nb_basic_dof_of_element(cv_y0);
1821  }
1822  const mesh_fem &mfu_y0 = (state == 1) ?
1823  cf.mfu_of_boundary(boundary_num_y0) : mfu;
1824 
1825  if (version & model::BUILD_RHS) {
1826  // Rhs term Lx
1827  gmm::resize(Velem, cvnbdofl); gmm::clear(Velem);
1828 
1829  // Rhs term Lx: (1/r)\int (\lambda - P(\zeta)).\mu
1830  base_small_vector vecaux(N);
1831  gmm::copy(zeta, vecaux);
1832  De_Saxce_projection(vecaux, n, scalar_type(0));
1833  gmm::scale(vecaux, -scalar_type(1));
1834  gmm::add(lambda, vecaux);
1835  for (size_type i = 0; i < cvnbdofl; ++i)
1836  Velem[i] = tl[i/N] * vecaux[i%N] * weight/r;
1837  vec_elem_assembly(cf.L_vector(boundary_num), Velem, mfl, cv);
1838 
1839  // Rhs terms Ux, Uy: \int \lambda.(\psi(x_0) - \psi(y_0))
1840  gmm::resize(Velem, cvnbdofu); gmm::clear(Velem);
1841  for (size_type i = 0; i < cvnbdofu; ++i)
1842  Velem[i] = tu[i/N] * lambda[i%N] * weight;
1843  vec_elem_assembly(cf.U_vector(boundary_num), Velem, mfu, cv);
1844 
1845  if (state == 1) {
1846  gmm::resize(Velem, cvnbdofu_y0); gmm::clear(Velem);
1847  for (size_type i = 0; i < cvnbdofu_y0; ++i)
1848  Velem[i] = -tu_y0[i/N] * lambda[i%N] * weight;
1849  vec_elem_assembly(cf.U_vector(boundary_num_y0), Velem, mfu_y0, cv_y0);
1850  }
1851  }
1852 
1853  if (version & model::BUILD_MATRIX) {
1854 
1855  base_small_vector gradinv_n(N);
1856  gmm::mult(gradinv, n, gradinv_n);
1857 
1858  // de Saxce projection gradient and normal gradient at zeta
1859  base_matrix pgrad(N,N), pgradn(N,N);
1860  De_Saxce_projection_grad(zeta, n, scalar_type(0), pgrad);
1861  De_Saxce_projection_gradn(zeta, n, scalar_type(0), pgradn);
1862 
1863  base_small_vector pgrad_n(N), pgradn_n(N);
1864  gmm::mult(pgrad, n, pgrad_n);
1865  gmm::mult(pgradn, n, pgradn_n);
1866  base_matrix gradinv_pgrad(N,N), gradinv_pgradn(N,N);
1867  gmm::mult(gradinv, gmm::transposed(pgrad), gradinv_pgrad);
1868  gmm::mult(gradinv, gmm::transposed(pgradn), gradinv_pgradn);
1869 
1870  // Tangent term LxLx
1871  gmm::resize(Melem, cvnbdofl, cvnbdofl); gmm::clear(Melem);
1872  // -(1/r) \int \delta\lambda.\mu
1873  for (size_type i = 0; i < cvnbdofl; i += N) {
1874  aux1 = -tl[i/N] * weight/r;
1875  for (size_type j = 0; j < cvnbdofl; j += N) {
1876  aux2 = aux1 * tl[j/N];
1877  for (size_type k = 0; k < N; k++) Melem(i+k,j+k) = aux2;
1878  } // Melem(i+k,j+k) = -tl[i/N] * tl[j/N] * weight/r;
1879  }
1880  // (1/r) \int \nabla P(\zeta) (d\zeta/d\lambda)(\delta\lambda) . \mu
1881  for (size_type i = 0, ii = 0; i < cvnbdofl; ++i, ii = i%N)
1882  for (size_type j = 0, jj = 0; j < cvnbdofl; ++j, jj = j%N)
1883  Melem(i,j) += tl[i/N] * tl[j/N] * pgrad(ii,jj) * weight/r;
1884  mat_elem_assembly(cf.LL_matrix(boundary_num, boundary_num),
1885  Melem, mfl, cv, mfl, cv);
1886 
1887  // Tangent term UxLx
1888  gmm::resize(Melem, cvnbdofu, cvnbdofl); gmm::clear(Melem);
1889  // \int -\delta\lambda.\psi(x_0)
1890  for (size_type i = 0; i < cvnbdofu; i += N) {
1891  aux1 = -tu[i/N] * weight;
1892  for (size_type j = 0; j < cvnbdofl; j += N) {
1893  aux2 = aux1 * tl[j/N];
1894  for (size_type k = 0; k < N; k++) Melem(i+k,j+k) = aux2;
1895  }
1896  }
1897  mat_elem_assembly(cf.UL_matrix(boundary_num, boundary_num),
1898  Melem, mfu, cv, mfl, cv);
1899 
1900  // Tangent term LxUx
1901  if (0) { // DISABLED
1902  gmm::resize(Melem, cvnbdofl, cvnbdofu); gmm::clear(Melem);
1903  // \int d_0(\nabla P(\zeta))(dn/du)(\delta u).\mu
1904  for (size_type i = 0, ii = 0; i < cvnbdofl; ++i, ii = i%N)
1905  for (size_type j = 0, jj = 0; j < cvnbdofu; ++j, jj = j%N) {
1906  aux1 = aux2 = scalar_type(0);
1907  for (size_type k = 0; k < N; ++k) {
1908  aux1 += tgradu[j/N+N*k] * gradinv_n[k];
1909  aux2 += tgradu[j/N+N*k] * gradinv_pgrad(k,ii);
1910  }
1911  Melem(i,j) = d0 * tl[i/N] * (pgrad_n[ii] * aux1 - aux2) * n[jj] * weight;
1912  }
1913 
1914  // (1/r)\int \nabla_n P(zeta) (dn/du)(\delta u) . \mu
1915  // On peut certainement factoriser d'avantage ce terme avec le
1916  // precedent. Attendre la version avec frottement.
1917  for (size_type i = 0, ii = 0; i < cvnbdofl; ++i, ii = i%N)
1918  for (size_type j = 0, jj = 0; j < cvnbdofu; ++j, jj = j%N) {
1919  aux1 = aux2 = scalar_type(0);
1920  for (size_type k = 0; k < N; ++k) {
1921  aux1 += tgradu[j/N+N*k] * gradinv_n[k];
1922  aux2 += tgradu[j/N+N*k] * gradinv_pgradn(k,ii);
1923  }
1924  Melem(i,j) += tl[i/N] * (pgradn_n[ii] * aux1 - aux2) * n[jj] * weight / r;
1925  }
1926  mat_elem_assembly(cf.LU_matrix(boundary_num, boundary_num),
1927  Melem, mfl, cv, mfu, cv);
1928  } // DISABLED
1929 
1930  if (state == 1) {
1931 
1932  base_tensor tgradu_y0;
1933  ctx_y0s[ibound].grad_base_value(tgradu_y0);
1934 
1935  base_matrix gradinv_y0(N,N);
1936  base_small_vector ntilde_y0(N);
1937  { // calculate gradinv_y0 and ntilde_y0
1938  base_matrix grad_y0(N,N);
1939  base_vector coeff_y0(cvnbdofu_y0);
1940  const model_real_plain_vector &U_y0
1941  = cf.disp_of_boundary(boundary_num_y0);
1942  slice_vector_on_basic_dof_of_element(mfu_y0, U_y0, cv_y0, coeff_y0);
1943  ctx_y0s[ibound].pf()->interpolation_grad(ctx_y0s[ibound], coeff_y0,
1944  grad_y0, dim_type(N));
1945  gmm::add(gmm::identity_matrix(), grad_y0);
1946 
1947  gmm::copy(grad_y0, gradinv_y0);
1948  gmm::lu_inverse(gradinv_y0); // a proteger contre la non-inversibilite
1949  gmm::mult(gmm::transposed(gradinv_y0), n0_y0s[ibound], ntilde_y0); // (not unit) normal vector
1950  }
1951 
1952  // Tangent term UyLx: \int \delta\lambda.\psi(y_0)
1953  gmm::resize(Melem, cvnbdofu_y0, cvnbdofl); gmm::clear(Melem);
1954  for (size_type i = 0; i < cvnbdofu_y0; i += N) {
1955  aux1 = tu_y0[i/N] * weight;
1956  for (size_type j = 0; j < cvnbdofl; j += N) {
1957  aux2 = aux1 * tl[j/N];
1958  for (size_type k = 0; k < N; k++) Melem(i+k,j+k) = aux2;
1959  }
1960  }
1961  mat_elem_assembly(cf.UL_matrix(boundary_num_y0, boundary_num),
1962  Melem, mfu_y0, cv_y0, mfl, cv);
1963 
1964  // Tangent terms UyUx, UyUy
1965  // \int \lambda.((\nabla \psi(y_0))(I+\nabla u(y_0))^{-1}(\delta u(x_0) - \delta u(y_0)))
1966 
1967  // Tangent term UyUx
1968  gmm::resize(Melem, cvnbdofu_y0, cvnbdofu); gmm::clear(Melem);
1969  // \int \lambda.((\nabla \psi(y_0))(I+\nabla u(y_0))^{-1}\delta u(x_0))
1970  for (size_type i = 0, ii = 0; i < cvnbdofu_y0; ++i, ii = i%N)
1971  for (size_type j = 0, jj = 0; j < cvnbdofu; ++j, jj = j%N) {
1972  aux1 = scalar_type(0);
1973  for (size_type k = 0; k < N; ++k)
1974  aux1 += tgradu_y0[i/N+N*k]* gradinv_y0(k,jj);
1975  Melem(i,j) = lambda[ii] * aux1 * tu[j/N] * weight;
1976  }
1977  mat_elem_assembly(cf.UU_matrix(boundary_num_y0, boundary_num),
1978  Melem, mfu_y0, cv_y0, mfu, cv);
1979 
1980  // Tangent term UyUy
1981  gmm::resize(Melem, cvnbdofu_y0, cvnbdofu_y0); gmm::clear(Melem);
1982  // -\int \lambda.((\nabla \psi(y_0))(I+\nabla u(y_0))^{-1}\delta u(y_0))
1983  for (size_type i = 0, ii = 0; i < cvnbdofu_y0; ++i, ii = i%N)
1984  for (size_type j = 0, jj = 0; j < cvnbdofu_y0; ++j, jj = j%N) {
1985  aux1 = scalar_type(0);
1986  for (size_type k = 0; k < N; ++k)
1987  aux1 += tgradu_y0[i/N+N*k] * gradinv_y0(k,jj);
1988  Melem(i,j) = - lambda[ii] * aux1 * tu_y0[j/N] * weight;
1989  }
1990  mat_elem_assembly(cf.UU_matrix(boundary_num_y0, boundary_num_y0),
1991  Melem, mfu_y0, cv_y0, mfu_y0, cv_y0);
1992 
1993  // Tangent term LxUy
1994  gmm::resize(Melem, cvnbdofl, cvnbdofu_y0); gmm::clear(Melem);
1995  // -\int (I+\nabla u(y_0))^{-T}\nabla \delta(y_0).\delta u(y_0)(\nabla P(\zeta) n . \mu)
1996  for (size_type i = 0; i < cvnbdofl; ++i) {
1997  aux1 = tl[i/N] * pgrad_n[i%N] * weight;
1998  for (size_type j = 0; j < cvnbdofu_y0; ++j)
1999  Melem(i,j) = - aux1 * tu_y0[j/N] * ntilde_y0[j%N];
2000  }
2001  mat_elem_assembly(cf.LU_matrix(boundary_num, boundary_num_y0),
2002  Melem, mfl, cv, mfu_y0, cv_y0);
2003 
2004  // Addition to tangent term LxUx
2005  gmm::resize(Melem, cvnbdofl, cvnbdofu); gmm::clear(Melem);
2006  // \int (I+\nabla u(y_0))^{-T}\nabla \delta(y_0).\delta u(x_0)(\nabla P(\zeta) n . \mu)
2007  for (size_type i = 0; i < cvnbdofl; ++i) {
2008  aux1 = tl[i/N] * pgrad_n[i%N] * weight;
2009  for (size_type j = 0; j < cvnbdofu; ++j)
2010  Melem(i,j) = aux1 * tu[j/N] * ntilde_y0[j%N];
2011  }
2012  }
2013  else {
2014  // Addition to tangent term LxUx
2015  gmm::resize(Melem, cvnbdofl, cvnbdofu); gmm::clear(Melem);
2016  // \int (I+\nabla u(y_0))^{-T}\nabla \delta(y_0).\delta u(x_0)(\nabla P(\zeta) n . \mu)
2017  for (size_type i = 0; i < cvnbdofl; ++i) {
2018  aux1 = tl[i/N] * pgrad_n[i%N] * weight;
2019  for (size_type j = 0; j < cvnbdofu; ++j)
2020  Melem(i,j) = aux1 * tu[j/N] * grad_obs[j%N];
2021  }
2022  }
2023  mat_elem_assembly(cf.LU_matrix(boundary_num, boundary_num),
2024  Melem, mfl, cv, mfu, cv);
2025 
2026  }
2027 
2028  } else { // state == 0
2029 
2030  // Rhs term Lx: (1/r)\int \lambda.\mu
2031  if (version & model::BUILD_RHS) {
2032  gmm::resize(Velem, cvnbdofl); gmm::clear(Velem);
2033  for (size_type i = 0; i < cvnbdofl; ++i)
2034  Velem[i] = tl[i/N] * lambda[i%N] * weight/r;
2035  vec_elem_assembly(cf.L_vector(boundary_num), Velem, mfl, cv);
2036  }
2037 
2038  // Tangent term LxLx: -(1/r)\int \delta\lambda.\mu
2039  if (version & model::BUILD_MATRIX) {
2040  gmm::resize(Melem, cvnbdofl, cvnbdofl); gmm::clear(Melem);
2041  for (size_type i = 0; i < cvnbdofl; i += N) {
2042  aux1 = -tl[i/N] * weight/r;
2043  for (size_type j = 0; j < cvnbdofl; j += N) {
2044  aux2 = aux1 * tl[j/N];
2045  for (size_type k = 0; k < N; k++) Melem(i+k,j+k) = aux2;
2046  } // Melem(i+k,j+k) = -tl[i/N] * tl[j/N] * weight/r;
2047  }
2048  mat_elem_assembly(cf.LL_matrix(boundary_num, boundary_num),
2049  Melem, mfl, cv, mfl, cv);
2050  }
2051  }
2052 
2053  return true;
2054  }
2055 
2056  //=========================================================================
2057  // 3)- Large sliding contact brick
2058  //=========================================================================
2059 
2060  struct integral_large_sliding_contact_brick_field_extension : public virtual_brick {
2061 
2062 
2063  struct contact_boundary {
2064  size_type region;
2065  std::string varname;
2066  std::string multname;
2067  const mesh_im *mim;
2068  };
2069 
2070  std::vector<contact_boundary> boundaries;
2071  std::vector<std::string> obstacles;
2072 
2073  void add_boundary(const std::string &varn, const std::string &multn,
2074  const mesh_im &mim, size_type region) {
2075  contact_boundary cb;
2076  cb.region = region; cb.varname = varn; cb.multname = multn; cb.mim=&mim;
2077  boundaries.push_back(cb);
2078  }
2079 
2080  void add_obstacle(const std::string &obs)
2081  { obstacles.push_back(obs); }
2082 
2083  void build_contact_frame(const model &md, contact_frame &cf) const {
2084  for (size_type i = 0; i < boundaries.size(); ++i) {
2085  const contact_boundary &cb = boundaries[i];
2086  cf.add_boundary(md.mesh_fem_of_variable(cb.varname),
2087  md.real_variable(cb.varname),
2088  md.mesh_fem_of_variable(cb.multname),
2089  md.real_variable(cb.multname), cb.region);
2090  }
2091  for (size_type i = 0; i < obstacles.size(); ++i)
2092  cf.add_obstacle(obstacles[i]);
2093  }
2094 
2095 
2096  virtual void asm_real_tangent_terms(const model &md, size_type /* ib */,
2097  const model::varnamelist &vl,
2098  const model::varnamelist &dl,
2099  const model::mimlist &mims,
2100  model::real_matlist &matl,
2101  model::real_veclist &vecl,
2102  model::real_veclist &,
2103  size_type region,
2104  build_version version) const;
2105 
2106  integral_large_sliding_contact_brick_field_extension() {
2107  set_flags("Integral large sliding contact brick",
2108  false /* is linear*/, false /* is symmetric */,
2109  false /* is coercive */, true /* is real */,
2110  false /* is complex */);
2111  }
2112 
2113  };
2114 
2115 
2116 
2117 
2118  void integral_large_sliding_contact_brick_field_extension::asm_real_tangent_terms
2119  (const model &md, size_type /* ib */, const model::varnamelist &vl,
2120  const model::varnamelist &dl, const model::mimlist &/* mims */,
2121  model::real_matlist &matl, model::real_veclist &vecl,
2122  model::real_veclist &, size_type /* region */,
2123  build_version version) const {
2124 
2125  fem_precomp_pool fppool;
2126  base_matrix G;
2127  size_type N = md.mesh_fem_of_variable(vl[0]).linked_mesh().dim();
2128  contact_frame cf(N);
2129  build_contact_frame(md, cf);
2130 
2131  size_type Nvar = vl.size(), Nu = cf.Urhs.size(), Nl = cf.Lrhs.size();
2132  GMM_ASSERT1(Nvar == Nu+Nl, "Wrong size of variable list for integral "
2133  "large sliding contact brick");
2134  GMM_ASSERT1(matl.size() == Nvar*Nvar, "Wrong size of terms for "
2135  "integral large sliding contact brick");
2136 
2137  if (version & model::BUILD_MATRIX) {
2138  for (size_type i = 0; i < Nvar; ++i)
2139  for (size_type j = 0; j < Nvar; ++j) {
2140  gmm::clear(matl[i*Nvar+j]);
2141  if (i < Nu && j < Nu) cf.UU(i,j) = &(matl[i*Nvar+j]);
2142  if (i >= Nu && j < Nu) cf.LU(i-Nu,j) = &(matl[i*Nvar+j]);
2143  if (i < Nu && j >= Nu) cf.UL(i,j-Nu) = &(matl[i*Nvar+j]);
2144  if (i >= Nu && j >= Nu) cf.LL(i-Nu,j-Nu) = &(matl[i*Nvar+j]);
2145  }
2146  }
2147  if (version & model::BUILD_RHS) {
2148  for (size_type i = 0; i < vl.size(); ++i) {
2149  if (i < Nu) cf.Urhs[i] = &(vecl[i*Nvar]);
2150  else cf.Lrhs[i-Nu] = &(vecl[i*Nvar]);
2151  }
2152  }
2153 
2154  // Data : r, [friction_coeff,]
2155  GMM_ASSERT1(dl.size() == 2, "Wrong number of data for integral large "
2156  "sliding contact brick");
2157 
2158  const model_real_plain_vector &vr = md.real_variable(dl[0]);
2159  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Parameter r should be a scalar");
2160 
2161  const model_real_plain_vector &f_coeff = md.real_variable(dl[1]);
2162  GMM_ASSERT1(gmm::vect_size(f_coeff) == 1,
2163  "Friction coefficient should be a scalar");
2164 
2165  contact_elements ce(cf);
2166  ce.init();
2167 
2168  for (size_type bnum = 0; bnum < boundaries.size(); ++bnum) {
2169  mesh_region rg(boundaries[bnum].region);
2170  const mesh_fem &mfu=md.mesh_fem_of_variable(boundaries[bnum].varname);
2171  const mesh_fem &mfl=md.mesh_fem_of_variable(boundaries[bnum].multname);
2172  const mesh_im &mim = *(boundaries[bnum].mim);
2173  const mesh &m = mfu.linked_mesh();
2174  mfu.linked_mesh().intersect_with_mpi_region(rg);
2175 
2176  for (getfem::mr_visitor v(rg, m); !v.finished(); ++v) {
2177  // cout << "boundary " << bnum << " element " << v.cv() << endl;
2178  size_type cv = v.cv();
2179  bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
2180  pfem pf_s = mfu.fem_of_element(cv);
2181  pfem pf_sl = mfl.fem_of_element(cv);
2182  pintegration_method pim = mim.int_method_of_element(cv);
2183  bgeot::vectors_to_base_matrix(G, m.points_of_convex(cv));
2184 
2185  pfem_precomp pfpu
2186  = fppool(pf_s, pim->approx_method()->pintegration_points());
2187  pfem_precomp pfpl
2188  = fppool(pf_sl, pim->approx_method()->pintegration_points());
2189  fem_interpolation_context ctxu(pgt, pfpu, size_type(-1), G, cv, v.f());
2190  fem_interpolation_context ctxl(pgt, pfpl, size_type(-1), G, cv, v.f());
2191 
2192  for (size_type k = 0;
2193  k < pim->approx_method()->nb_points_on_face(v.f()); ++k) {
2194  size_type ind
2195  = pim->approx_method()->ind_first_point_on_face(v.f()) + k;
2196  ctxu.set_ii(ind);
2197  ctxl.set_ii(ind);
2198  if (!(ce.add_point_contribution
2199  (bnum, ctxu, ctxl,pim->approx_method()->coeff(ind),
2200  f_coeff[0], vr[0], version))) return;
2201  }
2202  }
2203  }
2204  }
2205 
2206 
2207  // r ne peut pas etre variable pour le moment.
2208  // dataname_friction_coeff ne peut pas etre variable non plus ...
2209 
2210  size_type add_integral_large_sliding_contact_brick_field_extension
2211  (model &md, const mesh_im &mim, const std::string &varname_u,
2212  const std::string &multname, const std::string &dataname_r,
2213  const std::string &dataname_friction_coeff, size_type region) {
2214 
2215  auto pbr
2216  =std::make_shared<integral_large_sliding_contact_brick_field_extension>();
2217 
2218  pbr->add_boundary(varname_u, multname, mim, region);
2219 
2220  model::termlist tl;
2221  tl.push_back(model::term_description(varname_u, varname_u, false));
2222  tl.push_back(model::term_description(varname_u, multname, false));
2223  tl.push_back(model::term_description(multname, varname_u, false));
2224  tl.push_back(model::term_description(multname, multname, false));
2225 
2226  model::varnamelist dl(1, dataname_r);
2227  dl.push_back(dataname_friction_coeff);
2228 
2229  model::varnamelist vl(1, varname_u);
2230  vl.push_back(multname);
2231 
2232  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
2233  }
2234 
2235 
2236  void add_boundary_to_large_sliding_contact_brick
2237  (model &md, size_type indbrick, const mesh_im &mim,
2238  const std::string &varname_u, const std::string &multname,
2239  size_type region) {
2240  dim_type N = md.mesh_fem_of_variable(varname_u).linked_mesh().dim();
2241  pbrick pbr = md.brick_pointer(indbrick);
2242  md.touch_brick(indbrick);
2243  integral_large_sliding_contact_brick_field_extension *p
2244  = dynamic_cast<integral_large_sliding_contact_brick_field_extension *>
2245  (const_cast<virtual_brick *>(pbr.get()));
2246  GMM_ASSERT1(p, "Wrong type of brick");
2247  p->add_boundary(varname_u, multname, mim, region);
2248  md.add_mim_to_brick(indbrick, mim);
2249 
2250  contact_frame cf(N);
2251  p->build_contact_frame(md, cf);
2252 
2253  model::varnamelist vl;
2254  size_type nvaru = 0;
2255  for (size_type i = 0; i < cf.contact_boundaries.size(); ++i)
2256  if (cf.contact_boundaries[i].ind_U >= nvaru)
2257  { vl.push_back(p->boundaries[i].varname); ++nvaru; }
2258 
2259  size_type nvarl = 0;
2260  for (size_type i = 0; i < cf.contact_boundaries.size(); ++i)
2261  if (cf.contact_boundaries[i].ind_lambda >= nvarl)
2262  { vl.push_back(p->boundaries[i].multname); ++nvarl; }
2263  md.change_variables_of_brick(indbrick, vl);
2264 
2265  model::termlist tl;
2266  for (size_type i = 0; i < vl.size(); ++i)
2267  for (size_type j = 0; j < vl.size(); ++j)
2268  tl.push_back(model::term_description(vl[i], vl[j], false));
2269 
2270  md.change_terms_of_brick(indbrick, tl);
2271  }
2272 
2273  void add_rigid_obstacle_to_large_sliding_contact_brick
2274  (model &md, size_type indbrick, const std::string &obs) { // The velocity field should be added to an (optional) parameter ... (and optionally represented by a rigid motion only ... the velocity should be modifiable ...
2275  pbrick pbr = md.brick_pointer(indbrick);
2276  md.touch_brick(indbrick);
2277  integral_large_sliding_contact_brick_field_extension *p
2278  = dynamic_cast<integral_large_sliding_contact_brick_field_extension *>
2279  (const_cast<virtual_brick *>(pbr.get()));
2280  GMM_ASSERT1(p, "Wrong type of brick");
2281  p->add_obstacle(obs);
2282  }
2283 
2284 #endif
2285 
2286  // ----------------------------------------------------------------------
2287  //
2288  // Brick for large sliding contact with friction using raytracing contact
2289  // detection and the high-level generic assemly
2290  //
2291  // ----------------------------------------------------------------------
2292 
2293  struct intergral_large_sliding_contact_brick_raytracing
2294  : public virtual_brick {
2295 
2296  struct contact_boundary {
2297  size_type region;
2298  std::string varname_u;
2299  std::string varname_lambda;
2300  std::string varname_w;
2301  bool is_master;
2302  bool is_slave;
2303  const mesh_im *mim;
2304 
2305  std::string expr;
2306  };
2307 
2308  std::vector<contact_boundary> boundaries;
2309  std::string transformation_name;
2310  std::string u_group;
2311  std::string w_group;
2312  std::string friction_coeff;
2313  std::string alpha;
2314  std::string augmentation_param;
2315  model::varnamelist vl, dl;
2316  model::mimlist ml;
2317 
2318  bool sym_version, frame_indifferent;
2319 
2320  void add_contact_boundary(model &md, const mesh_im &mim, size_type region,
2321  bool is_master, bool is_slave,
2322  const std::string &u,
2323  const std::string &lambda,
2324  const std::string &w = "") {
2325  std::string test_u = "Test_" + sup_previous_and_dot_to_varname(u);
2326  std::string test_u_group = "Test_" + sup_previous_and_dot_to_varname(u_group);
2327  std::string test_lambda = "Test_" + sup_previous_and_dot_to_varname(lambda);
2328  GMM_ASSERT1(is_slave || is_master, "The contact boundary should be "
2329  "either master, slave or both");
2330  const mesh_fem *mf = md.pmesh_fem_of_variable(u);
2331  GMM_ASSERT1(mf, "The displacement variable should be a f.e.m. one");
2332  GMM_ASSERT1(&(mf->linked_mesh()) == &(mim.linked_mesh()),
2333  "The displacement variable and the integration method "
2334  "should share the same mesh");
2335  if (is_slave) {
2336  const mesh_fem *mf_l = md.pmesh_fem_of_variable(lambda);
2337  GMM_ASSERT1(mf, "The multiplier variable should be a f.e.m. one");
2338  GMM_ASSERT1(&(mf_l->linked_mesh()) == &(mim.linked_mesh()),
2339  "The displacement variable and the multiplier one "
2340  "should share the same mesh");
2341  }
2342 
2343  if (w.size()) {
2344  const mesh_fem *mf2 = md.pmesh_fem_of_variable(w);
2345  GMM_ASSERT1(!mf2 || &(mf2->linked_mesh()) == &(mf->linked_mesh()),
2346  "The data for the sliding velocity should be defined on "
2347  " the same mesh as the displacement variable");
2348  }
2349 
2350  for (size_type i = 0; i < boundaries.size(); ++i) {
2351  const contact_boundary &cb = boundaries[i];
2352  if (&(md.mesh_fem_of_variable(cb.varname_u).linked_mesh())
2353  == &(mf->linked_mesh()) && cb.region == region)
2354  GMM_ASSERT1(false, "This contact boundary has already been added");
2355  }
2356  if (is_master)
2357  add_master_contact_boundary_to_raytracing_transformation
2358  (md, transformation_name, mf->linked_mesh(), u_group, region);
2359  else
2360  add_slave_contact_boundary_to_raytracing_transformation
2361  (md, transformation_name, mf->linked_mesh(), u_group, region);
2362 
2363  boundaries.push_back(contact_boundary());
2364  contact_boundary &cb = boundaries.back();
2365  cb.region = region;
2366  cb.varname_u = u;
2367  if (is_slave) cb.varname_lambda = lambda;
2368  cb.varname_w = w;
2369  cb.is_master = is_master;
2370  cb.is_slave = is_slave;
2371  cb.mim = &mim;
2372  if (is_slave) {
2373  std::string n, n0, Vs, g, Y;
2374  n = "Transformed_unit_vector(Grad_"+u+", Normal)";
2375  n0 = "Transformed_unit_vector(Grad_"+w+", Normal)";
2376 
2377  // Coulomb_friction_coupled_projection(lambda,
2378  // Transformed_unit_vector(Grad_u, Normal),
2379  // (u-Interpolate(ug,trans)-(w-Interpolate(wg,trans)))*alpha,
2380  // (Interpolate(X,trans)+Interpolate(ug,trans)-X-u).
2381  // Transformed_unit_vector(Grad_u, Normal), f, r)
2382  Y = "Interpolate(X,"+transformation_name+")";
2383  g = "("+Y+"+Interpolate("+u_group+","+transformation_name+")-X-"+u+")."+n;
2384  Vs = "("+u+"-Interpolate("+u_group+","+transformation_name+")";
2385  if (w.size()) {
2386  Vs += "-"+w+"+Interpolate("+w_group+","+transformation_name+")";
2387  if (frame_indifferent)
2388  Vs += "+("+g+")*("+n+"-"+n0+")";
2389  }
2390  Vs += ")*"+alpha;
2391 
2392  std::string coupled_projection_def =
2393  "Coulomb_friction_coupled_projection("
2394  + lambda+","+n+","+Vs+","+g+","+friction_coeff+","
2395  + augmentation_param+")";
2396 
2397  // Coulomb_friction_coupled_projection(lambda,
2398  // Transformed_unit_vector(Grad_u, Normal), (u-w)*alpha,
2399  // (Interpolate(X,trans)-X-u).Transformed_unit_vector(Grad_u,
2400  // Normal), f, r)
2401  g = "(Interpolate(X,"+transformation_name+")-X-"+u+")."+n;
2402  if (frame_indifferent && w.size())
2403  Vs = "("+u+"-"+w+"+"+g+"*("+n+"-"+n0+"))*"+alpha;
2404  else
2405  Vs = "("+u+(w.size() ? ("-"+w):"")+")*"+alpha;
2406 
2407  std::string coupled_projection_rig =
2408  "Coulomb_friction_coupled_projection("
2409  + lambda+","+n+","+Vs+","+g+","+ friction_coeff+","
2410  + augmentation_param+")";
2411 
2412  cb.expr =
2413  // -lambda.Test_u for non-symmetric version
2414  (sym_version ? "" : ("-"+lambda+"."+test_u))
2415  // -coupled_projection_def.Test_u and -coupled_projection_rig.Test_u
2416  // for symmetric version
2417  + (sym_version ? ("+ Interpolate_filter("+transformation_name+",-"
2418  +coupled_projection_def+"."+test_u+",1)") : "")
2419  + (sym_version ? ("+ Interpolate_filter("+transformation_name+",-"
2420  +coupled_projection_rig+"."+test_u+",2)") : "")
2421  // Interpolate_filter(trans,
2422  // lambda.Interpolate(Test_ug, contact_trans), 1)
2423  // or
2424  // Interpolate_filter(trans,
2425  // coupled_projection_def.Interpolate(Test_ug, contact_trans), 1)
2426  + "+ Interpolate_filter("+transformation_name+","
2427  + (sym_version ? coupled_projection_def : lambda)
2428  + ".Interpolate("+test_u_group+"," + transformation_name+"), 1)"
2429  // -(1/r)*lambda.Test_lambda
2430  + "-(1/"+augmentation_param+")*"+lambda+"."+test_lambda
2431  // Interpolate_filter(trans,
2432  // (1/r)*coupled_projection_rig.Test_lambda, 2)
2433  + "+ Interpolate_filter("+transformation_name+","
2434  + "(1/"+augmentation_param+")*"+ coupled_projection_rig
2435  + "."+test_lambda+", 2)"
2436  // Interpolate_filter(trans,
2437  // (1/r)*coupled_projection_def.Test_lambda, 1)
2438  + "+ Interpolate_filter("+transformation_name+","
2439  + "(1/"+augmentation_param+")*" + coupled_projection_def + "."
2440  + test_lambda+", 1)";
2441  }
2442  }
2443 
2444  virtual void asm_real_tangent_terms(const model &md, size_type ,
2445  const model::varnamelist &,
2446  const model::varnamelist &,
2447  const model::mimlist &,
2448  model::real_matlist &,
2449  model::real_veclist &,
2450  model::real_veclist &,
2451  size_type,
2452  build_version) const {
2453  // GMM_ASSERT1(mims.size() == 1,
2454  // "Generic linear assembly brick needs one and only one "
2455  // "mesh_im"); // to be verified ...
2456 
2457  for (const contact_boundary &cb : boundaries) {
2458  if (cb.is_slave)
2459  md.add_generic_expression(cb.expr, *(cb.mim), cb.region);
2460  }
2461  }
2462 
2463 
2464  intergral_large_sliding_contact_brick_raytracing
2465  (const std::string &r,
2466  const std::string &f_coeff,const std::string &ug,
2467  const std::string &wg, const std::string &tr,
2468  const std::string &alpha_ = "1", bool sym_v = false,
2469  bool frame_indiff = false) {
2470  transformation_name = tr;
2471  u_group = ug; w_group = wg;
2472  friction_coeff = f_coeff;
2473  alpha = alpha_;
2474  augmentation_param = r;
2475  sym_version = sym_v;
2476  frame_indifferent = frame_indiff;
2477 
2478  set_flags("Integral large sliding contact bick raytracing",
2479  false /* is linear*/,
2480  false /* is symmetric */, false /* is coercive */,
2481  true /* is real */, false /* is complex */);
2482  }
2483 
2484  };
2485 
2486 
2487  const std::string &transformation_name_of_large_sliding_contact_brick
2488  (model &md, size_type indbrick) {
2489  pbrick pbr = md.brick_pointer(indbrick);
2490  intergral_large_sliding_contact_brick_raytracing *p
2491  = dynamic_cast<intergral_large_sliding_contact_brick_raytracing *>
2492  (const_cast<virtual_brick *>(pbr.get()));
2493  GMM_ASSERT1(p, "Wrong type of brick");
2494  return p->transformation_name;
2495  }
2496 
2497  const std::string &displacement_group_name_of_large_sliding_contact_brick
2498  (model &md, size_type indbrick) {
2499  pbrick pbr = md.brick_pointer(indbrick);
2500  intergral_large_sliding_contact_brick_raytracing *p
2501  = dynamic_cast<intergral_large_sliding_contact_brick_raytracing *>
2502  (const_cast<virtual_brick *>(pbr.get()));
2503  GMM_ASSERT1(p, "Wrong type of brick");
2504  return p->u_group;
2505  }
2506 
2507  const std::string &sliding_data_group_name_of_large_sliding_contact_brick
2508  (model &md, size_type indbrick) {
2509  pbrick pbr = md.brick_pointer(indbrick);
2510  intergral_large_sliding_contact_brick_raytracing *p
2511  = dynamic_cast<intergral_large_sliding_contact_brick_raytracing *>
2512  (const_cast<virtual_brick *>(pbr.get()));
2513  GMM_ASSERT1(p, "Wrong type of brick");
2514  return p->w_group;
2515  }
2516 
2517  void add_rigid_obstacle_to_large_sliding_contact_brick
2518  (model &md, size_type indbrick, const std::string &expr, size_type N) {
2519  pbrick pbr = md.brick_pointer(indbrick);
2520  intergral_large_sliding_contact_brick_raytracing *p
2521  = dynamic_cast<intergral_large_sliding_contact_brick_raytracing *>
2522  (const_cast<virtual_brick *>(pbr.get()));
2523  GMM_ASSERT1(p, "Wrong type of brick");
2524  add_rigid_obstacle_to_raytracing_transformation
2525  (md, p->transformation_name, expr, N);
2526  }
2527 
2528  void add_contact_boundary_to_large_sliding_contact_brick
2529  (model &md, size_type indbrick, const mesh_im &mim, size_type region,
2530  bool is_master, bool is_slave, const std::string &u,
2531  const std::string &lambda, const std::string &w) {
2532 
2533  pbrick pbr = md.brick_pointer(indbrick);
2534  intergral_large_sliding_contact_brick_raytracing *p
2535  = dynamic_cast<intergral_large_sliding_contact_brick_raytracing *>
2536  (const_cast<virtual_brick *>(pbr.get()));
2537  GMM_ASSERT1(p, "Wrong type of brick");
2538 
2539  bool found_u = false, found_lambda = false;
2540  for (const auto & v : p->vl) {
2541  if (v.compare(u) == 0) found_u = true;
2542  if (v.compare(lambda) == 0) found_lambda = true;
2543  }
2544  if (!found_u) p->vl.push_back(u);
2545  GMM_ASSERT1(!is_slave || lambda.size(),
2546  "You should define a multiplier on each slave boundary");
2547  if (is_slave && !found_lambda) p->vl.push_back(lambda);
2548  if (!found_u || (is_slave && !found_lambda))
2549  md.change_variables_of_brick(indbrick, p->vl);
2550 
2551  std::vector<std::string> ug = md.variable_group(p->u_group);
2552  found_u = false;
2553  for (const auto &uu : ug)
2554  if (uu.compare(u) == 0) { found_u = true; break; }
2555  if (!found_u) {
2556  ug.push_back(u);
2557  md.define_variable_group(p->u_group, ug);
2558  }
2559 
2560  if (w.size()) {
2561  bool found_w = false;
2562  for (const auto &ww : p->dl)
2563  if (ww.compare(w) == 0) { found_w = true; break; }
2564  if (!found_w) {
2565  p->dl.push_back(w);
2566  md.change_data_of_brick(indbrick, p->dl);
2567  }
2568  std::vector<std::string> wg = md.variable_group(p->w_group);
2569  found_w = false;
2570  for (const auto &ww : wg)
2571  if (ww.compare(w) == 0) { found_w = true; break; }
2572  if (!found_w) {
2573  wg.push_back(w);
2574  md.define_variable_group(p->w_group, wg);
2575  }
2576  }
2577 
2578  bool found_mim = false;
2579  for (const auto &pmim : p->ml)
2580  if (pmim == &mim) { found_mim = true; break; }
2581  if (!found_mim) {
2582  p->ml.push_back(&mim);
2583  md.change_mims_of_brick(indbrick, p->ml);
2584  }
2585 
2586  p->add_contact_boundary(md, mim, region, is_master, is_slave,
2587  u, lambda, w);
2588  }
2589 
2590  size_type add_integral_large_sliding_contact_brick_raytracing
2591  (model &md, const std::string &augm_param,
2592  scalar_type release_distance, const std::string &f_coeff,
2593  const std::string &alpha, bool sym_v, bool frame_indifferent) {
2594 
2595  char ugroupname[50], wgroupname[50], transname[50];
2596  for (int i = 0; i < 10000; ++i) {
2597  sprintf(ugroupname, "disp__group_raytracing_%d", i);
2598  if (!(md.variable_group_exists(ugroupname))
2599  && !(md.variable_exists(ugroupname)))
2600  break;
2601  }
2602  md.define_variable_group(ugroupname, std::vector<std::string>());
2603 
2604  for (int i = 0; i < 10000; ++i) {
2605  sprintf(wgroupname, "w__group_raytracing_%d", i);
2606  if (!(md.variable_group_exists(wgroupname))
2607  && !(md.variable_exists(wgroupname)))
2608  break;
2609  }
2610  md.define_variable_group(wgroupname, std::vector<std::string>());
2611 
2612  for (int i = 0; i < 10000; ++i) {
2613  sprintf(transname, "trans__raytracing_%d", i);
2614  if (!(md.interpolate_transformation_exists(transname)))
2615  break;
2616  }
2617  add_raytracing_transformation(md, transname, release_distance);
2618 
2619  model::varnamelist vl, dl;
2620  if (md.variable_exists(augm_param)) dl.push_back(augm_param);
2621  if (md.variable_exists(f_coeff)) dl.push_back(f_coeff);
2622  if (md.variable_exists(alpha)) dl.push_back(alpha);
2623 
2624  auto p = std::make_shared<intergral_large_sliding_contact_brick_raytracing>
2625  (augm_param, f_coeff, ugroupname, wgroupname, transname, alpha,
2626  sym_v, frame_indifferent);
2627  pbrick pbr(p);
2628  p->dl = dl;
2629 
2630  return md.add_brick(pbr, p->vl, p->dl, model::termlist(), model::mimlist(),
2631  size_type(-1));
2632  }
2633 
2634 
2635 
2636  struct Nitsche_large_sliding_contact_brick_raytracing
2637  : public virtual_brick {
2638 
2639  struct contact_boundary {
2640  size_type region;
2641  std::string varname_u;
2642  std::string sigma_u;
2643  std::string varname_w;
2644  bool is_master;
2645  bool is_slave;
2646  bool is_unbiased;
2647  const mesh_im *mim;
2648 
2649  std::string expr;
2650  };
2651 
2652  std::vector<contact_boundary> boundaries;
2653  std::string transformation_name;
2654  std::string u_group;
2655  std::string w_group;
2656  std::string friction_coeff;
2657  std::string alpha;
2658  std::string Nitsche_param;
2659  model::varnamelist vl, dl;
2660  model::mimlist ml;
2661 
2662  bool sym_version, frame_indifferent, unbiased;
2663 
2664  void add_contact_boundary(model &md, const mesh_im &mim, size_type region,
2665  bool is_master, bool is_slave, bool is_unbiased,
2666  const std::string &u,
2667  const std::string &sigma_u,
2668  const std::string &w = "") {
2669  std::string test_u = "Test_" + sup_previous_and_dot_to_varname(u);
2670  std::string test_u_group = "Test_" + sup_previous_and_dot_to_varname(u_group);
2671  GMM_ASSERT1(is_slave || is_master, "The contact boundary should be "
2672  "either master, slave or both");
2673  if (is_unbiased) {
2674  GMM_ASSERT1((is_slave && is_master), "The contact boundary should be "
2675  "both master and slave for the unbiased version");
2676  is_slave=true; is_master=true;
2677  }
2678  const mesh_fem *mf = md.pmesh_fem_of_variable(u);
2679  GMM_ASSERT1(mf, "The displacement variable should be a f.e.m. one");
2680  GMM_ASSERT1(&(mf->linked_mesh()) == &(mim.linked_mesh()),
2681  "The displacement variable and the integration method "
2682  "should share the same mesh");
2683 
2684  if (w.size()) {
2685  const mesh_fem *mf2 = md.pmesh_fem_of_variable(w);
2686  GMM_ASSERT1(!mf2 || &(mf2->linked_mesh()) == &(mf->linked_mesh()),
2687  "The data for the sliding velocity should be defined on "
2688  " the same mesh as the displacement variable");
2689  }
2690 
2691  for (size_type i = 0; i < boundaries.size(); ++i) {
2692  const contact_boundary &cb = boundaries[i];
2693  if (&(md.mesh_fem_of_variable(cb.varname_u).linked_mesh())
2694  == &(mf->linked_mesh()) && cb.region == region)
2695  GMM_ASSERT1(false, "This contact boundary has already been added");
2696  }
2697  if (is_master)
2698  add_master_contact_boundary_to_raytracing_transformation
2699  (md, transformation_name, mf->linked_mesh(), u_group, region);
2700  else
2701  add_slave_contact_boundary_to_raytracing_transformation
2702  (md, transformation_name, mf->linked_mesh(), u_group, region);
2703 
2704  boundaries.push_back(contact_boundary());
2705  contact_boundary &cb = boundaries.back();
2706  cb.region = region;
2707  cb.varname_u = u;
2708  if (is_slave) cb.sigma_u = sigma_u;
2709  cb.varname_w = w;
2710  cb.is_master = is_master;
2711  cb.is_slave = is_slave;
2712  cb.is_unbiased= is_unbiased;
2713  cb.mim = &mim;
2714  if (is_slave) {
2715  std::string n, n0, Vs, g, Y, gamma;
2716 
2717  gamma ="("+Nitsche_param+"/element_size)";
2718  n = "Transformed_unit_vector(Grad_"+u+", Normal)";
2719  n0 = "Transformed_unit_vector(Grad_"+w+", Normal)";
2720 
2721  // For deformable bodies:
2722  // Coulomb_friction_coupled_projection(sigma(u),
2723  // Transformed_unit_vector(Grad_u, Normal),
2724  // (u-Interpolate(ug,trans)-(w-Interpolate(wg,trans)))*alpha,
2725  // (Interpolate(X,trans)+Interpolate(ug,trans)-X-u).
2726  // Transformed_unit_vector(Grad_u, Normal), f, r)
2727  Y = "Interpolate(X,"+transformation_name+")";
2728  g = "("+Y+"+Interpolate("+u_group+","+transformation_name+")-X-"+u+")."+n;
2729  Vs = "("+u+"-Interpolate("+u_group+","+transformation_name+")";
2730  if (w.size()) {
2731  Vs += "-"+w+"+Interpolate("+w_group+","+transformation_name+")";
2732  if (frame_indifferent)
2733  Vs += "+("+g+")*("+n+"-"+n0+")";
2734  }
2735  Vs += ")*"+alpha;
2736 
2737  std::string coupled_projection_def =
2738  "Coulomb_friction_coupled_projection("
2739  + sigma_u +","+n+","+Vs+","+g+","+friction_coeff+","
2740  + gamma +")";
2741 
2742  // For regid obstacle:
2743  // Coulomb_friction_coupled_projection(sigma(u),
2744  // Transformed_unit_vector(Grad_u, Normal), (u-w)*alpha,
2745  // (Interpolate(X,trans)-X-u).Transformed_unit_vector(Grad_u,
2746  // Normal), f, r)
2747  g = "(Interpolate(X,"+transformation_name+")-X-"+u+")."+n;
2748  if (frame_indifferent && w.size())
2749  Vs = "("+u+"-"+w+"+"+g+"*("+n+"-"+n0+"))*"+alpha;
2750  else
2751  Vs = "("+u+(w.size() ? ("-"+w):"")+")*"+alpha;
2752 
2753  std::string coupled_projection_rig =
2754  "Coulomb_friction_coupled_projection("
2755  + sigma_u +","+n+","+Vs+","+g+","+ friction_coeff+","
2756  + gamma +")";
2757 
2758  cb.expr =
2759  // 0.5* for non-biaised version
2760  (is_unbiased ? "-0.5*" : "-")
2761  // -coupled_projection_def.Test_u
2762  + ("Interpolate_filter("+transformation_name+","
2763  +coupled_projection_def+"."+test_u+",1) ")
2764  +(is_unbiased ? "":"-Interpolate_filter("+transformation_name+","
2765  +coupled_projection_rig+"."+test_u+",2) ")
2766  // Interpolate_filter(trans,
2767  // lambda.Interpolate(Test_ug, contact_trans), 1)
2768  // or
2769  // Interpolate_filter(trans,
2770  // coupled_projection_def.Interpolate(Test_ug, contact_trans), 1)
2771  + (is_unbiased ? "+ 0.5*" : "+ ")
2772  +"Interpolate_filter("+transformation_name+","
2773  +coupled_projection_def +".Interpolate("+test_u_group+"," + transformation_name+"), 1)"
2774  +(is_unbiased ? "":"+ Interpolate_filter("+transformation_name+","
2775  +coupled_projection_rig +".Interpolate("+test_u_group+"," + transformation_name+"), 2)");
2776  }
2777  }
2778 
2779  virtual void asm_real_tangent_terms(const model &md, size_type ,
2780  const model::varnamelist &,
2781  const model::varnamelist &,
2782  const model::mimlist &,
2783  model::real_matlist &,
2784  model::real_veclist &,
2785  model::real_veclist &,
2786  size_type,
2787  build_version) const {
2788  // GMM_ASSERT1(mims.size() == 1,
2789  // "Generic linear assembly brick needs one and only one "
2790  // "mesh_im"); // to be verified ...
2791 
2792  for (const contact_boundary &cb : boundaries) {
2793  if (cb.is_slave)
2794  md.add_generic_expression(cb.expr, *(cb.mim), cb.region);
2795  }
2796  }
2797 
2798 
2799  Nitsche_large_sliding_contact_brick_raytracing
2800  ( bool unbia,
2801  const std::string &Nitsche_parameter,
2802  const std::string &f_coeff,const std::string &ug,
2803  const std::string &wg, const std::string &tr,
2804  const std::string &alpha_ = "1", bool sym_v = false,
2805  bool frame_indiff = false) {
2806  transformation_name = tr;
2807  u_group = ug; w_group = wg;
2808  friction_coeff = f_coeff;
2809  alpha = alpha_;
2810  Nitsche_param = Nitsche_parameter;
2811  sym_version = sym_v;
2812  frame_indifferent = frame_indiff;
2813  unbiased = unbia;
2814 
2815  set_flags("Integral large sliding contact bick raytracing",
2816  false /* is linear*/,
2817  false /* is symmetric */, false /* is coercive */,
2818  true /* is real */, false /* is complex */);
2819  }
2820 
2821  };
2822 
2823 
2824  const std::string &transformation_name_of_Nitsche_large_sliding_contact_brick
2825  (model &md, size_type indbrick) {
2826  pbrick pbr = md.brick_pointer(indbrick);
2827  Nitsche_large_sliding_contact_brick_raytracing *p
2828  = dynamic_cast<Nitsche_large_sliding_contact_brick_raytracing *>
2829  (const_cast<virtual_brick *>(pbr.get()));
2830  GMM_ASSERT1(p, "Wrong type of brick");
2831  return p->transformation_name;
2832  }
2833 
2834  const std::string &displacement_group_name_of_Nitsche_large_sliding_contact_brick
2835  (model &md, size_type indbrick) {
2836  pbrick pbr = md.brick_pointer(indbrick);
2837  Nitsche_large_sliding_contact_brick_raytracing *p
2838  = dynamic_cast<Nitsche_large_sliding_contact_brick_raytracing *>
2839  (const_cast<virtual_brick *>(pbr.get()));
2840  GMM_ASSERT1(p, "Wrong type of brick");
2841  return p->u_group;
2842  }
2843 
2844  const std::string &sliding_data_group_name_of_Nitsche_large_sliding_contact_brick
2845  (model &md, size_type indbrick) {
2846  pbrick pbr = md.brick_pointer(indbrick);
2847  Nitsche_large_sliding_contact_brick_raytracing *p
2848  = dynamic_cast<Nitsche_large_sliding_contact_brick_raytracing *>
2849  (const_cast<virtual_brick *>(pbr.get()));
2850  GMM_ASSERT1(p, "Wrong type of brick");
2851  return p->w_group;
2852  }
2853 
2854  void add_rigid_obstacle_to_Nitsche_large_sliding_contact_brick
2855  (model &md, size_type indbrick, const std::string &expr, size_type N) {
2856  pbrick pbr = md.brick_pointer(indbrick);
2857  Nitsche_large_sliding_contact_brick_raytracing *p
2858  = dynamic_cast<Nitsche_large_sliding_contact_brick_raytracing *>
2859  (const_cast<virtual_brick *>(pbr.get()));
2860  GMM_ASSERT1(p, "Wrong type of brick");
2861  add_rigid_obstacle_to_raytracing_transformation
2862  (md, p->transformation_name, expr, N);
2863  }
2864 
2865  void add_contact_boundary_to_Nitsche_large_sliding_contact_brick
2866  (model &md, size_type indbrick, const mesh_im &mim, size_type region,
2867  bool is_master, bool is_slave, bool is_unbiased, const std::string &u,
2868  const std::string &sigma_u, const std::string &w) {
2869 
2870  pbrick pbr = md.brick_pointer(indbrick);
2871  Nitsche_large_sliding_contact_brick_raytracing *p
2872  = dynamic_cast<Nitsche_large_sliding_contact_brick_raytracing *>
2873  (const_cast<virtual_brick *>(pbr.get()));
2874  GMM_ASSERT1(p, "Wrong type of brick");
2875 
2876  bool found_u = false, found_sigma = false;
2877  for (const auto & v : p->vl) {
2878  if (v.compare(u) == 0) found_u = true;
2879  if (v.compare(sigma_u) == 0) found_sigma = true;
2880  }
2881  if (!found_u) p->vl.push_back(u);
2882  GMM_ASSERT1(!is_slave || sigma_u.size(),
2883  "You should define an expression of sigma(u) on each slave boundary");
2884  // if (is_slave && !found_sigma) p->vl.push_back(sigma_u);
2885  if (!found_u)
2886  md.change_variables_of_brick(indbrick, p->vl);
2887 
2888  std::vector<std::string> ug = md.variable_group(p->u_group);
2889  found_u = false;
2890  for (const auto &uu : ug)
2891  if (uu.compare(u) == 0) { found_u = true; break; }
2892  if (!found_u) {
2893  ug.push_back(u);
2894  md.define_variable_group(p->u_group, ug);
2895  }
2896 
2897  if (w.size()) {
2898  bool found_w = false;
2899  for (const auto &ww : p->dl)
2900  if (ww.compare(w) == 0) { found_w = true; break; }
2901  if (!found_w) {
2902  p->dl.push_back(w);
2903  md.change_data_of_brick(indbrick, p->dl);
2904  }
2905  std::vector<std::string> wg = md.variable_group(p->w_group);
2906  found_w = false;
2907  for (const auto &ww : wg)
2908  if (ww.compare(w) == 0) { found_w = true; break; }
2909  if (!found_w) {
2910  wg.push_back(w);
2911  md.define_variable_group(p->w_group, wg);
2912  }
2913  }
2914 
2915  bool found_mim = false;
2916  for (const auto &pmim : p->ml)
2917  if (pmim == &mim) { found_mim = true; break; }
2918  if (!found_mim) {
2919  p->ml.push_back(&mim);
2920  md.change_mims_of_brick(indbrick, p->ml);
2921  }
2922 
2923  p->add_contact_boundary(md, mim, region, is_master, is_slave,is_unbiased,
2924  u, sigma_u, w);
2925  }
2926 
2927  size_type add_Nitsche_large_sliding_contact_brick_raytracing
2928  (model &md, bool is_unbiased, const std::string &Nitsche_param,
2929  scalar_type release_distance, const std::string &f_coeff,
2930  const std::string &alpha, bool sym_v, bool frame_indifferent) {
2931 
2932  char ugroupname[50], wgroupname[50], transname[50];
2933  for (int i = 0; i < 10000; ++i) {
2934  sprintf(ugroupname, "disp__group_raytracing_%d", i);
2935  if (!(md.variable_group_exists(ugroupname))
2936  && !(md.variable_exists(ugroupname)))
2937  break;
2938  }
2939  md.define_variable_group(ugroupname, std::vector<std::string>());
2940 
2941  for (int i = 0; i < 10000; ++i) {
2942  sprintf(wgroupname, "w__group_raytracing_%d", i);
2943  if (!(md.variable_group_exists(wgroupname))
2944  && !(md.variable_exists(wgroupname)))
2945  break;
2946  }
2947  md.define_variable_group(wgroupname, std::vector<std::string>());
2948 
2949  for (int i = 0; i < 10000; ++i) {
2950  sprintf(transname, "trans__raytracing_%d", i);
2951  if (!(md.interpolate_transformation_exists(transname)))
2952  break;
2953  }
2954  add_raytracing_transformation(md, transname, release_distance);
2955 
2956  model::varnamelist vl, dl;
2957  if (md.variable_exists(Nitsche_param)) dl.push_back(Nitsche_param);
2958  if (md.variable_exists(f_coeff)) dl.push_back(f_coeff);
2959  if (md.variable_exists(alpha)) dl.push_back(alpha);
2960 
2961  auto p = std::make_shared<Nitsche_large_sliding_contact_brick_raytracing>
2962  (is_unbiased, Nitsche_param , f_coeff, ugroupname, wgroupname, transname, alpha,
2963  sym_v, frame_indifferent);
2964  pbrick pbr(p);
2965  p->dl = dl;
2966 
2967  return md.add_brick(pbr, p->vl, p->dl, model::termlist(), model::mimlist(),
2968  size_type(-1));
2969  }
2970 
2971 }
2972 
2973 
2974 
2975 
2976 
2977  /* end of namespace getfem. */
getfem::add_rigid_obstacle_to_Nitsche_large_sliding_contact_brick
void add_rigid_obstacle_to_Nitsche_large_sliding_contact_brick(model &md, size_type indbrick, const std::string &expr, size_type N)
Adds a rigid obstacle to an existing large sliding contact with friction brick.
Definition: getfem_contact_and_friction_large_sliding.cc:2855
gmm::lu_factor
size_type lu_factor(DenseMatrix &A, lapack_ipvt &ipvt)
LU Factorization of a general (dense) matrix (real or complex).
Definition: gmm_dense_lu.h:129
getfem::model::variable_exists
bool variable_exists(const std::string &name) const
States if a name corresponds to a declared variable.
Definition: getfem_models.cc:947
bgeot::compute_normal
base_small_vector compute_normal(const geotrans_interpolation_context &c, size_type face)
norm of returned vector is the ratio between the face surface on the real element and the face surfac...
Definition: bgeot_geometric_trans.cc:1082
bgeot::geotrans_interpolation_context::xreal
const base_node & xreal() const
coordinates of the current point, in the real convex.
Definition: bgeot_geometric_trans.cc:299
getfem::displacement_group_name_of_large_sliding_contact_brick
const std::string & displacement_group_name_of_large_sliding_contact_brick(model &md, size_type indbrick)
Gives the name of the group of variables corresponding to the displacement for an existing large slid...
Definition: getfem_contact_and_friction_large_sliding.cc:2498
getfem::add_Nitsche_large_sliding_contact_brick_raytracing
size_type add_Nitsche_large_sliding_contact_brick_raytracing(model &md, bool is_unbiased, const std::string &Nitsche_param, scalar_type release_distance, const std::string &f_coeff="0", const std::string &alpha="1", bool sym_v=false, bool frame_indifferent=false)
Adds a large sliding contact with friction brick to the model based on Nitsche's method.
Definition: getfem_contact_and_friction_large_sliding.cc:2928
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::add_contact_boundary_to_Nitsche_large_sliding_contact_brick
void add_contact_boundary_to_Nitsche_large_sliding_contact_brick(model &md, size_type indbrick, const mesh_im &mim, size_type region, bool is_master, bool is_slave, bool is_unbiased, const std::string &u, const std::string &lambda="", const std::string &w="")
Adds a contact boundary to an existing large sliding contact with friction brick.
Definition: getfem_contact_and_friction_large_sliding.cc:2866
getfem::sliding_data_group_name_of_Nitsche_large_sliding_contact_brick
const std::string & sliding_data_group_name_of_Nitsche_large_sliding_contact_brick(model &md, size_type indbrick)
Gives the name of the group of variables corresponding to the sliding data for an existing large slid...
Definition: getfem_contact_and_friction_large_sliding.cc:2845
getfem::add_integral_large_sliding_contact_brick_raytracing
size_type add_integral_large_sliding_contact_brick_raytracing(model &md, const std::string &augm_param, scalar_type release_distance, const std::string &f_coeff="0", const std::string &alpha="1", bool sym_v=false, bool frame_indifferent=false)
Adds a large sliding contact with friction brick to the model.
Definition: getfem_contact_and_friction_large_sliding.cc:2591
getfem::add_contact_boundary_to_large_sliding_contact_brick
void add_contact_boundary_to_large_sliding_contact_brick(model &md, size_type indbrick, const mesh_im &mim, size_type region, bool is_master, bool is_slave, const std::string &u, const std::string &lambda="", const std::string &w="")
Adds a contact boundary to an existing large sliding contact with friction brick.
Definition: getfem_contact_and_friction_large_sliding.cc:2529
getfem::add_rigid_obstacle_to_large_sliding_contact_brick
void add_rigid_obstacle_to_large_sliding_contact_brick(model &md, size_type indbrick, const std::string &expr, size_type N)
Adds a rigid obstacle to an existing large sliding contact with friction brick.
Definition: getfem_contact_and_friction_large_sliding.cc:2518
getfem::asm_mass_matrix
void asm_mass_matrix(const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_region &rg=mesh_region::all_convexes())
generic mass matrix assembly (on the whole mesh or on the specified convex set or boundary)
Definition: getfem_assembling.h:697
gmm::vect_dist2
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist2(const V1 &v1, const V2 &v2)
Euclidean distance between two vectors.
Definition: gmm_blas.h:597
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::fem_interpolation_context::convex_num
size_type convex_num() const
get the current convex number
Definition: getfem_fem.cc:52
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
getfem::add_slave_contact_boundary_to_raytracing_transformation
void add_slave_contact_boundary_to_raytracing_transformation(model &md, const std::string &transname, const mesh &m, const std::string &dispname, size_type region)
Add a slave boundary with corresponding displacement variable 'dispname' on a specific boundary 'regi...
Definition: getfem_contact_and_friction_common.cc:2450
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::rtree
Balanced tree of n-dimensional rectangles.
Definition: bgeot_rtree.h:98
getfem::model::change_mims_of_brick
void change_mims_of_brick(size_type ib, const mimlist &ml)
Change the mim list of a brick.
Definition: getfem_models.cc:1120
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
bgeot_rtree.h
region-tree for window/point search on a set of rectangles.
getfem::model::interpolate_transformation_exists
bool interpolate_transformation_exists(const std::string &name) const
Tests if name corresponds to an interpolate transformation.
Definition: getfem_models.h:1129
getfem::fem_interpolation_context::face_num
short_type face_num() const
get the current face number
Definition: getfem_fem.cc:61
getfem::displacement_group_name_of_Nitsche_large_sliding_contact_brick
const std::string & displacement_group_name_of_Nitsche_large_sliding_contact_brick(model &md, size_type indbrick)
Gives the name of the group of variables corresponding to the displacement for an existing large slid...
Definition: getfem_contact_and_friction_large_sliding.cc:2835
getfem::mesh_fem::nb_basic_dof
virtual size_type nb_basic_dof() const
Return the total number of basic degrees of freedom (before the optional reduction).
Definition: getfem_mesh_fem.h:556
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
getfem::asm_source_term
void asm_source_term(const VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg=mesh_region::all_convexes())
source term (for both volumic sources and boundary (Neumann) sources).
Definition: getfem_assembling.h:877
getfem::add_raytracing_transformation
void add_raytracing_transformation(model &md, const std::string &transname, scalar_type release_distance)
Add a raytracing interpolate transformation called 'transname' to a model to be used by the generic a...
Definition: getfem_contact_and_friction_common.cc:2426
getfem::add_integral_large_sliding_contact_brick_raytrace
size_type add_integral_large_sliding_contact_brick_raytrace(model &md, multi_contact_frame &mcf, const std::string &dataname_r, const std::string &dataname_friction_coeff=std::string(), const std::string &dataname_alpha=std::string())
Adds a large sliding contact with friction brick to the model.
Definition: getfem_contact_and_friction_large_sliding.cc:1137
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
bgeot::geotrans_interpolation_context::set_ii
void set_ii(size_type ii__)
change the current point (assuming a geotrans_precomp_ is used)
Definition: bgeot_geometric_trans.h:463
getfem_contact_and_friction_integral.h
Unilateral contact and Coulomb friction condition brick.
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem::model::add_brick
size_type add_brick(pbrick pbr, const varnamelist &varnames, const varnamelist &datanames, const termlist &terms, const mimlist &mims, size_type region)
Add a brick to the model.
Definition: getfem_models.cc:1033
getfem::model::change_data_of_brick
void change_data_of_brick(size_type ib, const varnamelist &vl)
Change the data list of a brick.
Definition: getfem_models.cc:1112
getfem::add_master_contact_boundary_to_raytracing_transformation
void add_master_contact_boundary_to_raytracing_transformation(model &md, const std::string &transname, const mesh &m, const std::string &dispname, size_type region)
Add a master boundary with corresponding displacement variable 'dispname' on a specific boundary 'reg...
Definition: getfem_contact_and_friction_common.cc:2440
getfem::model::change_variables_of_brick
void change_variables_of_brick(size_type ib, const varnamelist &vl)
Change the variable list of a brick.
Definition: getfem_models.cc:1104
getfem::transformation_name_of_Nitsche_large_sliding_contact_brick
const std::string & transformation_name_of_Nitsche_large_sliding_contact_brick(model &md, size_type indbrick)
Gives the name of the raytracing interpolate transformation for an existing large sliding contact bri...
Definition: getfem_contact_and_friction_large_sliding.cc:2825
getfem::sliding_data_group_name_of_large_sliding_contact_brick
const std::string & sliding_data_group_name_of_large_sliding_contact_brick(model &md, size_type indbrick)
Gives the name of the group of variables corresponding to the sliding data for an existing large slid...
Definition: getfem_contact_and_friction_large_sliding.cc:2508
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:54
getfem::slice_vector_on_basic_dof_of_element
void slice_vector_on_basic_dof_of_element(const mesh_fem &mf, const VEC1 &vec, size_type cv, VEC2 &coeff, size_type qmult1=size_type(-1), size_type qmult2=size_type(-1))
Given a mesh_fem.
Definition: getfem_mesh_fem.h:659
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::virtual_brick
The virtual brick has to be derived to describe real model bricks.
Definition: getfem_models.h:1440
getfem::fem_interpolation_context::pf
const pfem pf() const
get the current FEM descriptor
Definition: getfem_fem.h:782
getfem::fem_interpolation_context::grad_base_value
void grad_base_value(base_tensor &t, bool withM=true) const
fill the tensor with the gradient of the base functions (taken at point this->xref())
Definition: getfem_fem.cc:202
getfem_contact_and_friction_common.h
Comomon tools for unilateral contact and Coulomb friction bricks.
gmm_condition_number.h
computation of the condition number of dense matrices.
getfem::fem_interpolation_context::base_value
void base_value(base_tensor &t, bool withM=true) const
fill the tensor with the values of the base functions (taken at point this->xref())
Definition: getfem_fem.cc:120
getfem::pbrick
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
Definition: getfem_models.h:49
getfem::transformation_name_of_large_sliding_contact_brick
const std::string & transformation_name_of_large_sliding_contact_brick(model &md, size_type indbrick)
Gives the name of the raytracing interpolate transformation for an existing large sliding contact bri...
Definition: getfem_contact_and_friction_large_sliding.cc:2488
gmm::mult_add
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1781
getfem_assembling.h
Miscelleanous assembly routines for common terms. Use the low-level generic assembly....
getfem::add_rigid_obstacle_to_raytracing_transformation
void add_rigid_obstacle_to_raytracing_transformation(model &md, const std::string &transname, const std::string &expr, size_type N)
Add a rigid obstacle whose geometry corresponds to the zero level-set of the high-level generic assem...
Definition: getfem_contact_and_friction_common.cc:2480