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