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