GetFEM  5.5
getfem_contact_and_friction_nodal.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-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 
22 #include "getfem/getfem_contact_and_friction_nodal.h"
23 #include "getfem/getfem_contact_and_friction_common.h"
24 #include "getfem/getfem_assembling.h"
25 
26 #ifndef GETFEM_HAVE_QHULL_QHULL_H
27 #include <getfem/bgeot_kdtree.h>
28 #endif
29 
30 namespace getfem {
31 
32  typedef bgeot::convex<base_node>::dref_convex_pt_ct dref_convex_pt_ct;
33  typedef bgeot::basic_mesh::ref_mesh_face_pt_ct ref_mesh_face_pt_ct;
34 
35 
36  // Computation of an orthonormal basis to a unit vector.
37  static void orthonormal_basis_to_unit_vec(size_type d, const base_node &un,
38  base_node *ut) {
39  size_type n = 0;
40  for (size_type k = 0; k <= d && n < d; ++k) {
41  gmm::resize(ut[n], d+1);
42  gmm::clear(ut[n]);
43  ut[n][k] = scalar_type(1);
44 
45  ut[n] -= gmm::vect_sp(un, ut[n]) * un;
46  for (size_type nn = 0; nn < n; ++nn)
47  ut[n] -= gmm::vect_sp(ut[nn], ut[n]) * ut[nn];
48 
49  if (gmm::vect_norm2(ut[n]) < 1e-3) continue;
50  ut[n] /= gmm::vect_norm2(ut[n]);
51  ++n;
52  }
53  GMM_ASSERT1(n == d, "Gram-Schmidt algorithm to find an "
54  "orthonormal basis for the tangential displacement failed");
55  }
56 
57  // "contact_node" is an object which contains data about nodes expected
58  // to participate in a contact condition. A contact node refers to a
59  // specific mesh_fem.
60  struct contact_node {
61  const mesh_fem *mf; // Pointer to the mesh_fem the contact node is
62  // associated with
63  size_type dof; // first dof id of the node in the considered mesh_fem
64  std::vector<size_type> cvs; // list of ids of neigbouring convexes
65  std::vector<short_type> fcs; // list of local ids of neigbouring faces
66 
67  contact_node() : mf(0), cvs(0), fcs(0) {}
68  contact_node(const mesh_fem &mf_) {mf = &mf_;}
69  };
70 
71  // contact_node's pair
72  struct contact_node_pair {
73  contact_node cn_s, cn_m; // Slave and master contact_node's
74  scalar_type dist2; // Square of distance between slave and master nodes
75  bool is_active;
76  contact_node_pair(scalar_type threshold=10.) : cn_s(), cn_m()
77  {dist2 = threshold * threshold; is_active = false;}
78  };
79 
80  // contact_node's pair list
81  class contact_node_pair_list : public std::vector<contact_node_pair> {
82 
83  void contact_node_list_from_region
84  (const mesh_fem &mf, size_type contact_region,
85  std::vector<contact_node> &cnl) {
86 
87  cnl.clear();
88  const mesh &m = mf.linked_mesh();
89  size_type qdim = mf.get_qdim();
90  std::map<size_type, size_type> dof_to_cnid;
91  size_type cnid = 0;
92  dal::bit_vector dofs = mf.basic_dof_on_region(contact_region);
93  for (dal::bv_visitor dof(dofs); !dof.finished(); ++dof)
94  if ( dof % qdim == 0) {
95  dof_to_cnid[dof] = cnid++;
96  contact_node new_cn(mf);
97  new_cn.dof = dof;
98  cnl.push_back(new_cn);
99  }
100  for (mr_visitor face(m.region(contact_region));
101  !face.finished(); ++face) {
102  assert(face.is_face());
103  mesh_fem::ind_dof_face_ct
104  face_dofs = mf.ind_basic_dof_of_face_of_element(face.cv(),face.f());
105  for (size_type it=0; it < face_dofs.size(); it += qdim ) {
106  size_type dof = face_dofs[it];
107  cnid = dof_to_cnid[dof];
108  cnl[cnid].cvs.push_back(face.cv());
109  cnl[cnid].fcs.push_back(face.f());
110  } // for:it
111  } // for:face
112  } // append
113 
114  public:
115  contact_node_pair_list() : std::vector<contact_node_pair>() {}
116 
117  void append_min_dist_cn_pairs(const mesh_fem &mf1, const mesh_fem &mf2,
118  size_type rg1, size_type rg2,
119  bool slave1=true, bool slave2=false) {
120 
121  std::vector<contact_node> cnl1(0), cnl2(0);
122  contact_node_list_from_region(mf1, rg1, cnl1);
123  contact_node_list_from_region(mf2, rg2, cnl2);
124 
125  // Find minimum distance node pairs
126  size_type size0 = this->size();
127  size_type size1 = slave1 ? cnl1.size() : 0;
128  size_type size2 = slave2 ? cnl2.size() : 0;
129  this->resize( size0 + size1 + size2 );
130 # ifndef GETFEM_HAVE_QHULL_QHULL_H
131  bgeot::kdtree tree1, tree2;
132  for (size_type i1 = 0; i1 < cnl1.size(); ++i1) {
133  contact_node *cn1 = &cnl1[i1];
134  tree1.add_point_with_id(cn1->mf->point_of_basic_dof(cn1->dof), i1);
135  }
136  for (size_type i2 = 0; i2 < cnl2.size(); ++i2) {
137  contact_node *cn2 = &cnl2[i2];
138  tree2.add_point_with_id(cn2->mf->point_of_basic_dof(cn2->dof), i2);
139  }
140  if (slave1) {
141  size_type ii1=size0;
142  for (size_type i1 = 0; i1 < cnl1.size(); ++i1, ++ii1) {
143  contact_node *cn1 = &cnl1[i1];
144  base_node node1 = cn1->mf->point_of_basic_dof(cn1->dof);
145  bgeot::index_node_pair ipt;
146  scalar_type dist2 = tree2.nearest_neighbor(ipt, node1);
147  if (ipt.i != size_type(-1) && dist2 < (*this)[ii1].dist2) {
148  (*this)[ii1].cn_s = *cn1;
149  (*this)[ii1].cn_m = cnl2[ipt.i];
150  (*this)[ii1].dist2 = dist2;
151  (*this)[ii1].is_active = true;
152  }
153  }
154  }
155  if (slave2) {
156  size_type ii2=size0+size1;
157  for (size_type i2 = 0; i2 < cnl2.size(); ++i2, ++ii2) {
158  contact_node *cn2 = &cnl2[i2];
159  base_node node2 = cn2->mf->point_of_basic_dof(cn2->dof);
160  bgeot::index_node_pair ipt;
161  scalar_type dist2 = tree1.nearest_neighbor(ipt, node2);
162  if (ipt.i != size_type(-1) && dist2 < (*this)[ii2].dist2) {
163  (*this)[ii2].cn_s = *cn2;
164  (*this)[ii2].cn_m = cnl1[ipt.i];
165  (*this)[ii2].dist2 = dist2;
166  (*this)[ii2].is_active = true;
167  }
168  }
169  }
170 # else
171  std::vector<base_node> pts;
172  for (size_type i1 = 0; i1 < cnl1.size(); ++i1) {
173  contact_node *cn1 = &cnl1[i1];
174  pts.push_back(cn1->mf->point_of_basic_dof(cn1->dof));
175  }
176  for (size_type i2 = 0; i2 < cnl2.size(); ++i2) {
177  contact_node *cn2 = &cnl2[i2];
178  pts.push_back(cn2->mf->point_of_basic_dof(cn2->dof));
179  }
180  gmm::dense_matrix<size_type> simplexes;
181 
182  getfem::delaunay(pts, simplexes);
183 
184  size_type nb_vertices = gmm::mat_nrows(simplexes);
185  std::vector<size_type> facet_vertices(nb_vertices);
186  std::vector< std::vector<size_type> > pt1_neighbors(size1);
187  for (size_type i = 0; i < gmm::mat_ncols(simplexes); ++i) {
188  gmm::copy(gmm::mat_col(simplexes, i), facet_vertices);
189  for (size_type iv1 = 0; iv1 < nb_vertices-1; ++iv1) {
190  size_type v1 = facet_vertices[iv1];
191  bool v1_on_surface1 = (v1 < size1);
192  for (size_type iv2 = iv1 + 1; iv2 < nb_vertices; ++iv2) {
193  size_type v2 = facet_vertices[iv2];
194  bool v2_on_surface1 = (v2 < size1);
195  if (v1_on_surface1 ^ v2_on_surface1) {
196  bool already_in = false;
197  size_type vv1 = (v1_on_surface1 ? v1 : v2);
198  size_type vv2 = (v2_on_surface1 ? v1 : v2);
199  for (size_type j = 0; j < pt1_neighbors[vv1].size(); ++j)
200  if (pt1_neighbors[vv1][j] == vv2) {
201  already_in = true;
202  break;
203  }
204  if (!already_in) pt1_neighbors[vv1].push_back(vv2);
205  }
206  }
207  }
208  }
209 
210  for (size_type i1 = 0; i1 < size1; ++i1)
211  for (size_type j = 0; j < pt1_neighbors[i1].size(); ++j) {
212  size_type i2 = pt1_neighbors[i1][j] - size1;
213  size_type ii1 = size0 + i1;
214  size_type ii2 = size0 + size1 + i2;
215  contact_node *cn1 = &cnl1[i1];
216  base_node node1 = cn1->mf->point_of_basic_dof(cn1->dof);
217  contact_node *cn2 = &cnl2[i2];
218  base_node node2 = cn2->mf->point_of_basic_dof(cn2->dof);
219  scalar_type dist2 = gmm::vect_norm2_sqr(node1-node2);
220  if (slave1 && dist2 < (*this)[ii1].dist2) {
221  (*this)[ii1].cn_s = *cn1;
222  (*this)[ii1].cn_m = *cn2;
223  (*this)[ii1].dist2 = dist2;
224  (*this)[ii1].is_active = true;
225  }
226  if (slave2 && dist2 < (*this)[ii2].dist2) {
227  (*this)[ii2].cn_s = *cn2;
228  (*this)[ii2].cn_m = *cn1;
229  (*this)[ii2].dist2 = dist2;
230  (*this)[ii2].is_active = true;
231  }
232  }
233 #endif
234  }
235 
236  void append_min_dist_cn_pairs(const mesh_fem &mf,
237  size_type rg1, size_type rg2,
238  bool slave1=true, bool slave2=false) {
239  append_min_dist_cn_pairs(mf, mf, rg1, rg2, slave1, slave2);
240  }
241  };
242 
243  scalar_type projection_on_convex_face
244  (const mesh &m, const size_type cv, const short_type fc,
245  const base_node &master_node, const base_node &slave_node,
246  base_node &un, base_node &proj_node, base_node &proj_node_ref) {
247 
248  bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
249 
250  if (pgt->is_linear()) { //this condition is practically too strict
251 
252  un = m.normal_of_face_of_convex(cv,fc);
253  un /= gmm::vect_norm2(un);
254  //proj_node = slave_node - [(slave_node-master_node)*n] * n
255  gmm::add(master_node, gmm::scaled(slave_node, -1.), proj_node);
256  gmm::copy(gmm::scaled(un, gmm::vect_sp(proj_node, un)), proj_node);
257  gmm::add(slave_node, proj_node);
258 
259  bgeot::geotrans_inv_convex gic;
260  gic.init(m.points_of_convex(cv), pgt);
261  gic.invert(proj_node, proj_node_ref);
262  return pgt->convex_ref()->is_in(proj_node_ref);
263 
264  } else {
265 
266  size_type N = m.dim();
267  size_type P = pgt->structure()->dim();
268 
269  size_type nb_pts_cv = pgt->nb_points();
270  size_type nb_pts_fc = pgt->structure()->nb_points_of_face(fc);
271 
272  bgeot::convex_ind_ct ind_pts_fc = pgt->structure()->ind_points_of_face(fc);
273  ref_mesh_face_pt_ct pts_fc = m.points_of_face_of_convex(cv, fc);
274 
275  // Local base on reference face
276  base_matrix base_ref_fc(P-1,N);
277  {
278  dref_convex_pt_ct dref_pts_fc = pgt->convex_ref()->dir_points_of_face(fc);
279  GMM_ASSERT1( dref_pts_fc.size() == P, "Dimensions mismatch");
280  base_node vec(dref_pts_fc[0].size());
281  for (size_type i = 0; i < P-1; ++i) {
282  vec = dref_pts_fc[i+1] - dref_pts_fc[0];
283  gmm::copy(vec,gmm::mat_row(base_ref_fc,i));
284  }
285  }
286 
287  GMM_ASSERT1( slave_node.size() == N, "Dimensions mismatch");
288  const base_node &xx = slave_node;
289  base_node &xxp = proj_node; xxp.resize(N);
290  base_node &xp = proj_node_ref; xp.resize(P);
291  base_node vres(P);
292  scalar_type res= 1.;
293 
294  // initial guess
295  xp = gmm::mean_value(pgt->convex_ref()->points_of_face(fc));
296 
297  gmm::clear(xxp);
298  base_vector val(nb_pts_fc);
299  pgt->poly_vector_val(xp, ind_pts_fc, val);
300  for (size_type l = 0; l < nb_pts_fc; ++l)
301  gmm::add(gmm::scaled(pts_fc[l], val[l] ), xxp);
302 
303  base_matrix G(N, nb_pts_fc);
304  vectors_to_base_matrix(G, pts_fc);
305 
306  base_matrix K(N,P-1);
307 
308  base_matrix grad_fc(nb_pts_fc, P);
309  base_matrix grad_fc1(nb_pts_fc, P-1);
310  base_matrix B(N,P-1), BB(N,P), CS(P-1,P-1);
311 
312  scalar_type EPS = 10E-12;
313  unsigned cnt = 50;
314  while (res > EPS && --cnt) {
315  // computation of the pseudo inverse matrix B at point xp
316  pgt->poly_vector_grad(xp, ind_pts_fc, grad_fc);
317  gmm::mult(grad_fc, gmm::transposed(base_ref_fc), grad_fc1);
318  gmm::mult(G, grad_fc1, K);
319  gmm::mult(gmm::transposed(K), K, CS);
320  bgeot::lu_inverse(&(*(CS.begin())), P-1);
321  gmm::mult(K, CS, B);
322  gmm::mult(B, base_ref_fc, BB);
323 
324  // Projection onto the face of the convex
325  gmm::mult_add(gmm::transposed(BB), xx-xxp, xp);
326  gmm::clear(xxp);
327  pgt->poly_vector_val(xp, ind_pts_fc, val);
328  for (size_type l = 0; l < nb_pts_fc; ++l)
329  gmm::add(gmm::scaled(pts_fc[l], val[l]), xxp);
330 
331  gmm::mult(gmm::transposed(BB), xx - xxp, vres);
332  res = gmm::vect_norm2(vres);
333  }
334  GMM_ASSERT1( res <= EPS,
335  "Iterative pojection on convex face did not converge");
336  { // calculate K at the final point
337  pgt->poly_vector_grad(xp, ind_pts_fc, grad_fc);
338  gmm::mult(grad_fc, gmm::transposed(base_ref_fc), grad_fc1);
339  gmm::mult(G, grad_fc1, K);
340  }
341 
342  // computation of normal vector
343  un.resize(N);
344  // un = xx - xxp;
345  // gmm::scale(un, 1/gmm::vect_norm2(un));
346  gmm::clear(un);
347  {
348  base_matrix KK(N,P);
349  { // calculate KK
350  base_matrix grad_cv(nb_pts_cv, P);
351  pgt->poly_vector_grad(xp, grad_cv);
352 
353  base_matrix GG(N, nb_pts_cv);
354  m.points_of_convex(cv, GG);
355 
356  gmm::mult(GG, grad_cv, KK);
357  }
358 
359  base_matrix bases_product(P-1, P);
360  gmm::mult(gmm::transposed(K), KK, bases_product);
361 
362  for (size_type i = 0; i < P; ++i) {
363  std::vector<size_type> ind(0);
364  for (size_type j = 0; j < P; ++j)
365  if (j != i ) ind.push_back(j);
366  gmm::sub_index SUBI(ind);
367  scalar_type det
368  = gmm::lu_det(gmm::sub_matrix(bases_product,
369  gmm::sub_interval(0, P-1), SUBI));
370  gmm::add(gmm::scaled(gmm::mat_col(KK, i), (i % 2) ? -det : +det ),
371  un);
372  }
373  }
374  // normalizing
375  gmm::scale(un, 1/gmm::vect_norm2(un));
376  // ensure that normal points outwards
377  if (gmm::vect_sp(un, gmm::mean_value(pts_fc) -
378  gmm::mean_value(m.points_of_convex(cv))) < 0)
379  gmm::scale(un,scalar_type(-1));
380 
381  return pgt->convex_ref()->is_in(proj_node_ref);
382  }
383  }
384 
385  void compute_contact_matrices
386  (const mesh_fem &mf_disp1, const mesh_fem &mf_disp2,
387  contact_node_pair_list &cnpl, model_real_plain_vector &gap,
388  CONTACT_B_MATRIX *BN1, CONTACT_B_MATRIX *BN2 = 0,
389  CONTACT_B_MATRIX *BT1 = 0, CONTACT_B_MATRIX *BT2 = 0) {
390 
391  GMM_ASSERT1(gmm::vect_size(gap) == cnpl.size(),
392  "Wrong number of contact node pairs or wrong size of gap");
393  gmm::clear(*BN1);
394  GMM_ASSERT1( gmm::mat_nrows(*BN1) == cnpl.size(), "Wrong size of BN1");
395  if (BN2) {
396  gmm::clear(*BN2);
397  GMM_ASSERT1( gmm::mat_nrows(*BN2) == cnpl.size(), "Wrong size of BN2");
398  }
399  dim_type qdim = mf_disp1.get_qdim();
400  size_type d = qdim - 1;
401  if (BT1) {
402  gmm::clear(*BT1);
403  GMM_ASSERT1( gmm::mat_nrows(*BT1) == cnpl.size() * d, "Wrong size of BT1");
404  }
405  if (BT2) {
406  gmm::clear(*BT2);
407  GMM_ASSERT1( gmm::mat_nrows(*BT2) == cnpl.size() * d, "Wrong size of BT2");
408  }
409  gmm::fill(gap, scalar_type(10)); //FIXME: Needs a threshold value
410  for (size_type row = 0; row < cnpl.size(); ++row) {
411  contact_node_pair *cnp = &cnpl[row];
412  if (cnp->is_active) {
413  contact_node *cn_s = &cnp->cn_s; //slave contact node
414  contact_node *cn_m = &cnp->cn_m; //master contact node
415  const mesh &mesh_m = cn_m->mf->linked_mesh();
416  base_node slave_node = cn_s->mf->point_of_basic_dof(cn_s->dof);
417  base_node master_node = cn_m->mf->point_of_basic_dof(cn_m->dof);
418  GMM_ASSERT1(slave_node.size() == qdim && master_node.size() == qdim,
419  "Internal error");
420  base_node un_sel(qdim), proj_node_sel(qdim), proj_node_ref_sel(qdim);
421  scalar_type is_in_min = 1e5; //FIXME
422  size_type cv_sel = 0;
423  short_type fc_sel = 0;
424  std::vector<size_type>::iterator cv;
425  std::vector<short_type>::iterator fc;
426  for (cv = cn_m->cvs.begin(), fc = cn_m->fcs.begin();
427  cv != cn_m->cvs.end() && fc != cn_m->fcs.end(); cv++, fc++) {
428  base_node un(qdim), proj_node(qdim), proj_node_ref(qdim);
429  scalar_type is_in = projection_on_convex_face
430  (mesh_m, *cv, *fc, master_node, slave_node, un, proj_node, proj_node_ref);
431  if (is_in < is_in_min) {
432  is_in_min = is_in;
433  cv_sel = *cv;
434  fc_sel = *fc;
435  un_sel = un;
436  proj_node_sel = proj_node;
437  proj_node_ref_sel = proj_node_ref;
438  }
439  }
440  if (is_in_min < 0.05) { //FIXME
441  gap[row] = gmm::vect_sp(slave_node-proj_node_sel, un_sel);
442 
443  std::vector<base_node> ut(d);
444  if (BT1) orthonormal_basis_to_unit_vec(d, un_sel, &(ut[0]));
445 
446  CONTACT_B_MATRIX *BN = 0;
447  CONTACT_B_MATRIX *BT = 0;
448  if (cn_s->mf == &mf_disp1) {
449  BN = BN1;
450  BT = BT1;
451  } else if (cn_s->mf == &mf_disp2) {
452  BN = BN2;
453  BT = BT2;
454  }
455  if (BN)
456  for (size_type k = 0; k <= d; ++k)
457  (*BN)(row, cn_s->dof + k) -= un_sel[k];
458  if (BT)
459  for (size_type k = 0; k <= d; ++k)
460  for (size_type n = 0; n < d; ++n)
461  (*BT)(row * d + n, cn_s->dof + k) -= ut[n][k];
462 
463  BN = 0;
464  const mesh_fem *mf_disp = 0;
465  if (cn_m->mf == &mf_disp1) {
466  BN = BN1;
467  BT = BT1;
468  mf_disp = &mf_disp1;
469  } else if (cn_m->mf == &mf_disp2) {
470  BN = BN2;
471  BT = BT2;
472  mf_disp = &mf_disp2;
473  }
474  if (BN) {
475  base_matrix G;
476  base_matrix M(qdim, mf_disp->nb_basic_dof_of_element(cv_sel));
477  mesh_m.points_of_convex(cv_sel, G);
478  pfem pf = mf_disp->fem_of_element(cv_sel);
479  bgeot::pgeometric_trans pgt = mesh_m.trans_of_convex(cv_sel);
480  fem_interpolation_context
481  ctx(pgt, pf, proj_node_ref_sel, G, cv_sel, fc_sel);
482  pf->interpolation (ctx, M, int(qdim));
483 
484  mesh_fem::ind_dof_ct
485  master_dofs = mf_disp->ind_basic_dof_of_element(cv_sel);
486 
487  model_real_plain_vector MT_u(mf_disp->nb_basic_dof_of_element(cv_sel));
488  gmm::mult(gmm::transposed(M), un_sel, MT_u);
489  for (size_type j = 0; j < master_dofs.size(); ++j)
490  (*BN)(row, master_dofs[j]) += MT_u[j];
491 
492  if (BT) {
493  for (size_type n = 0; n < d; ++n) {
494  gmm::mult(gmm::transposed(M), ut[n], MT_u);
495  for (size_type j = 0; j < master_dofs.size(); ++j)
496  (*BT)(row * d + n, master_dofs[j]) += MT_u[j];
497  }
498  }
499  } // BN
500 
501  }
502  } // if:cnp->cn_s
503  } // cnp
504 
505  } // compute_contact_matrices
506 
507 
508 
509  //=========================================================================
510  //
511  // Basic Brick (with given BN, BT, gap) and possibly two bodies
512  //
513  //=========================================================================
514 
515  struct Coulomb_friction_brick : public virtual_brick {
516 
517  mutable CONTACT_B_MATRIX BN1, BT1, BN2, BT2;
518  mutable CONTACT_B_MATRIX DN, DDN, DT, DDT; // For Hughes stabilization
519  mutable CONTACT_B_MATRIX BBN1, BBT1, BBN2, BBT2;
520  mutable model_real_plain_vector gap, threshold, friction_coeff, alpha;
521  mutable model_real_plain_vector RLN, RLT;
522  mutable scalar_type r, gamma;
523  mutable bool is_init;
524  bool Tresca_version, contact_only;
525  bool really_stationary, friction_dynamic_term;
526  bool two_variables, Hughes_stabilized;
527  int augmentation_version; // 0 for non-symmetric Alart-Curnier version
528  // 1 for symmetric Alart-Curnier version
529  // 2 for new version (augmented multipliers)
530  // 3 for new version with De Saxcé projection
531 
532  void init_BBN_BBT(void) const {
533  gmm::resize(BBN1, gmm::mat_nrows(BN1), gmm::mat_ncols(BN1));
534  gmm::copy(BN1, BBN1);
535  if (Hughes_stabilized) {
536  gmm::resize(DDN, gmm::mat_nrows(DN), gmm::mat_ncols(DN));
537  gmm::copy(DN, DDN);
538  }
539  if (two_variables) {
540  gmm::resize(BBN2, gmm::mat_nrows(BN2), gmm::mat_ncols(BN2));
541  gmm::copy(BN2, BBN2);
542  }
543  if (!contact_only) {
544  if (Hughes_stabilized) {
545  gmm::resize(DDT, gmm::mat_nrows(DT), gmm::mat_ncols(DT));
546  gmm::copy(DT, DDT);
547  }
548  gmm::resize(BBT1, gmm::mat_nrows(BT1), gmm::mat_ncols(BT1));
549  gmm::copy(BT1, BBT1);
550  if (two_variables) {
551  gmm::resize(BBT2, gmm::mat_nrows(BT2), gmm::mat_ncols(BT2));
552  gmm::copy(BT2, BBT2);
553  }
554  }
555  size_type nbc = gmm::mat_nrows(BN1);
556  size_type d = gmm::mat_nrows(BT1)/nbc;
557  for (size_type i = 0; i < nbc; ++i) {
558  gmm::scale(gmm::mat_row(BBN1, i), alpha[i]);
559  if (Hughes_stabilized) gmm::scale(gmm::mat_row(DDN, i), alpha[i]);
560  if (two_variables)
561  gmm::scale(gmm::mat_row(BBN2, i), alpha[i]);
562  if (!contact_only)
563  for (size_type k = 0; k < d; ++k) {
564  if (Hughes_stabilized)
565  gmm::scale(gmm::mat_row(DDT, d*i+k), alpha[i]);
566  gmm::scale(gmm::mat_row(BBT1, d*i+k), alpha[i]);
567  if (two_variables)
568  gmm::scale(gmm::mat_row(BBT2, d*i+k), alpha[i]);
569  }
570  }
571  is_init = true;
572  }
573 
574  void set_stationary(void) { really_stationary = true; }
575 
576  void precomp(const model_real_plain_vector &u1,
577  const model_real_plain_vector &u2,
578  const model_real_plain_vector &lambda_n,
579  const model_real_plain_vector &lambda_t,
580  const model_real_plain_vector &wt1,
581  const model_real_plain_vector &wt2) const {
582  gmm::copy(gmm::scaled(gap, r), RLN);
583  for (size_type i = 0; i < gmm::mat_nrows(BN1); ++i) RLN[i] *= alpha[i];
584  gmm::add(lambda_n, RLN);
585  gmm::mult_add(BBN1, gmm::scaled(u1, -r), RLN);
586  if (Hughes_stabilized)
587  gmm::mult_add(DDN, gmm::scaled(lambda_n, -r), RLN);
588  if (two_variables) gmm::mult_add(BBN2, gmm::scaled(u2, -r), RLN);
589  if (!contact_only) {
590  gmm::copy(lambda_t, RLT);
591  if (friction_dynamic_term) {
592  gmm::mult_add(BBT1, gmm::scaled(wt1, -r*gamma), RLT);
593  if (two_variables)
594  gmm::mult_add(BBT2, gmm::scaled(wt2, -r*gamma), RLT);
595  }
596  if (!really_stationary) {
597  gmm::mult_add(BBT1, gmm::scaled(u1, -r), RLT);
598  if (two_variables) gmm::mult_add(BBT2, gmm::scaled(u2, -r), RLT);
599  }
600  if (Hughes_stabilized)
601  gmm::mult_add(DDT, gmm::scaled(lambda_t, -r), RLT);
602  }
603  }
604 
605  // Common part for all contact with friction bricks
606  void basic_asm_real_tangent_terms(const model_real_plain_vector &u1,
607  const model_real_plain_vector &u2,
608  const model_real_plain_vector &lambda_n,
609  const model_real_plain_vector &lambda_t,
610  const model_real_plain_vector &wt1,
611  const model_real_plain_vector &wt2,
612  model::real_matlist &matl,
613  model::real_veclist &vecl,
614  build_version version) const {
615  size_type nbt = 4 + (contact_only ? 0 : 4) + (two_variables ? 3 : 0)
616  + (two_variables && !contact_only ? 2 : 0);
617  GMM_ASSERT1(matl.size() == nbt,
618  "Wrong number of terms for the contact brick");
619 
620  const scalar_type vt1 = scalar_type(1), vt0 = scalar_type(0);
621  size_type nbc = gmm::mat_nrows(BN1);
622  size_type d = gmm::mat_nrows(BT1)/nbc;
623 
624  // Matrices to be filled
625  size_type nt = 0;
626  model_real_sparse_matrix &T_u1_u1 = matl[nt++], &T_u2_u2 = matl[nt++];
627  if (!two_variables) nt--;
628  model_real_sparse_matrix &T_u1_n = matl[nt++], &T_n_u1 = matl[nt++];
629  if (!two_variables) nt -= 2;
630  model_real_sparse_matrix &T_u2_n = matl[nt++], &T_n_u2 = matl[nt++];
631  size_type nvec_lambda_n = nt;
632  model_real_sparse_matrix &T_n_n = matl[nt++];
633  if (contact_only) nt -= 2;
634  model_real_sparse_matrix &T_u1_t = matl[nt++], &T_t_u1 = matl[nt++];
635  if (contact_only || !two_variables) nt -= 2;
636  model_real_sparse_matrix &T_u2_t = matl[nt++], &T_t_u2 = matl[nt++];
637  if (contact_only) nt -= 2;
638  size_type nvec_lambda_t = nt;
639  model_real_sparse_matrix &T_t_t = matl[nt++], &T_t_n = matl[nt++];
640 
641  // Rhs to be filled
642  model_real_plain_vector &ru1 = vecl[0];
643  model_real_plain_vector &ru2 = vecl[1];
644  model_real_plain_vector &rlambda_n = vecl[nvec_lambda_n];
645  model_real_plain_vector &rlambda_t = vecl[nvec_lambda_t];
646 
647  // pre-computations
648  if (!is_init) init_BBN_BBT();
649  gmm::resize(RLN, nbc);
650  if (!contact_only) gmm::resize(RLT, nbc*d);
651  if (augmentation_version <= 2)
652  precomp(u1, u2, lambda_n, lambda_t, wt1, wt2);
653 
654  if (version & model::BUILD_MATRIX) {
655  base_matrix pg(d, d);
656  base_vector vg(d);
657 
658  gmm::clear(T_n_n); gmm::clear(T_n_u1);
659  gmm::clear(T_u1_n); gmm::clear(T_u1_u1);
660  if (two_variables)
661  { gmm::clear(T_u2_u2); gmm::clear(T_n_u2); gmm::clear(T_u2_n); }
662  if (!contact_only) {
663  gmm::clear(T_u1_t); gmm::clear(T_t_n); gmm::clear(T_t_t);
664  if (two_variables) gmm::clear(T_u2_t);
665  }
666 
667  switch (augmentation_version) {
668  case 1: case 2:
669  gmm::copy(gmm::scaled(gmm::transposed(BN1), -vt1), T_u1_n);
670  if (two_variables)
671  gmm::copy(gmm::scaled(gmm::transposed(BN2), -vt1), T_u2_n);
672  for (size_type i=0; i < nbc; ++i) {
673  if (RLN[i] > vt0) {
674  gmm::clear(gmm::mat_col(T_u1_n, i));
675  if (two_variables) gmm::clear(gmm::mat_col(T_u2_n, i));
676  T_n_n(i, i) = -vt1/(r*alpha[i]);
677  }
678  if (Hughes_stabilized && RLN[i] <= vt0)
679  gmm::copy(gmm::scaled(gmm::mat_row(DN, i), -vt1),
680  gmm::mat_col(T_n_n, i));
681  }
682  if (Hughes_stabilized) {
683  model_real_sparse_matrix aux(nbc, nbc);
684  gmm::copy(gmm::transposed(T_n_n), aux);
685  gmm::copy(aux, T_n_n);
686  }
687  gmm::copy(gmm::transposed(T_u1_n), T_n_u1);
688  if (two_variables) gmm::copy(gmm::transposed(T_u2_n), T_n_u2);
689  if (!contact_only) {
690  for (size_type i=0; i < nbc; ++i) {
691  gmm::sub_interval SUBI(i*d, d);
692  scalar_type th = Tresca_version ? threshold[i]
693  : - (std::min(vt0, RLN[i])) * friction_coeff[i];
694  ball_projection_grad(gmm::sub_vector(RLT, SUBI), th, pg);
695  if (!really_stationary)
696  for (size_type k1 = 0; k1 < d; ++k1)
697  for (size_type k2 = 0; k2 < d; ++k2) {
698  gmm::add(gmm::scaled(gmm::mat_row(BT1,i*d+k1),-pg(k2,k1)),
699  gmm::mat_col(T_u1_t, i*d+k2));
700  if (two_variables)
701  gmm::add(gmm::scaled(gmm::mat_row(BT2,i*d+k1),
702  -pg(k2,k1)),
703  gmm::mat_col(T_u2_t, i*d+k2));
704  }
705 
706  if (!Tresca_version) {
707  if (RLN[i] <= vt0) {
708  ball_projection_grad_r(gmm::sub_vector(RLT, SUBI), th, vg);
709  for (size_type k = 0; k < d; ++k)
710  T_t_n(i*d+k, i) = - friction_coeff[i] * vg[k]/(r*alpha[i]);
711  }
712  }
713  for (size_type k = 0; k < d; ++k) pg(k,k) -= vt1;
714  gmm::copy(gmm::scaled(pg, vt1/(r*alpha[i])),
715  gmm::sub_matrix(T_t_t,SUBI));
716  if (Hughes_stabilized) {
717  for (size_type k = 0; k < d; ++k)
718  for (size_type l = 0; l < d; ++l) {
719  gmm::add(gmm::scaled(gmm::mat_row(DT, d*i+l), -pg(k,l)),
720  gmm::mat_col(T_t_t, d*i+k));
721  }
722  }
723 
724  }
725  if (Hughes_stabilized) {
726  model_real_sparse_matrix aux(gmm::mat_nrows(T_t_t),
727  gmm::mat_nrows(T_t_t));
728  gmm::copy(gmm::transposed(T_t_t), aux);
729  gmm::copy(aux, T_t_t);
730  }
731  gmm::copy(gmm::transposed(T_u1_t), T_t_u1);
732  if (two_variables) gmm::copy(gmm::transposed(T_u2_t), T_t_u2);
733  }
734 
735  if (augmentation_version == 1) {
736  gmm::copy(gmm::scaled(gmm::transposed(BN1), -vt1), T_u1_n);
737  if (two_variables)
738  gmm::copy(gmm::scaled(gmm::transposed(BN2), -vt1), T_u2_n);
739  if (!contact_only) {
740  gmm::copy(gmm::scaled(gmm::transposed(BT1), -vt1), T_u1_t);
741  if (two_variables)
742  gmm::copy(gmm::scaled(gmm::transposed(BT2), -vt1), T_u2_t);
743  }
744  } else {
745  model_real_sparse_matrix tmp1(gmm::mat_ncols(BN1),
746  gmm::mat_ncols(BN1));
747  gmm::mult(gmm::transposed(gmm::scaled(BBN1,-r)), T_n_u1, tmp1);
748  gmm::add(tmp1, T_u1_u1);
749  if (two_variables) {
750  gmm::mult(gmm::transposed(gmm::scaled(BBN2,-r)), T_n_u2, tmp1);
751  gmm::add(tmp1, T_u2_u2);
752  }
753 
754  if (!contact_only) {
755  gmm::mult(gmm::transposed(gmm::scaled(BBT1,-r)), T_t_u1, tmp1);
756  gmm::add(tmp1, T_u1_u1);
757  if (two_variables) {
758  gmm::mult(gmm::transposed(gmm::scaled(BBT2,-r)), T_t_u2, tmp1);
759  gmm::add(tmp1, T_u2_u2);
760  }
761  }
762  }
763 
764  if (!contact_only && !Tresca_version) {
765  // should be simplified ... !
766  model_real_sparse_matrix tmp5(gmm::mat_ncols(BT1),
767  gmm::mat_ncols(BT1));
768  model_real_sparse_matrix tmp6(gmm::mat_ncols(BT1),
769  gmm::mat_ncols(BT1));
770  model_real_sparse_matrix tmp7(gmm::mat_ncols(BT2),
771  gmm::mat_ncols(BT2));
772  model_real_sparse_matrix tmp8(gmm::mat_ncols(BT2),
773  gmm::mat_ncols(BT2));
774  model_real_sparse_matrix tmp3(gmm::mat_ncols(T_t_u1),
775  gmm::mat_nrows(T_t_u1));
776  model_real_sparse_matrix tmp4(gmm::mat_ncols(T_t_u2),
777  gmm::mat_nrows(T_t_u2));
778 
779  for (size_type i=0; i < nbc; ++i) {
780  gmm::sub_interval SUBI(i*d, d);
781  scalar_type th = - (std::min(vt0, RLN[i])) * friction_coeff[i];
782  if (RLN[i] <= vt0) {
783  ball_projection_grad_r(gmm::sub_vector(RLT, SUBI), th, vg);
784  for (size_type k = 0; k < d; ++k) {
785  gmm::add(gmm::scaled(gmm::mat_row(BN1,i),
786  vg[k]*friction_coeff[i]), gmm::mat_col(tmp3, i*d+k));
787  if (two_variables)
788  gmm::add(gmm::scaled(gmm::mat_row(BN2,i),
789  vg[k]*friction_coeff[i]), gmm::mat_col(tmp4, i*d+k));
790 
791  if (augmentation_version == 2) {
792 
793  gmm::add(gmm::scaled(gmm::mat_row(BT1,i*d+k),
794  vg[k]*friction_coeff[i]), gmm::mat_col(T_u1_n, i));
795  if (two_variables)
796  gmm::add(gmm::scaled(gmm::mat_row(BT2,i*d+k),
797  vg[k]*friction_coeff[i]), gmm::mat_col(T_u2_n, i));
798 
799  gmm::copy(gmm::scaled(gmm::mat_row(BBT1,i*d+k),
800  -r*friction_coeff[i]*vg[k]),
801  gmm::mat_col(tmp5, i*d+k));
802  gmm::copy(gmm::mat_row(BN1,i),
803  gmm::mat_col(tmp6, i*d+k));
804  if (two_variables) {
805  gmm::copy(gmm::scaled(gmm::mat_row(BBT2,i*d+k),
806  -r*friction_coeff[i]*vg[k]),
807  gmm::mat_col(tmp7, i*d+k));
808  gmm::copy(gmm::mat_row(BN2,i),
809  gmm::mat_col(tmp8, i*d+k));
810  }
811  }
812  }
813  }
814  }
815 
816  gmm::add(gmm::transposed(tmp3), T_t_u1);
817  if (two_variables)
818  gmm::add(gmm::transposed(tmp4), T_t_u2);
819 
820  if (augmentation_version == 2) {
821  model_real_sparse_matrix tmp1(gmm::mat_ncols(BN1),
822  gmm::mat_ncols(BN1));
823  gmm::mult(tmp5, gmm::transposed(tmp6), gmm::transposed(tmp1));
824  gmm::add(gmm::transposed(tmp1), T_u1_u1);
825  if (two_variables) {
826  gmm::mult(tmp7, gmm::transposed(tmp8),gmm::transposed(tmp1));
827  gmm::add(gmm::transposed(tmp1), T_u2_u2);
828  }
829  }
830  }
831  break;
832 
833  case 3:
834  gmm::copy(gmm::scaled(gmm::transposed(BN1), -vt1), T_u1_n);
835  if (two_variables)
836  gmm::copy(gmm::scaled(gmm::transposed(BN2), -vt1), T_u2_n);
837  for (size_type i=0; i < nbc; ++i) {
838  if (lambda_n[i] > vt0) {
839  gmm::clear(gmm::mat_col(T_u1_n, i));
840  if (two_variables) gmm::clear(gmm::mat_col(T_u2_n, i));
841  T_n_n(i, i) = -vt1/r;
842  }
843  }
844  gmm::copy(gmm::scaled(BN1, -vt1), T_n_u1);
845  if (two_variables) gmm::copy(gmm::scaled(BN2, -r), T_n_u2);
846  if (!contact_only) {
847  for (size_type i=0; i < nbc; ++i) {
848  gmm::sub_interval SUBI(i*d, d);
849  scalar_type th = Tresca_version ? threshold[i]
850  : gmm::neg(lambda_n[i]) * friction_coeff[i];
851  ball_projection_grad(gmm::sub_vector(lambda_t, SUBI), th, pg);
852  if (!really_stationary)
853  for (size_type k1 = 0; k1 < d; ++k1)
854  for (size_type k2 = 0; k2 < d; ++k2) {
855  gmm::add(gmm::scaled(gmm::mat_row(BT1,i*d+k1),-pg(k2,k1)),
856  gmm::mat_col(T_u1_t, i*d+k2));
857  if (two_variables)
858  gmm::add(gmm::scaled(gmm::mat_row(BT2,i*d+k1),
859  -pg(k2,k1)),
860  gmm::mat_col(T_u2_t, i*d+k2));
861  }
862  if (!Tresca_version) {
863  ball_projection_grad_r(gmm::sub_vector(lambda_t, SUBI),th,vg);
864  for (size_type k1 = 0; k1 < d; ++k1) {
865  gmm::add(gmm::scaled(gmm::mat_row(BT1,i*d+k1),
866  friction_coeff[i]*vg[k1]),
867  gmm::mat_col(T_u1_n, i));
868  if (two_variables)
869  gmm::add(gmm::scaled(gmm::mat_row(BT2,i*d+k1),
870  friction_coeff[i]*vg[k1]),
871  gmm::mat_col(T_u2_n, i));
872  T_t_n(i*d+k1, i) = friction_coeff[i] * vg[k1] / (r*alpha[i]);
873  }
874  }
875  for (size_type k = 0; k < d; ++k) pg(k,k) -= vt1;
876 
877  gmm::copy(gmm::scaled(pg, vt1/(r*alpha[i])),
878  gmm::sub_matrix(T_t_t, SUBI));
879 
880  }
881  gmm::copy(gmm::scaled(BT1, -vt1), T_t_u1);
882  if (two_variables) gmm::copy(gmm::scaled(BT2, -r), T_t_u2);
883  }
884  break;
885 
886  case 4: // Desaxce projection
887  base_small_vector x(d+1), n(d+1), u(d); n[0] = vt1;
888  base_matrix g(d+1, d+1);
889  model_real_sparse_matrix T_n_u1_transp(gmm::mat_ncols(T_n_u1), nbc);
890  model_real_sparse_matrix T_n_u2_transp(gmm::mat_ncols(T_n_u2), nbc);
891 
892  gmm::mult(BT1, u1, RLT);
893  if (two_variables) gmm::mult_add(BT2, u2, RLT);
894 
895  for (size_type i=0; i < nbc; ++i) {
896  x[0] = lambda_n[i];
897  for (size_type j=0; j < d; ++j) x[1+j] = lambda_t[i*d+j];
898  De_Saxce_projection_grad(x, n, friction_coeff[i], g);
899 
900  gmm::add(gmm::scaled(gmm::mat_row(BN1, i), -g(0,0)),
901  gmm::mat_col(T_u1_n, i));
902  if (two_variables)
903  gmm::add(gmm::scaled(gmm::mat_row(BN2, i), -g(0,0)),
904  gmm::mat_col(T_u2_n, i));
905  T_n_n(i, i) = (g(0,0) - vt1)/(r*alpha[i]);
906 
907  gmm::copy(gmm::sub_vector(RLT, gmm::sub_interval(i*d,d)), u);
908  scalar_type nu = gmm::vect_norm2(u);
909  if (nu != vt0)
910  for (size_type j=0; j < d; ++j) {
911  gmm::add(gmm::scaled(gmm::mat_row(BT1, i*d+j),
912  friction_coeff[i] * u[j] / nu),
913  gmm::mat_col(T_n_u1_transp, i));
914  if (two_variables)
915  gmm::add(gmm::scaled(gmm::mat_row(BT2, i*d+j),
916  friction_coeff[i] * u[j] / nu),
917  gmm::mat_col(T_n_u2_transp, i));
918  }
919 
920  for (size_type j=0; j < d; ++j) {
921  gmm::add(gmm::scaled(gmm::mat_row(BT1, i*d+j), -g(0,j+1)),
922  gmm::mat_col(T_u1_n, i));
923  if (two_variables)
924  gmm::add(gmm::scaled(gmm::mat_row(BT2, i*d+j), -g(0,j+1)),
925  gmm::mat_col(T_u2_n, i));
926 
927  gmm::add(gmm::scaled(gmm::mat_row(BN1, i), -g(j+1,0)),
928  gmm::mat_col(T_u1_t, i*d+j));
929  if (two_variables)
930  gmm::add(gmm::scaled(gmm::mat_row(BN2, i), -g(j+1,0)),
931  gmm::mat_col(T_u2_t, i*d+j));
932 
933  for (size_type k=0; k < d; ++k) {
934  gmm::add(gmm::scaled(gmm::mat_row(BT1, i*d+k), -g(1+j,1+k)),
935  gmm::mat_col(T_u1_t, i*d+j));
936  if (two_variables)
937  gmm::add(gmm::scaled(gmm::mat_row(BT2, i*d+k), -g(1+j,1+k)),
938  gmm::mat_col(T_u2_t, i*d+j));
939  T_t_t(i*d+j, i*d+k) = g(1+j, 1+k)/r;
940  }
941  T_t_t(i*d+j, i*d+j) -= vt1/(r*alpha[i]);
942  T_t_n(i*d+j, i) = g(1+j,0)/(r*alpha[i]);
943  // T_n_t(i, i*d+j) = g(0,1+j)/(r*alpha[i]);
944  }
945  }
946  gmm::copy(gmm::scaled(BN1, -vt1), T_n_u1);
947  if (two_variables) gmm::copy(gmm::scaled(BN2, -vt1), T_n_u2);
948  gmm::add(gmm::transposed(T_n_u1_transp), T_n_u1);
949  if (two_variables) gmm::add(gmm::transposed(T_n_u2_transp), T_n_u2);
950  gmm::copy(gmm::scaled(BT1, -vt1), T_t_u1);
951  if (two_variables) gmm::copy(gmm::scaled(BT2, -vt1), T_t_u2);
952  break;
953  }
954  }
955 
956  if (version & model::BUILD_RHS) {
957 
958  switch (augmentation_version) {
959  case 1: // unsymmetric Alart-Curnier
960  for (size_type i=0; i < nbc; ++i) {
961  RLN[i] = std::min(scalar_type(0), RLN[i]);
962  if (!contact_only) {
963  scalar_type radius = Tresca_version ? threshold[i]
964  : -friction_coeff[i]*RLN[i];
965  ball_projection
966  (gmm::sub_vector(RLT, gmm::sub_interval(i*d,d)), radius);
967  }
968  }
969  gmm::mult_add(gmm::transposed(BN1), lambda_n, ru1);
970  if (two_variables)
971  gmm::mult_add(gmm::transposed(BN2), lambda_n, ru2);
972  if (!contact_only) {
973  gmm::mult_add(gmm::transposed(BT1), lambda_t, ru1);
974  if (two_variables)
975  gmm::mult_add(gmm::transposed(BT2), lambda_t, ru2);
976  }
977  for (size_type i = 0; i < nbc; ++i) {
978  rlambda_n[i] = (lambda_n[i] - RLN[i]) / (r * alpha[i]);
979  if (!contact_only)
980  for (size_type k = 0; k < d; ++k)
981  rlambda_t[i*d+k]
982  = (lambda_t[i*d+k] - RLT[i*d+k]) / (r * alpha[i]);
983  }
984  break;
985  case 2: // symmetric Alart-Curnier
986  for (size_type i=0; i < nbc; ++i) {
987  RLN[i] = std::min(vt0, RLN[i]);
988  if (!contact_only) {
989  scalar_type radius = Tresca_version ? threshold[i]
990  : -friction_coeff[i]*RLN[i];
991  ball_projection
992  (gmm::sub_vector(RLT, gmm::sub_interval(i*d,d)), radius);
993  }
994  }
995  gmm::mult_add(gmm::transposed(BN1), RLN, ru1);
996  if (two_variables) gmm::mult_add(gmm::transposed(BN2), RLN, ru2);
997  if (!contact_only) {
998  gmm::mult_add(gmm::transposed(BT1), RLT, ru1);
999  if (two_variables) gmm::mult_add(gmm::transposed(BT2), RLT, ru2);
1000  }
1001  for (size_type i = 0; i < nbc; ++i) {
1002  rlambda_n[i] = (lambda_n[i] - RLN[i]) / (r * alpha[i]);
1003  if (!contact_only)
1004  for (size_type k = 0; k < d; ++k)
1005  rlambda_t[i*d+k]
1006  = (lambda_t[i*d+k] - RLT[i*d+k]) / (r * alpha[i]);
1007  }
1008  break;
1009  case 3: // New unsymmetric method
1010  if (!contact_only) gmm::copy(lambda_t, RLT);
1011  for (size_type i=0; i < nbc; ++i) {
1012  RLN[i] = -gmm::neg(lambda_n[i]);
1013  rlambda_n[i] = gmm::pos(lambda_n[i])/r - alpha[i]*gap[i];
1014 
1015  if (!contact_only) {
1016  scalar_type radius = Tresca_version ? threshold[i]
1017  : friction_coeff[i]*gmm::neg(lambda_n[i]);
1018  ball_projection
1019  (gmm::sub_vector(RLT, gmm::sub_interval(i*d,d)), radius);
1020  }
1021  }
1022  gmm::mult(gmm::transposed(BN1), RLN, ru1);
1023  if (two_variables) gmm::mult(gmm::transposed(BN2), RLN, ru2);
1024  gmm::mult_add(BBN1, u1, rlambda_n);
1025  if (two_variables) gmm::mult_add(BBN2, u2, rlambda_n);
1026  if (!contact_only) {
1027  gmm::mult_add(gmm::transposed(BT1), RLT, ru1);
1028  if (two_variables) gmm::mult_add(gmm::transposed(BT2), RLT, ru2);
1029  gmm::add(gmm::scaled(lambda_t, vt1/r), gmm::scaled(RLT,-vt1/r),
1030  rlambda_t);
1031  gmm::mult_add(BBT1, u1, rlambda_t);
1032  if (two_variables) gmm::mult_add(BBT2, u2, rlambda_t);
1033  }
1034  for (size_type i = 0; i < nbc; ++i) {
1035  rlambda_n[i] /= alpha[i];
1036  if (!contact_only)
1037  for (size_type k = 0; k < d; ++k) rlambda_t[i*d+k] /= alpha[i];
1038  }
1039  break;
1040  case 4: // New unsymmetric method with De Saxce projection
1041  base_small_vector x(d+1), n(d+1);
1042  n[0] = vt1;
1043  GMM_ASSERT1(!Tresca_version,
1044  "Augmentation version incompatible with Tresca friction law");
1045  gmm::mult(BBT1, u1, rlambda_t);
1046  if (two_variables)
1047  gmm::mult_add(BBT2, u2, rlambda_t);
1048  for (size_type i=0; i < nbc; ++i) {
1049  x[0] = lambda_n[i];
1050  gmm::copy(gmm::sub_vector(lambda_t, gmm::sub_interval(i*d,d)),
1051  gmm::sub_vector(x, gmm::sub_interval(1, d)));
1052  De_Saxce_projection(x, n, friction_coeff[i]);
1053  RLN[i] = x[0];
1054  gmm::copy(gmm::sub_vector(x, gmm::sub_interval(1, d)),
1055  gmm::sub_vector(RLT, gmm::sub_interval(i*d,d)));
1056  rlambda_n[i] = lambda_n[i]/r - x[0]/r - alpha[i]*gap[i]
1057  - friction_coeff[i] * gmm::vect_norm2(gmm::sub_vector(rlambda_t,
1058  gmm::sub_interval(i*d,d)));
1059  }
1060  gmm::mult_add(gmm::transposed(BT1), RLT, ru1);
1061  if (two_variables) gmm::mult_add(gmm::transposed(BT2), RLT, ru2);
1062  gmm::mult_add(gmm::transposed(BN1), RLN, ru1);
1063  if (two_variables) gmm::mult_add(gmm::transposed(BN2), RLN, ru2);
1064  gmm::add(gmm::scaled(lambda_t, vt1/r), rlambda_t);
1065  gmm::add(gmm::scaled(RLT, -vt1/r), rlambda_t);
1066  gmm::mult_add(BBN1, u1, rlambda_n);
1067  if (two_variables) gmm::mult_add(BBN2, u2, rlambda_n);
1068  for (size_type i = 0; i < nbc; ++i) {
1069  rlambda_n[i] /= alpha[i];
1070  if (!contact_only)
1071  for (size_type k = 0; k < d; ++k) rlambda_t[i*d+k] /= alpha[i];
1072  }
1073  break;
1074  }
1075  }
1076  }
1077 
1078  // specific part for the basic bricks : BN, BT, gap, r, alpha are given.
1079  virtual void asm_real_tangent_terms(const model &md, size_type ib,
1080  const model::varnamelist &vl,
1081  const model::varnamelist &dl,
1082  const model::mimlist &mims,
1083  model::real_matlist &matl,
1084  model::real_veclist &vecl,
1085  model::real_veclist &,
1086  size_type /* region */,
1087  build_version version) const {
1088  if (MPI_IS_MASTER()) {
1089  GMM_ASSERT1(mims.size() == 0, "Contact brick need no mesh_im");
1090  size_type nbvar = 2 + (contact_only ? 0 : 1) + (two_variables ? 1 : 0);
1091  GMM_ASSERT1(vl.size() == nbvar,
1092  "Wrong number of variables for contact brick");
1093  size_type nbdl = 3 + (contact_only ? 0 : 1) + (Tresca_version ? 1 : 0)
1094  + (friction_dynamic_term ? 2 : 0);
1095  GMM_ASSERT1(dl.size() == nbdl, "Wrong number of data for contact brick, "
1096  << dl.size() << " should be " << nbdl);
1097 
1098  size_type nbc = gmm::mat_nrows(BN1);
1099 
1100  // Variables
1101  // Without friction and one displacement : u1, lambda_n
1102  // Without friction and two displacements : u1, u2, lambda_n
1103  // With friction and one displacement : u1, lambda_n, lambda_t
1104  // With friction and two displacements : u1, u2, lambda_n, lambda_t
1105  size_type nv = 0;
1106  const model_real_plain_vector &u1 = md.real_variable(vl[nv++]);
1107  const model_real_plain_vector &u2 = md.real_variable(vl[nv++]);
1108  if (!two_variables) nv--;
1109  const model_real_plain_vector &lambda_n = md.real_variable(vl[nv++]);
1110  if (contact_only) nv--;
1111  const model_real_plain_vector &lambda_t = md.real_variable(vl[nv]);
1112 
1113  // Parameters
1114  // (order : r, gap, alpha, friction_coeff, gamma, wt, threshold)
1115  size_type np = 0, np_wt1 = 0, np_wt2 = 0, np_alpha = 0;
1116  const model_real_plain_vector &vr = md.real_variable(dl[np++]);
1117  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Parameter r should be a scalar");
1118  r = vr[0];
1119  const model_real_plain_vector &vgap = md.real_variable(dl[np++]);
1120  GMM_ASSERT1(gmm::vect_size(vgap) == 1 || gmm::vect_size(vgap) == nbc,
1121  "Parameter gap has a wrong size");
1122  gmm::resize(gap, nbc);
1123  if (gmm::vect_size(vgap) == 1)
1124  gmm::fill(gap, vgap[0]);
1125  else
1126  gmm::copy(vgap, gap);
1127  np_alpha = np++;
1128  const model_real_plain_vector &valpha = md.real_variable(dl[np_alpha]);
1129  GMM_ASSERT1(gmm::vect_size(valpha)== 1 || gmm::vect_size(valpha) == nbc,
1130  "Parameter alpha has a wrong size");
1131  gmm::resize(alpha, nbc);
1132  if (gmm::vect_size(valpha) == 1)
1133  gmm::fill(alpha, valpha[0]);
1134  else
1135  gmm::copy(valpha, alpha);
1136  if (!contact_only) {
1137  const model_real_plain_vector &vfr = md.real_variable(dl[np++]);
1138  GMM_ASSERT1(gmm::vect_size(vfr)==1 || gmm::vect_size(vfr) == nbc,
1139  "Parameter friction_coeff has a wrong size");
1140  gmm::resize(friction_coeff, nbc);
1141  if (gmm::vect_size(vfr) == 1)
1142  gmm::fill(friction_coeff, vfr[0]);
1143  else
1144  gmm::copy(vfr, friction_coeff);
1145  if (friction_dynamic_term) {
1146  const model_real_plain_vector &vg = md.real_variable(dl[np++]);
1147  GMM_ASSERT1(gmm::vect_size(vg) == 1,
1148  "Parameter gamma should be a scalar");
1149  gamma = vg[0];
1150  np_wt1 = np++;
1151  if (two_variables) np_wt2 = np++;
1152  }
1153  if (Tresca_version) {
1154  const model_real_plain_vector &vth = md.real_variable(dl[np++]);
1155  GMM_ASSERT1(gmm::vect_size(vth) == 1 || gmm::vect_size(vth) == nbc,
1156  "Parameter threshold has a wrong size");
1157  gmm::resize(threshold, nbc);
1158  if (gmm::vect_size(vth) == 1)
1159  gmm::fill(threshold, vth[0]);
1160  else
1161  gmm::copy(vth, threshold);
1162  }
1163  }
1164 
1165  if (md.is_var_newer_than_brick(dl[np_alpha], ib)) is_init = false;
1166 
1167  basic_asm_real_tangent_terms
1168  (u1, u2, lambda_n, lambda_t, md.real_variable(dl[np_wt1]),
1169  md.real_variable(dl[np_wt2]), matl, vecl, version);
1170  }
1171 
1172  }
1173 
1174  Coulomb_friction_brick(int aug_version, bool contact_only_,
1175  bool two_variables_=false,
1176  bool Tresca_version_=false,
1177  bool Hughes_stabilized_=false,
1178  bool friction_dynamic_term_=false) {
1179 
1180  if (aug_version == 4 && contact_only_) aug_version = 3;
1181  augmentation_version = aug_version;
1182  GMM_ASSERT1(aug_version >= 1 && aug_version <= 4,
1183  "Wrong augmentation version");
1184  GMM_ASSERT1(!Hughes_stabilized_ || aug_version <= 2,
1185  "The Hughes stabilized version is only for Alart-Curnier "
1186  "version");
1187  contact_only = contact_only_;
1188  is_init = false;
1189  Tresca_version = Tresca_version_;
1190  really_stationary = false; // for future version ...
1191  friction_dynamic_term = friction_dynamic_term_;
1192  two_variables = two_variables_;
1193  Hughes_stabilized = Hughes_stabilized_;
1194  set_flags("Coulomb friction brick", false /* is linear*/,
1195  /* is symmetric */
1196  (augmentation_version == 2) && (contact_only||Tresca_version),
1197  false /* is coercive */, true /* is real */,
1198  false /* is complex */);
1199  }
1200 
1201  void set_BN1(CONTACT_B_MATRIX &BN1_) {
1202  gmm::resize(BN1, gmm::mat_nrows(BN1_), gmm::mat_ncols(BN1_));
1203  gmm::copy(BN1_, BN1);
1204  is_init = false;
1205  }
1206 
1207  void set_BN2(CONTACT_B_MATRIX &BN2_) {
1208  gmm::resize(BN2, gmm::mat_nrows(BN2_), gmm::mat_ncols(BN2_));
1209  gmm::copy(BN2_, BN2);
1210  is_init = false;
1211  }
1212 
1213  void set_DN(CONTACT_B_MATRIX &DN_) {
1214  gmm::resize(DN, gmm::mat_nrows(DN_), gmm::mat_ncols(DN_));
1215  gmm::copy(DN_, DN);
1216  is_init = false;
1217  }
1218 
1219  void set_DT(CONTACT_B_MATRIX &DT_) {
1220  gmm::resize(DT, gmm::mat_nrows(DT_), gmm::mat_ncols(DT_));
1221  gmm::copy(DT_, DT);
1222  is_init = false;
1223  }
1224 
1225  void set_BT1(CONTACT_B_MATRIX &BT1_) {
1226  gmm::resize(BT1, gmm::mat_nrows(BT1_), gmm::mat_ncols(BT1_));
1227  gmm::copy(BT1_, BT1);
1228  is_init = false;
1229  }
1230 
1231  CONTACT_B_MATRIX &get_BN1(void) { return BN1; }
1232  CONTACT_B_MATRIX &get_DN(void) { return DN; }
1233  CONTACT_B_MATRIX &get_DT(void) { return DT; }
1234  CONTACT_B_MATRIX &get_BT1(void) { return BT1; }
1235  const CONTACT_B_MATRIX &get_BN1(void) const { return BN1; }
1236  const CONTACT_B_MATRIX &get_DN(void) const { return DN; }
1237  const CONTACT_B_MATRIX &get_BT1(void) const { return BT1; }
1238 
1239  };
1240 
1241 
1242  void contact_brick_set_stationary(model &md, size_type indbrick) {
1243  pbrick pbr = md.brick_pointer(indbrick);
1244  md.touch_brick(indbrick);
1245  Coulomb_friction_brick *p = dynamic_cast<Coulomb_friction_brick *>
1246  (const_cast<virtual_brick *>(pbr.get()));
1247  GMM_ASSERT1(p, "Wrong type of brick");
1248  p->set_stationary();
1249  }
1250 
1251  CONTACT_B_MATRIX &contact_brick_set_BN
1252  (model &md, size_type indbrick) {
1253  pbrick pbr = md.brick_pointer(indbrick);
1254  md.touch_brick(indbrick);
1255  Coulomb_friction_brick *p = dynamic_cast<Coulomb_friction_brick *>
1256  (const_cast<virtual_brick *>(pbr.get()));
1257  GMM_ASSERT1(p, "Wrong type of brick");
1258  return p->get_BN1();
1259  }
1260 
1261 
1262  CONTACT_B_MATRIX &contact_brick_set_DN
1263  (model &md, size_type indbrick) {
1264  pbrick pbr = md.brick_pointer(indbrick);
1265  md.touch_brick(indbrick);
1266  Coulomb_friction_brick *p = dynamic_cast<Coulomb_friction_brick *>
1267  (const_cast<virtual_brick *>(pbr.get()));
1268  GMM_ASSERT1(p, "Wrong type of brick");
1269  return p->get_DN();
1270  }
1271 
1272  CONTACT_B_MATRIX &contact_brick_set_DT
1273  (model &md, size_type indbrick) {
1274  pbrick pbr = md.brick_pointer(indbrick);
1275  md.touch_brick(indbrick);
1276  Coulomb_friction_brick *p = dynamic_cast<Coulomb_friction_brick *>
1277  (const_cast<virtual_brick *>(pbr.get()));
1278  GMM_ASSERT1(p, "Wrong type of brick");
1279  return p->get_DT();
1280  }
1281 
1282 
1283  CONTACT_B_MATRIX &contact_brick_set_BT
1284  (model &md, size_type indbrick) {
1285  pbrick pbr = md.brick_pointer(indbrick);
1286  md.touch_brick(indbrick);
1287  Coulomb_friction_brick *p = dynamic_cast<Coulomb_friction_brick *>
1288  (const_cast<virtual_brick *>(pbr.get()));
1289  GMM_ASSERT1(p, "Wrong type of brick");
1290  return p->get_BT1();
1291  }
1292 
1293  //=========================================================================
1294  // Add a frictionless contact condition with BN, r, alpha given.
1295  //=========================================================================
1296 
1297  size_type add_basic_contact_brick
1298  (model &md, const std::string &varname_u, const std::string &multname_n,
1299  const std::string &dataname_r, CONTACT_B_MATRIX &BN,
1300  std::string dataname_gap, std::string dataname_alpha,
1301  int aug_version, bool Hughes_stabilized) {
1302 
1303  auto pbr_ = std::make_shared<Coulomb_friction_brick>
1304  (aug_version, true, false, false, Hughes_stabilized);
1305  pbrick pbr(pbr_);
1306  pbr_->set_BN1(BN);
1307 
1308  model::termlist tl;
1309  tl.push_back(model::term_description(varname_u, varname_u, false));
1310  tl.push_back(model::term_description(varname_u, multname_n, false));
1311  tl.push_back(model::term_description(multname_n, varname_u, false));
1312  tl.push_back(model::term_description(multname_n, multname_n, false));
1313  model::varnamelist dl(1, dataname_r);
1314 
1315  if (dataname_gap.size() == 0) {
1316  dataname_gap = md.new_name("contact_gap_on_" + varname_u);
1317  md.add_initialized_fixed_size_data
1318  (dataname_gap, model_real_plain_vector(1, scalar_type(0)));
1319  }
1320  dl.push_back(dataname_gap);
1321 
1322  if (dataname_alpha.size() == 0) {
1323  dataname_alpha = md.new_name("contact_parameter_alpha_on_"+ multname_n);
1324  md.add_initialized_fixed_size_data
1325  (dataname_alpha, model_real_plain_vector(1, scalar_type(1)));
1326  }
1327  dl.push_back(dataname_alpha);
1328 
1329  model::varnamelist vl(1, varname_u);
1330  vl.push_back(multname_n);
1331 
1332  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), size_type(-1));
1333  }
1334 
1335  //=========================================================================
1336  // Add a frictionless contact condition with BN, r, alpha given for
1337  // two deformable bodies
1338  //=========================================================================
1339 
1340  size_type add_basic_contact_brick_two_deformable_bodies
1341  (model &md, const std::string &varname_u1, const std::string &varname_u2,
1342  const std::string &multname_n,
1343  const std::string &dataname_r, CONTACT_B_MATRIX &BN1, CONTACT_B_MATRIX &BN2,
1344  std::string dataname_gap, std::string dataname_alpha,
1345  int aug_version, bool Hughes_stabilized) {
1346 
1347  auto pbr_ = std::make_shared<Coulomb_friction_brick>
1348  (aug_version, true, true, false, Hughes_stabilized);
1349  pbrick pbr(pbr_);
1350  pbr_->set_BN1(BN1);
1351  pbr_->set_BN2(BN2);
1352 
1353  model::termlist tl;
1354  tl.push_back(model::term_description(varname_u1, varname_u1, false));
1355  tl.push_back(model::term_description(varname_u2, varname_u2, false));
1356  tl.push_back(model::term_description(varname_u1, multname_n, false));
1357  tl.push_back(model::term_description(multname_n, varname_u1, false));
1358  tl.push_back(model::term_description(varname_u2, multname_n, false));
1359  tl.push_back(model::term_description(multname_n, varname_u2, false));
1360  tl.push_back(model::term_description(multname_n, multname_n, false));
1361  model::varnamelist dl(1, dataname_r);
1362 
1363  if (dataname_gap.size() == 0) {
1364  dataname_gap = md.new_name("contact_gap_on_" + varname_u1);
1365  md.add_initialized_fixed_size_data
1366  (dataname_gap, model_real_plain_vector(1, scalar_type(0)));
1367  }
1368  dl.push_back(dataname_gap);
1369 
1370  if (dataname_alpha.size() == 0) {
1371  dataname_alpha = md.new_name("contact_parameter_alpha_on_"+ multname_n);
1372  md.add_initialized_fixed_size_data
1373  (dataname_alpha, model_real_plain_vector(1, scalar_type(1)));
1374  }
1375  dl.push_back(dataname_alpha);
1376 
1377  model::varnamelist vl(1, varname_u1);
1378  vl.push_back(varname_u2);
1379  vl.push_back(multname_n);
1380 
1381  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), size_type(-1));
1382  }
1383 
1384 
1385  //=========================================================================
1386  // Add a contact with friction condition with BN, r, alpha given.
1387  //=========================================================================
1388 
1389  size_type add_basic_contact_brick
1390  (model &md, const std::string &varname_u, const std::string &multname_n,
1391  const std::string &multname_t, const std::string &dataname_r,
1392  CONTACT_B_MATRIX &BN, CONTACT_B_MATRIX &BT,
1393  std::string dataname_friction_coeff,
1394  std::string dataname_gap, std::string dataname_alpha,
1395  int aug_version, bool Tresca_version, const std::string dataname_threshold,
1396  std::string dataname_gamma, std::string dataname_wt, bool Hughes_stabilized) {
1397 
1398  bool dynamic_terms = (dataname_gamma.size() > 0);
1399 
1400  auto pbr_ = std::make_shared<Coulomb_friction_brick>
1401  (aug_version,false, false, Tresca_version, Hughes_stabilized,
1402  dynamic_terms);
1403  pbrick pbr(pbr_);
1404  pbr_->set_BN1(BN);
1405  pbr_->set_BT1(BT);
1406 
1407  model::termlist tl;
1408  tl.push_back(model::term_description(varname_u, varname_u, false));
1409  tl.push_back(model::term_description(varname_u, multname_n, false));
1410  tl.push_back(model::term_description(multname_n, varname_u, false));
1411  tl.push_back(model::term_description(multname_n, multname_n, false));
1412  tl.push_back(model::term_description(varname_u, multname_t, false));
1413  tl.push_back(model::term_description(multname_t, varname_u, false));
1414  tl.push_back(model::term_description(multname_t, multname_t, false));
1415  tl.push_back(model::term_description(multname_t, multname_n,
1416  (aug_version == 4)));
1417  model::varnamelist dl(1, dataname_r);
1418  if (dataname_gap.size() == 0) {
1419  dataname_gap = md.new_name("contact_gap_on_" + varname_u);
1420  md.add_initialized_fixed_size_data
1421  (dataname_gap, model_real_plain_vector(1, scalar_type(0)));
1422  }
1423  dl.push_back(dataname_gap);
1424 
1425  if (dataname_alpha.size() == 0) {
1426  dataname_alpha = md.new_name("contact_parameter_alpha_on_"+ multname_n);
1427  md.add_initialized_fixed_size_data
1428  (dataname_alpha, model_real_plain_vector(1, scalar_type(1)));
1429  }
1430  dl.push_back(dataname_alpha);
1431  dl.push_back(dataname_friction_coeff);
1432  if (dataname_gamma.size()) {
1433  dl.push_back(dataname_gamma);
1434  dl.push_back(dataname_wt);
1435  }
1436  if (Tresca_version)
1437  dl.push_back(dataname_threshold);
1438 
1439  model::varnamelist vl(1, varname_u);
1440  vl.push_back(multname_n);
1441  vl.push_back(multname_t);
1442 
1443  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), size_type(-1));
1444  }
1445 
1446 
1447  //=========================================================================
1448  //
1449  // Brick with a given rigid obstacle (one body, build BN, BT, gap, alpha)
1450  //
1451  //=========================================================================
1452  // TODO : add an option for a weak contact condition
1453 
1454  struct Coulomb_friction_brick_rigid_obstacle
1455  : public Coulomb_friction_brick {
1456 
1457  std::string obstacle; // obstacle given with a signed distance expression.
1458 
1459  public :
1460 
1461  virtual void asm_real_tangent_terms(const model &md, size_type ib,
1462  const model::varnamelist &vl,
1463  const model::varnamelist &dl,
1464  const model::mimlist &mims,
1465  model::real_matlist &matl,
1466  model::real_veclist &vecl,
1467  model::real_veclist &,
1468  size_type region,
1469  build_version version) const {
1470  if (MPI_IS_MASTER()) {
1471  GMM_ASSERT1(mims.size() == 1, "This contact brick needs one mesh_im");
1472  size_type nbvar = 2 + (contact_only ? 0 : 1);
1473  GMM_ASSERT1(vl.size() == nbvar,
1474  "Wrong number of variables for contact brick: "
1475  << vl.size() << " should be " << nbvar);
1476  size_type nbdl = 1 + (contact_only ? 0 : 1) + (Tresca_version ? 1 : 0)
1477  + (friction_dynamic_term ? 1 : 0);
1478  GMM_ASSERT1(dl.size() == nbdl,
1479  "Wrong number of data for contact brick: "
1480  << dl.size() << " should be " << nbdl);
1481  GMM_ASSERT1(!two_variables, "internal error");
1482  const mesh_im &mim = *mims[0];
1483 
1484  // Variables
1485  // Without friction and one displacement : u1, lambda_n
1486  // With friction and one displacement : u1, lambda_n, lambda_t
1487  size_type nv = 0;
1488  const model_real_plain_vector &u1 = md.real_variable(vl[nv++]);
1489  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(vl[0]);
1490  const model_real_plain_vector &lambda_n = md.real_variable(vl[nv++]);
1491  if (contact_only) nv--;
1492  const model_real_plain_vector &lambda_t = md.real_variable(vl[nv]);
1493 
1494 
1495  // Parameters (order : r, friction_coeff, gamma, wt, threshold)
1496  size_type np = 0, np_wt1 = 0, nbc;
1497  const model_real_plain_vector &vr = md.real_variable(dl[np++]);
1498  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Parameter r should be a scalar");
1499  r = vr[0];
1500 
1501  // Computation of BN, BT, gap and alpha
1502  if (md.is_var_mf_newer_than_brick(vl[0], ib)) {
1503 
1504  // Verification that mf_u1 is a pure Lagrange fem.
1505  GMM_ASSERT1(!(mf_u1.is_reduced()),
1506  "This contact brick works only for pure Lagrange fems");
1507  dal::bit_vector dofs = mf_u1.basic_dof_on_region(region);
1508  for (dal::bv_visitor id(dofs); !id.finished(); ++id) {
1509  size_type cv = mf_u1.first_convex_of_basic_dof(id);
1510  GMM_ASSERT1(mf_u1.fem_of_element(cv)->is_lagrange(),
1511  "This contact brick works only for pure Lagrange fems");
1512  }
1513  size_type d = mf_u1.get_qdim() - 1, i = 0, j = 0;
1514  nbc = dofs.card() / (d+1);
1515 
1516  // computation of alpha vector.
1517  base_node Pmin, Pmax;
1518  mf_u1.linked_mesh().bounding_box(Pmin, Pmax);
1519  scalar_type l = scalar_type(0);
1520  for (i = 0; i < Pmin.size(); ++i)
1521  l = std::max(l, gmm::abs(Pmax[i] - Pmin[i]));
1522 
1523  CONTACT_B_MATRIX MM(mf_u1.nb_dof(), mf_u1.nb_dof());
1524  asm_mass_matrix(MM, mim, mf_u1, region);
1525  gmm::resize(alpha, nbc);
1526  i = 0; j = 0;
1527  for (dal::bv_visitor id(dofs); !id.finished(); ++id, ++i)
1528  if ((i % (d+1)) == 0) alpha[j++] = MM(id, id) / l;
1529 
1530  getfem::ga_workspace gw;
1531  size_type N = d+1;
1532  getfem::model_real_plain_vector pt(N);
1533  gw.add_fixed_size_constant("X", pt);
1534  if (N >= 1) gw.add_macro("x", "X(1)");
1535  if (N >= 2) gw.add_macro("y", "X(2)");
1536  if (N >= 3) gw.add_macro("z", "X(3)");
1537  if (N >= 4) gw.add_macro("w", "X(4)");
1538 
1539  getfem::ga_function f(gw, obstacle);
1540  f.compile();
1541 
1542  gmm::resize(gap, nbc);
1543  gmm::resize(BN1, nbc, mf_u1.nb_dof());
1544  gmm::clear(BN1);
1545  if (!contact_only) {
1546  gmm::resize(BT1, d*nbc, mf_u1.nb_dof());
1547  gmm::clear(BT1);
1548  }
1549  base_node grad(d+1), ut[3];
1550 
1551  i = 0; j = 0;
1552  for (dal::bv_visitor id(dofs); !id.finished(); ++id, ++i) {
1553  if ((i % (d+1)) == 0) {
1554  gmm::copy(mf_u1.point_of_basic_dof(id), pt);
1555 
1556  // Computation of gap
1557  gap[j] = (f.eval())[0];
1558 
1559  // computation of BN
1560  size_type cv = mf_u1.first_convex_of_basic_dof(id);
1561  scalar_type eps
1562  = mf_u1.linked_mesh().convex_radius_estimate(cv) * 1E-3;
1563  for (size_type k = 0; k <= d; ++k) {
1564  pt[k] += eps;
1565  grad[k] = ((f.eval())[0] - gap[j]) / eps;
1566  pt[k] -= eps;
1567  }
1568  // unit normal vector
1569  base_node un = - grad / gmm::vect_norm2(grad);
1570 
1571  for (size_type k = 0; k <= d; ++k)
1572  BN1(j, id + k) = un[k];
1573 
1574  // computation of BT
1575  if (!contact_only) {
1576 
1577  orthonormal_basis_to_unit_vec(d, un, ut);
1578 
1579  for (size_type k = 0; k <= d; ++k)
1580  for (size_type nn = 0; nn < d; ++nn)
1581  BT1(j*d+nn, id + k) = ut[nn][k];
1582  }
1583  ++j;
1584  }
1585 
1586  }
1587 
1588  GMM_ASSERT1(gmm::vect_size(md.real_variable(vl[1])) == nbc,
1589  "Wrong size of multiplier for the contact condition");
1590 
1591  if (!contact_only)
1592  GMM_ASSERT1(gmm::vect_size(md.real_variable(vl[2])) == nbc*d,
1593  "Wrong size of multiplier for the friction condition");
1594 
1595  is_init = false;
1596  }
1597  else
1598  nbc = gmm::mat_nrows(BN1);
1599 
1600  if (!contact_only) {
1601  const model_real_plain_vector &vfr = md.real_variable(dl[np++]);
1602  GMM_ASSERT1(gmm::vect_size(vfr)==1 || gmm::vect_size(vfr) == nbc,
1603  "Parameter friction_coeff has a wrong size");
1604  gmm::resize(friction_coeff, nbc);
1605  if (gmm::vect_size(vfr) == 1)
1606  gmm::fill(friction_coeff, vfr[0]);
1607  else
1608  gmm::copy(vfr, friction_coeff);
1609  if (friction_dynamic_term) {
1610  const model_real_plain_vector &vg = md.real_variable(dl[np++]);
1611  GMM_ASSERT1(gmm::vect_size(vg) == 1,
1612  "Parameter gamma should be a scalar");
1613  gamma = vg[0];
1614  np_wt1 = np++;
1615  }
1616  if (Tresca_version) {
1617  const model_real_plain_vector &vth = md.real_variable(dl[np++]);
1618  GMM_ASSERT1(gmm::vect_size(vth) == 1 || gmm::vect_size(vth) == nbc,
1619  "Parameter threshold has a wrong size");
1620  gmm::resize(threshold, nbc);
1621  if (gmm::vect_size(vth) == 1)
1622  gmm::fill(threshold, vth[0]);
1623  else
1624  gmm::copy(vth, threshold);
1625  }
1626  }
1627 
1628  basic_asm_real_tangent_terms
1629  (u1, u1, lambda_n, lambda_t, md.real_variable(dl[np_wt1]),
1630  md.real_variable(dl[np_wt1]), matl, vecl, version);
1631  }
1632 
1633  }
1634 
1635  Coulomb_friction_brick_rigid_obstacle
1636  (int aug_version, bool contact_only_, const std::string &obs)
1637  : Coulomb_friction_brick(aug_version, contact_only_), obstacle(obs) {}
1638 
1639  };
1640 
1641 
1642  //=========================================================================
1643  // Add a frictionless contact condition with a rigid obstacle given
1644  // by a signed distance.
1645  //=========================================================================
1646 
1647  size_type add_nodal_contact_with_rigid_obstacle_brick
1648  (model &md, const mesh_im &mim, const std::string &varname_u,
1649  const std::string &multname_n, const std::string &dataname_r,
1650  size_type region, const std::string &obstacle, int aug_version) {
1651  pbrick pbr = std::make_shared<Coulomb_friction_brick_rigid_obstacle>
1652  (aug_version, true, obstacle);
1653 
1654  model::termlist tl;
1655  tl.push_back(model::term_description(varname_u, varname_u, false));
1656  tl.push_back(model::term_description(varname_u, multname_n, false));
1657  tl.push_back(model::term_description(multname_n, varname_u, false));
1658  tl.push_back(model::term_description(multname_n, multname_n, false));
1659  model::varnamelist dl(1, dataname_r);
1660 
1661  model::varnamelist vl(1, varname_u);
1662  vl.push_back(multname_n);
1663 
1664  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1665  }
1666 
1667 
1668  //=========================================================================
1669  // Add a contact with friction condition with a rigid obstacle given
1670  // by a signed distance.
1671  //=========================================================================
1672 
1673  size_type add_nodal_contact_with_rigid_obstacle_brick
1674  (model &md, const mesh_im &mim, const std::string &varname_u,
1675  const std::string &multname_n, const std::string &multname_t,
1676  const std::string &dataname_r, const std::string &dataname_friction_coeff,
1677  size_type region, const std::string &obstacle, int aug_version) {
1678  pbrick pbr = std::make_shared<Coulomb_friction_brick_rigid_obstacle>
1679  (aug_version, false, obstacle);
1680 
1681  model::termlist tl;
1682  tl.push_back(model::term_description(varname_u, varname_u, false));
1683  tl.push_back(model::term_description(varname_u, multname_n, false));
1684  tl.push_back(model::term_description(multname_n, varname_u, false));
1685  tl.push_back(model::term_description(multname_n, multname_n, false));
1686  tl.push_back(model::term_description(varname_u, multname_t, false));
1687  tl.push_back(model::term_description(multname_t, varname_u, false));
1688  tl.push_back(model::term_description(multname_t, multname_t, false));
1689  tl.push_back(model::term_description(multname_t, multname_n,
1690  (aug_version == 4)));
1691  model::varnamelist dl(1, dataname_r);
1692  dl.push_back(dataname_friction_coeff);
1693 
1694  model::varnamelist vl(1, varname_u);
1695  vl.push_back(multname_n);
1696  vl.push_back(multname_t);
1697 
1698  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1699  }
1700 
1701 
1702  //=========================================================================
1703  //
1704  // Brick with elastic bodies (one or two bodies, build BN, BT, Gap, alpha)
1705  //
1706  //=========================================================================
1707  // To be done:
1708  // - Large deformations: what happens when cnpl and nbc change during
1709  // the iterative solution?
1710 
1711  struct Coulomb_friction_brick_nonmatching_meshes
1712  : public Coulomb_friction_brick {
1713 
1714  std::vector<size_type> rg1, rg2; // ids of mesh regions expected to come in
1715  // contact. For one displacement they refer
1716  // both to u1. For two displacements they
1717  // respectively refer to u1, u2.
1718  bool slave1, slave2; // if true, then rg1 or respectively rg2 are treated
1719  // as slave surfaces (the contact multipliers are
1720  // defined on these surfaces)
1721 
1722  virtual void asm_real_tangent_terms(const model &md, size_type ib,
1723  const model::varnamelist &vl,
1724  const model::varnamelist &dl,
1725  const model::mimlist &mims,
1726  model::real_matlist &matl,
1727  model::real_veclist &vecl,
1728  model::real_veclist &,
1729  size_type region,
1730  build_version version) const {
1731  if (MPI_IS_MASTER()) {
1732  GMM_ASSERT1(mims.size() == 2, "This contact brick needs two mesh_im");
1733  const mesh_im &mim1 = *mims[0];
1734  const mesh_im &mim2 = *mims[1];
1735 
1736  // Variables
1737  // Without friction and one displacement : u1, lambda_n
1738  // With friction and one displacement : u1, lambda_n, lambda_t
1739  // Without friction and two displacements : u1, u2, lambda_n
1740  // With friction and two displacements : u1, u2, lambda_n, lambda_t
1741  size_type nv = 0;
1742  std::string varname_u1 = vl[nv];
1743  const model_real_plain_vector &u1 = md.real_variable(varname_u1);
1744  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(varname_u1);
1745  if (two_variables) nv++;
1746  std::string varname_u2 = vl[nv++];
1747  const model_real_plain_vector &u2 = md.real_variable(varname_u2);
1748  const mesh_fem &mf_u2 = md.mesh_fem_of_variable(varname_u2);
1749  const model_real_plain_vector &lambda_n = md.real_variable(vl[nv]);
1750  if (!contact_only) nv++;
1751  const model_real_plain_vector &lambda_t = md.real_variable(vl[nv]);
1752 
1753  size_type nbc = lambda_n.size();
1754 
1755  // Parameters (order: r, friction_coeff)
1756  const model_real_plain_vector &vr = md.real_variable(dl[0]);
1757  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Parameter r should be a scalar");
1758  r = vr[0];
1759  if (!contact_only) {
1760  const model_real_plain_vector &vfr = md.real_variable(dl[1]);
1761  GMM_ASSERT1(gmm::vect_size(vfr)==1 || gmm::vect_size(vfr) == nbc,
1762  "Parameter friction_coeff has a wrong size");
1763  gmm::resize(friction_coeff, nbc);
1764  if (gmm::vect_size(vfr) == 1)
1765  gmm::fill(friction_coeff, vfr[0]);
1766  else
1767  gmm::copy(vfr, friction_coeff);
1768  }
1769 
1770  // Computation of BN, BT, gap and alpha
1771  if ( md.is_var_mf_newer_than_brick(varname_u1, ib)
1772  || md.is_var_mf_newer_than_brick(varname_u2, ib)) {
1773 
1774  for (size_type i = 0; i <= size_type(two_variables ? 1 : 0); i++) {
1775  const mesh_fem &mf_u = i ? mf_u2 : mf_u1;
1776  // Verification that mf_u is a pure Lagrange fem.
1777  GMM_ASSERT1(!(mf_u.is_reduced()),
1778  "This contact brick works only for pure Lagrange fems");
1779  dal::bit_vector dofs = mf_u.basic_dof_on_region(region);
1780  for (dal::bv_visitor id(dofs); !id.finished(); ++id) {
1781  size_type cv = mf_u.first_convex_of_basic_dof(id);
1782  GMM_ASSERT1(mf_u.fem_of_element(cv)->is_lagrange(),
1783  "This contact brick works only for pure Lagrange fems");
1784  }
1785  }
1786 
1787  contact_node_pair_list cnpl;
1788  for (size_type it = 0; it < rg1.size() && it < rg2.size(); ++it)
1789  cnpl.append_min_dist_cn_pairs
1790  (mf_u1, mf_u2, rg1[it], rg2[it], slave1, slave2);
1791 
1792  // Computation of gap, BN and BT
1793  gmm::resize(gap, nbc);
1794  gmm::resize(BN1, nbc, mf_u1.nb_dof());
1795  if (contact_only) {
1796  if (!two_variables) {
1797  compute_contact_matrices(mf_u1, mf_u2, cnpl, gap, &BN1);
1798  } else {
1799  gmm::resize(BN2, nbc, mf_u2.nb_dof());
1800  compute_contact_matrices(mf_u1, mf_u2, cnpl, gap, &BN1, &BN2);
1801  }
1802  } else {
1803  size_type d = mf_u1.get_qdim() - 1;
1804  gmm::resize(BT1, nbc * d, mf_u1.nb_dof());
1805  if (!two_variables) {
1806  compute_contact_matrices(mf_u1, mf_u2, cnpl, gap, &BN1, 0, &BT1);
1807  } else {
1808  // d == mf_u2.get_qdim() - 1;
1809  gmm::resize(BN2, nbc, mf_u2.nb_dof());
1810  gmm::resize(BT2, nbc * d, mf_u2.nb_dof());
1811  compute_contact_matrices(mf_u1, mf_u2, cnpl, gap, &BN1, &BN2, &BT1, &BT2);
1812  }
1813  }
1814 
1815  // computation of alpha vector.
1816  scalar_type l = scalar_type(0);
1817  for (size_type i = 0; i <= size_type(two_variables ? 1 : 0); i++) {
1818  const mesh_fem &mf_u = i ? mf_u2 : mf_u1;
1819  base_node Pmin, Pmax;
1820  mf_u.linked_mesh().bounding_box(Pmin, Pmax);
1821  for (size_type j = 0; j < Pmin.size(); ++j)
1822  l = std::max(l, gmm::abs(Pmax[j] - Pmin[j]));
1823  }
1824  gmm::resize(alpha, nbc);
1825  size_type mult_id = 0;
1826  for (size_type it = 0; it < rg1.size() && it < rg2.size(); ++it) {
1827  for (size_type swap = 0; swap <= 1; ++swap) {
1828  if (swap ? slave2 : slave1) {
1829  size_type rg = swap ? rg2[it] : rg1[it];
1830  const mesh_fem &mf_u = swap ? mf_u2 : mf_u1;
1831  const mesh_im &mim = swap ? mim2 : mim1;
1832  CONTACT_B_MATRIX MM(mf_u.nb_dof(), mf_u.nb_dof());
1833  asm_mass_matrix(MM, mim, mf_u, rg);
1834  size_type qdim = mf_u.get_qdim();
1835  dal::bit_vector rg_dofs = mf_u.basic_dof_on_region(rg);
1836  for (dal::bv_visitor id(rg_dofs); !id.finished(); ++id)
1837  if (id % qdim == 0) alpha[mult_id++] = MM(id, id) / l;
1838  }
1839  }
1840  }
1841  }
1842 
1843  const model_real_plain_vector dummy_wt;
1844  basic_asm_real_tangent_terms
1845  (u1, u2, lambda_n, lambda_t, dummy_wt, dummy_wt, matl, vecl, version);
1846  }
1847  }
1848 
1849  Coulomb_friction_brick_nonmatching_meshes
1850  (int aug_version, bool contact_only_, bool two_variables_,
1851  const std::vector<size_type> &rg1_, const std::vector<size_type> &rg2_,
1852  bool slave1_=true, bool slave2_=false)
1853  : Coulomb_friction_brick(aug_version, contact_only_, two_variables_),
1854  rg1(rg1_), rg2(rg2_), slave1(slave1_), slave2(slave2_) {}
1855 
1856  };
1857 
1858 
1859  //=========================================================================
1860  // Add a frictionless contact condition between two faces of one or two
1861  // elastic bodies.
1862  //=========================================================================
1863 
1864  size_type add_nodal_contact_between_nonmatching_meshes_brick
1865  (model &md, const mesh_im &mim1, const mesh_im &mim2,
1866  const std::string &varname_u1, const std::string &varname_u2,
1867  std::string &multname_n, const std::string &dataname_r,
1868  const std::vector<size_type> &rg1, const std::vector<size_type> &rg2,
1869  bool slave1, bool slave2, int aug_version) {
1870 
1871  bool two_variables = (varname_u1.compare(varname_u2) != 0);
1872 
1873  pbrick pbr = std::make_shared<Coulomb_friction_brick_nonmatching_meshes>
1874  (aug_version, true, two_variables, rg1, rg2, slave1, slave2);
1875 
1876  // Calculate multipliers size
1877  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(varname_u1);
1878  const mesh_fem &mf_u2 = md.mesh_fem_of_variable(varname_u2);
1879  size_type nbc = 0;
1880  for (size_type it = 0; it < rg1.size() && it < rg2.size(); ++it) {
1881  for (size_type swap = 0; swap <= 1; ++swap) {
1882  if (swap ? slave2 : slave1) {
1883  const mesh_fem &mf = swap ? mf_u2 : mf_u1;
1884  size_type rg = swap ? rg2[it] : rg1[it];
1885  dal::bit_vector rg_dofs = mf.basic_dof_on_region(rg);
1886  nbc += rg_dofs.card() / mf.get_qdim();
1887  }
1888  }
1889  }
1890 
1891  if (multname_n.size() == 0)
1892  multname_n = md.new_name("contact_multiplier");
1893  else
1894  GMM_ASSERT1(multname_n.compare(md.new_name(multname_n)) == 0,
1895  "The given name for the multiplier is already reserved in the model");
1896  md.add_fixed_size_variable(multname_n, nbc);
1897 
1898  model::termlist tl;
1899  tl.push_back(model::term_description(varname_u1, varname_u1, false));
1900  if (two_variables) {
1901  tl.push_back(model::term_description(varname_u2, varname_u2, false));
1902  }
1903  tl.push_back(model::term_description(varname_u1, multname_n, false));
1904  tl.push_back(model::term_description(multname_n, varname_u1, false));
1905  if (two_variables) {
1906  tl.push_back(model::term_description(varname_u2, multname_n, false));
1907  tl.push_back(model::term_description(multname_n, varname_u2, false));
1908  }
1909  tl.push_back(model::term_description(multname_n, multname_n, false));
1910 
1911  // Variables (order: varname_u, multname_n)
1912  model::varnamelist vl;
1913  vl.push_back(varname_u1);
1914  if (two_variables) vl.push_back(varname_u2);
1915  vl.push_back(multname_n);
1916 
1917  // Parameters (order: r, ...)
1918  model::varnamelist dl;
1919  dl.push_back(dataname_r);
1920 
1921  model::mimlist ml;
1922  ml.push_back(&mim1);
1923  ml.push_back(&mim2);
1924 
1925  return md.add_brick(pbr, vl, dl, tl, ml, size_type(-1));
1926  }
1927 
1928 
1929  //=========================================================================
1930  // Add a contact with friction condition between two faces of one or two
1931  // elastic bodies.
1932  //=========================================================================
1933 
1934  size_type add_nodal_contact_between_nonmatching_meshes_brick
1935  (model &md, const mesh_im &mim1, const mesh_im &mim2,
1936  const std::string &varname_u1, const std::string &varname_u2,
1937  std::string &multname_n, std::string &multname_t,
1938  const std::string &dataname_r, const std::string &dataname_friction_coeff,
1939  const std::vector<size_type> &rg1, const std::vector<size_type> &rg2,
1940  bool slave1, bool slave2, int aug_version) {
1941 
1942  bool two_variables = (varname_u1.compare(varname_u2) != 0);
1943 
1944  pbrick pbr = std::make_shared<Coulomb_friction_brick_nonmatching_meshes>
1945  (aug_version, false, two_variables, rg1, rg2, slave1, slave2);
1946 
1947  // Calculate multipliers size
1948  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(varname_u1);
1949  const mesh_fem &mf_u2 = md.mesh_fem_of_variable(varname_u2);
1950  size_type nbc = 0;
1951  for (size_type it = 0; it < rg1.size() && it < rg2.size(); ++it) {
1952  for (size_type swap = 0; swap <= 1; ++swap) {
1953  if (swap ? slave2 : slave1) {
1954  const mesh_fem &mf = swap ? mf_u2 : mf_u1;
1955  size_type rg = swap ? rg2[it] : rg1[it];
1956  dal::bit_vector rg_dofs = mf.basic_dof_on_region(rg);
1957  nbc += rg_dofs.card() / mf.get_qdim();
1958  }
1959  }
1960  }
1961 
1962  if (multname_n.size() == 0)
1963  multname_n = md.new_name("contact_normal_multiplier");
1964  else
1965  GMM_ASSERT1(multname_n.compare(md.new_name(multname_n)) == 0,
1966  "The given name for the multiplier is already reserved in the model");
1967  md.add_fixed_size_variable(multname_n, nbc);
1968  if (multname_t.size() == 0)
1969  multname_t = md.new_name("contact_tangent_multiplier");
1970  else
1971  GMM_ASSERT1(multname_t.compare(md.new_name(multname_t)) == 0,
1972  "The given name for the multiplier is already reserved in the model");
1973  md.add_fixed_size_variable(multname_t, nbc * (mf_u1.get_qdim() - 1) ); // ??
1974 
1975  model::termlist tl;
1976  tl.push_back(model::term_description(varname_u1, varname_u1, false));
1977  if (two_variables) {
1978  tl.push_back(model::term_description(varname_u2, varname_u2, false));
1979  }
1980 
1981  tl.push_back(model::term_description(varname_u1, multname_n, false));
1982  tl.push_back(model::term_description(multname_n, varname_u1, false));
1983  if (two_variables) {
1984  tl.push_back(model::term_description(varname_u2, multname_n, false));
1985  tl.push_back(model::term_description(multname_n, varname_u2, false));
1986  }
1987  tl.push_back(model::term_description(multname_n, multname_n, false));
1988 
1989  tl.push_back(model::term_description(varname_u1, multname_t, false));
1990  tl.push_back(model::term_description(multname_t, varname_u1, false));
1991  if (two_variables) {
1992  tl.push_back(model::term_description(varname_u2, multname_t, false));
1993  tl.push_back(model::term_description(multname_t, varname_u2, false));
1994  }
1995  tl.push_back(model::term_description(multname_t, multname_t, false));
1996  tl.push_back(model::term_description(multname_t, multname_n,
1997  (aug_version == 4)));
1998 
1999  // Variables (order: varname_u, multname_n, multname_t)
2000  model::varnamelist vl;
2001  vl.push_back(varname_u1);
2002  if (two_variables) vl.push_back(varname_u2);
2003  vl.push_back(multname_n);
2004  vl.push_back(multname_t);
2005 
2006  // Parameters (order: r, friction_coeff)
2007  model::varnamelist dl;
2008  dl.push_back(dataname_r);
2009  dl.push_back(dataname_friction_coeff);
2010 
2011  model::mimlist ml;
2012  ml.push_back(&mim1);
2013  ml.push_back(&mim2);
2014 
2015  return md.add_brick(pbr, vl, dl, tl, ml, size_type(-1));
2016  }
2017 
2018 } /* end of namespace getfem. */
bgeot_kdtree.h
Simple implementation of a KD-tree.
bgeot::convex::dir_points_of_face
dref_convex_pt_ct dir_points_of_face(short_type i) const
Direct points for a given face.
Definition: bgeot_convex.h:84
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:63
bgeot::geotrans_inv_convex::invert
bool invert(const base_node &n, base_node &n_ref, scalar_type IN_EPS=1e-12, bool project_into_element=false)
given the node on the real element, returns the node on the reference element (even if it is outside ...
Definition: bgeot_geotrans_inv.cc:60
bgeot::kdtree
Balanced tree over a set of points.
Definition: bgeot_kdtree.h:101
bgeot::kdtree::add_point_with_id
void add_point_with_id(const base_node &n, size_type i)
insert a new point, with an associated number.
Definition: bgeot_kdtree.h:119
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:151
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:315
getfem::mesh_fem::basic_dof_on_region
virtual dal::bit_vector basic_dof_on_region(const mesh_region &b) const
Get a list of basic dof lying on a given mesh_region.
Definition: getfem_mesh_fem.cc:76
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:46
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
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::mesh_fem_of_variable
const mesh_fem & mesh_fem_of_variable(const std::string &name) const
Gives the access to the mesh_fem of a variable if any.
Definition: getfem_models.cc:2888
getfem::model::add_initialized_fixed_size_data
void add_initialized_fixed_size_data(const std::string &name, const VECT &v)
Add a fixed size data (assumed to be a vector) to the model and initialized with v.
Definition: getfem_models.h:740
getfem::model::add_fixed_size_variable
void add_fixed_size_variable(const std::string &name, size_type size, size_type niter=1)
Add a fixed size variable to the model assumed to be a vector.
Definition: getfem_models.cc:725
getfem::model::new_name
std::string new_name(const std::string &name)
Gives a non already existing variable name begining by name.
Definition: getfem_models.cc:174
getfem::model::touch_brick
void touch_brick(size_type ib)
Force the re-computation of a brick for the next assembly.
Definition: getfem_models.h:1028
getfem::virtual_brick
The virtual brick has to be derived to describe real model bricks.
Definition: getfem_models.h:1474
gmm::tab_ref_index_ref
indexed array reference (given a container X, and a set of indexes I, this class provides a pseudo-co...
Definition: gmm_ref.h:303
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_nodal.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::fill
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*‍/
Definition: gmm_blas.h:103
gmm::vect_norm2_sqr
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2_sqr(const V &v)
squared Euclidean norm of a vector.
Definition: gmm_blas.h:543
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
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::lu_det
linalg_traits< DenseMatrixLU >::value_type lu_det(const DenseMatrixLU &LU, const Pvector &pvector)
Compute the matrix determinant (via a LU factorization)
Definition: gmm_dense_lu.h:240
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::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:243
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
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_basic_contact_brick
size_type add_basic_contact_brick(model &md, const std::string &varname_u, const std::string &multname_n, const std::string &dataname_r, CONTACT_B_MATRIX &BN, std::string dataname_gap="", std::string dataname_alpha="", int aug_version=1, bool Hughes_stabilized=false)
Add a frictionless contact condition to the model.
Definition: getfem_contact_and_friction_nodal.cc:1298
getfem::add_basic_contact_brick_two_deformable_bodies
size_type add_basic_contact_brick_two_deformable_bodies(model &md, const std::string &varname_u1, const std::string &varname_u2, const std::string &multname_n, const std::string &dataname_r, CONTACT_B_MATRIX &BN1, CONTACT_B_MATRIX &BN2, std::string dataname_gap="", std::string dataname_alpha="", int aug_version=1, bool Hughes_stabilized=false)
Add a frictionless contact condition to the model between two deformable bodies.
Definition: getfem_contact_and_friction_nodal.cc:1341
getfem::contact_brick_set_BN
CONTACT_B_MATRIX & contact_brick_set_BN(model &md, size_type indbrick)
Can be used to change the matrix BN of a basic contact/friction brick.
Definition: getfem_contact_and_friction_nodal.cc:1252
getfem::contact_brick_set_stationary
void contact_brick_set_stationary(model &md, size_type indbrick)
Can be used to set the stationary option.
Definition: getfem_contact_and_friction_nodal.cc:1242
getfem::contact_brick_set_DT
CONTACT_B_MATRIX & contact_brick_set_DT(model &md, size_type indbrick)
Can be used to change the matrix DT of a basic contact/friction brick.
Definition: getfem_contact_and_friction_nodal.cc:1273
getfem::contact_brick_set_DN
CONTACT_B_MATRIX & contact_brick_set_DN(model &md, size_type indbrick)
Can be used to change the matrix DN of a basic contact/friction brick.
Definition: getfem_contact_and_friction_nodal.cc:1263
getfem::pbrick
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
Definition: getfem_models.h:48
getfem::add_nodal_contact_between_nonmatching_meshes_brick
size_type add_nodal_contact_between_nonmatching_meshes_brick(model &md, const mesh_im &mim1, const mesh_im &mim2, const std::string &varname_u1, const std::string &varname_u2, std::string &multname_n, const std::string &dataname_r, const std::vector< size_type > &rg1, const std::vector< size_type > &rg2, bool slave1=true, bool slave2=false, int aug_version=1)
Add a frictionless contact condition between two faces of one or two elastic bodies.
Definition: getfem_contact_and_friction_nodal.cc:1865
getfem::add_nodal_contact_with_rigid_obstacle_brick
size_type add_nodal_contact_with_rigid_obstacle_brick(model &md, const mesh_im &mim, const std::string &varname_u, const std::string &multname_n, const std::string &dataname_r, size_type region, const std::string &obstacle, int aug_version=1)
Add a frictionless contact condition with a rigid obstacle to the model.
Definition: getfem_contact_and_friction_nodal.cc:1648
getfem::contact_brick_set_BT
CONTACT_B_MATRIX & contact_brick_set_BT(model &md, size_type indbrick)
Can be used to change the matrix BT of a basic contact/friction brick.
Definition: getfem_contact_and_friction_nodal.cc:1284
bgeot::index_node_pair
store a point and the associated index for the kdtree.
Definition: bgeot_kdtree.h:58