GetFEM  5.4.2
getfem_contact_and_friction_nodal.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-2020 Yves Renard, Konstantinos Poulios.
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 
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_neighbors(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_neighbors[vv1].size(); ++j)
201  if (pt1_neighbors[vv1][j] == vv2) {
202  already_in = true;
203  break;
204  }
205  if (!already_in) pt1_neighbors[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_neighbors[i1].size(); ++j) {
213  size_type i2 = pt1_neighbors[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  getfem::ga_workspace gw;
1532  size_type N = d+1;
1533  getfem::model_real_plain_vector pt(N);
1534  gw.add_fixed_size_constant("X", pt);
1535  if (N >= 1) gw.add_macro("x", "X(1)");
1536  if (N >= 2) gw.add_macro("y", "X(2)");
1537  if (N >= 3) gw.add_macro("z", "X(3)");
1538  if (N >= 4) gw.add_macro("w", "X(4)");
1539 
1540  getfem::ga_function f(gw, obstacle);
1541  f.compile();
1542 
1543  gmm::resize(gap, nbc);
1544  gmm::resize(BN1, nbc, mf_u1.nb_dof());
1545  gmm::clear(BN1);
1546  if (!contact_only) {
1547  gmm::resize(BT1, d*nbc, mf_u1.nb_dof());
1548  gmm::clear(BT1);
1549  }
1550  base_node grad(d+1), ut[3];
1551 
1552  i = 0; j = 0;
1553  for (dal::bv_visitor id(dofs); !id.finished(); ++id, ++i) {
1554  if ((i % (d+1)) == 0) {
1555  gmm::copy(mf_u1.point_of_basic_dof(id), pt);
1556 
1557  // Computation of gap
1558  gap[j] = (f.eval())[0];
1559 
1560  // computation of BN
1561  size_type cv = mf_u1.first_convex_of_basic_dof(id);
1562  scalar_type eps
1563  = mf_u1.linked_mesh().convex_radius_estimate(cv) * 1E-3;
1564  for (size_type k = 0; k <= d; ++k) {
1565  pt[k] += eps;
1566  grad[k] = ((f.eval())[0] - gap[j]) / eps;
1567  pt[k] -= eps;
1568  }
1569  // unit normal vector
1570  base_node un = - grad / gmm::vect_norm2(grad);
1571 
1572  for (size_type k = 0; k <= d; ++k)
1573  BN1(j, id + k) = un[k];
1574 
1575  // computation of BT
1576  if (!contact_only) {
1577 
1578  orthonormal_basis_to_unit_vec(d, un, ut);
1579 
1580  for (size_type k = 0; k <= d; ++k)
1581  for (size_type nn = 0; nn < d; ++nn)
1582  BT1(j*d+nn, id + k) = ut[nn][k];
1583  }
1584  ++j;
1585  }
1586 
1587  }
1588 
1589  GMM_ASSERT1(gmm::vect_size(md.real_variable(vl[1])) == nbc,
1590  "Wrong size of multiplier for the contact condition");
1591 
1592  if (!contact_only)
1593  GMM_ASSERT1(gmm::vect_size(md.real_variable(vl[2])) == nbc*d,
1594  "Wrong size of multiplier for the friction condition");
1595 
1596  is_init = false;
1597  }
1598  else
1599  nbc = gmm::mat_nrows(BN1);
1600 
1601  if (!contact_only) {
1602  const model_real_plain_vector &vfr = md.real_variable(dl[np++]);
1603  GMM_ASSERT1(gmm::vect_size(vfr)==1 || gmm::vect_size(vfr) == nbc,
1604  "Parameter friction_coeff has a wrong size");
1605  gmm::resize(friction_coeff, nbc);
1606  if (gmm::vect_size(vfr) == 1)
1607  gmm::fill(friction_coeff, vfr[0]);
1608  else
1609  gmm::copy(vfr, friction_coeff);
1610  if (friction_dynamic_term) {
1611  const model_real_plain_vector &vg = md.real_variable(dl[np++]);
1612  GMM_ASSERT1(gmm::vect_size(vg) == 1,
1613  "Parameter gamma should be a scalar");
1614  gamma = vg[0];
1615  np_wt1 = np++;
1616  }
1617  if (Tresca_version) {
1618  const model_real_plain_vector &vth = md.real_variable(dl[np++]);
1619  GMM_ASSERT1(gmm::vect_size(vth) == 1 || gmm::vect_size(vth) == nbc,
1620  "Parameter threshold has a wrong size");
1621  gmm::resize(threshold, nbc);
1622  if (gmm::vect_size(vth) == 1)
1623  gmm::fill(threshold, vth[0]);
1624  else
1625  gmm::copy(vth, threshold);
1626  }
1627  }
1628 
1629  basic_asm_real_tangent_terms
1630  (u1, u1, lambda_n, lambda_t, md.real_variable(dl[np_wt1]),
1631  md.real_variable(dl[np_wt1]), matl, vecl, version);
1632  }
1633 
1634  }
1635 
1636  Coulomb_friction_brick_rigid_obstacle
1637  (int aug_version, bool contact_only_, const std::string &obs)
1638  : Coulomb_friction_brick(aug_version, contact_only_), obstacle(obs) {}
1639 
1640  };
1641 
1642 
1643  //=========================================================================
1644  // Add a frictionless contact condition with a rigid obstacle given
1645  // by a signed distance.
1646  //=========================================================================
1647 
1648  size_type add_nodal_contact_with_rigid_obstacle_brick
1649  (model &md, const mesh_im &mim, const std::string &varname_u,
1650  const std::string &multname_n, const std::string &dataname_r,
1651  size_type region, const std::string &obstacle, int aug_version) {
1652  pbrick pbr = std::make_shared<Coulomb_friction_brick_rigid_obstacle>
1653  (aug_version, true, obstacle);
1654 
1655  model::termlist tl;
1656  tl.push_back(model::term_description(varname_u, varname_u, false));
1657  tl.push_back(model::term_description(varname_u, multname_n, false));
1658  tl.push_back(model::term_description(multname_n, varname_u, false));
1659  tl.push_back(model::term_description(multname_n, multname_n, false));
1660  model::varnamelist dl(1, dataname_r);
1661 
1662  model::varnamelist vl(1, varname_u);
1663  vl.push_back(multname_n);
1664 
1665  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1666  }
1667 
1668 
1669  //=========================================================================
1670  // Add a contact with friction condition with a rigid obstacle given
1671  // by a signed distance.
1672  //=========================================================================
1673 
1674  size_type add_nodal_contact_with_rigid_obstacle_brick
1675  (model &md, const mesh_im &mim, const std::string &varname_u,
1676  const std::string &multname_n, const std::string &multname_t,
1677  const std::string &dataname_r, const std::string &dataname_friction_coeff,
1678  size_type region, const std::string &obstacle, int aug_version) {
1679  pbrick pbr = std::make_shared<Coulomb_friction_brick_rigid_obstacle>
1680  (aug_version, false, obstacle);
1681 
1682  model::termlist tl;
1683  tl.push_back(model::term_description(varname_u, varname_u, false));
1684  tl.push_back(model::term_description(varname_u, multname_n, false));
1685  tl.push_back(model::term_description(multname_n, varname_u, false));
1686  tl.push_back(model::term_description(multname_n, multname_n, false));
1687  tl.push_back(model::term_description(varname_u, multname_t, false));
1688  tl.push_back(model::term_description(multname_t, varname_u, false));
1689  tl.push_back(model::term_description(multname_t, multname_t, false));
1690  tl.push_back(model::term_description(multname_t, multname_n,
1691  (aug_version == 4)));
1692  model::varnamelist dl(1, dataname_r);
1693  dl.push_back(dataname_friction_coeff);
1694 
1695  model::varnamelist vl(1, varname_u);
1696  vl.push_back(multname_n);
1697  vl.push_back(multname_t);
1698 
1699  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1700  }
1701 
1702 
1703  //=========================================================================
1704  //
1705  // Brick with elastic bodies (one or two bodies, build BN, BT, Gap, alpha)
1706  //
1707  //=========================================================================
1708  // To be done:
1709  // - Large deformations: what happens when cnpl and nbc change during
1710  // the iterative solution?
1711 
1712  struct Coulomb_friction_brick_nonmatching_meshes
1713  : public Coulomb_friction_brick {
1714 
1715  std::vector<size_type> rg1, rg2; // ids of mesh regions expected to come in
1716  // contact. For one displacement they refer
1717  // both to u1. For two displacements they
1718  // respectively refer to u1, u2.
1719  bool slave1, slave2; // if true, then rg1 or respectively rg2 are treated
1720  // as slave surfaces (the contact multipliers are
1721  // defined on these surfaces)
1722 
1723  virtual void asm_real_tangent_terms(const model &md, size_type ib,
1724  const model::varnamelist &vl,
1725  const model::varnamelist &dl,
1726  const model::mimlist &mims,
1727  model::real_matlist &matl,
1728  model::real_veclist &vecl,
1729  model::real_veclist &,
1730  size_type region,
1731  build_version version) const {
1732  if (MPI_IS_MASTER()) {
1733  GMM_ASSERT1(mims.size() == 2, "This contact brick needs two mesh_im");
1734  const mesh_im &mim1 = *mims[0];
1735  const mesh_im &mim2 = *mims[1];
1736 
1737  // Variables
1738  // Without friction and one displacement : u1, lambda_n
1739  // With friction and one displacement : u1, lambda_n, lambda_t
1740  // Without friction and two displacements : u1, u2, lambda_n
1741  // With friction and two displacements : u1, u2, lambda_n, lambda_t
1742  size_type nv = 0;
1743  std::string varname_u1 = vl[nv];
1744  const model_real_plain_vector &u1 = md.real_variable(varname_u1);
1745  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(varname_u1);
1746  if (two_variables) nv++;
1747  std::string varname_u2 = vl[nv++];
1748  const model_real_plain_vector &u2 = md.real_variable(varname_u2);
1749  const mesh_fem &mf_u2 = md.mesh_fem_of_variable(varname_u2);
1750  const model_real_plain_vector &lambda_n = md.real_variable(vl[nv]);
1751  if (!contact_only) nv++;
1752  const model_real_plain_vector &lambda_t = md.real_variable(vl[nv]);
1753 
1754  size_type nbc = lambda_n.size();
1755 
1756  // Parameters (order: r, friction_coeff)
1757  const model_real_plain_vector &vr = md.real_variable(dl[0]);
1758  GMM_ASSERT1(gmm::vect_size(vr) == 1, "Parameter r should be a scalar");
1759  r = vr[0];
1760  if (!contact_only) {
1761  const model_real_plain_vector &vfr = md.real_variable(dl[1]);
1762  GMM_ASSERT1(gmm::vect_size(vfr)==1 || gmm::vect_size(vfr) == nbc,
1763  "Parameter friction_coeff has a wrong size");
1764  gmm::resize(friction_coeff, nbc);
1765  if (gmm::vect_size(vfr) == 1)
1766  gmm::fill(friction_coeff, vfr[0]);
1767  else
1768  gmm::copy(vfr, friction_coeff);
1769  }
1770 
1771  // Computation of BN, BT, gap and alpha
1772  if ( md.is_var_mf_newer_than_brick(varname_u1, ib)
1773  || md.is_var_mf_newer_than_brick(varname_u2, ib)) {
1774 
1775  for (size_type i = 0; i <= size_type(two_variables ? 1 : 0); i++) {
1776  const mesh_fem &mf_u = i ? mf_u2 : mf_u1;
1777  // Verification that mf_u is a pure Lagrange fem.
1778  GMM_ASSERT1(!(mf_u.is_reduced()),
1779  "This contact brick works only for pure Lagrange fems");
1780  dal::bit_vector dofs = mf_u.basic_dof_on_region(region);
1781  for (dal::bv_visitor id(dofs); !id.finished(); ++id) {
1782  size_type cv = mf_u.first_convex_of_basic_dof(id);
1783  GMM_ASSERT1(mf_u.fem_of_element(cv)->is_lagrange(),
1784  "This contact brick works only for pure Lagrange fems");
1785  }
1786  }
1787 
1788  contact_node_pair_list cnpl;
1789  for (size_type it = 0; it < rg1.size() && it < rg2.size(); ++it)
1790  cnpl.append_min_dist_cn_pairs
1791  (mf_u1, mf_u2, rg1[it], rg2[it], slave1, slave2);
1792 
1793  // Computation of gap, BN and BT
1794  gmm::resize(gap, nbc);
1795  gmm::resize(BN1, nbc, mf_u1.nb_dof());
1796  if (contact_only) {
1797  if (!two_variables) {
1798  compute_contact_matrices(mf_u1, mf_u2, cnpl, gap, &BN1);
1799  } else {
1800  gmm::resize(BN2, nbc, mf_u2.nb_dof());
1801  compute_contact_matrices(mf_u1, mf_u2, cnpl, gap, &BN1, &BN2);
1802  }
1803  } else {
1804  size_type d = mf_u1.get_qdim() - 1;
1805  gmm::resize(BT1, nbc * d, mf_u1.nb_dof());
1806  if (!two_variables) {
1807  compute_contact_matrices(mf_u1, mf_u2, cnpl, gap, &BN1, 0, &BT1);
1808  } else {
1809  // d == mf_u2.get_qdim() - 1;
1810  gmm::resize(BN2, nbc, mf_u2.nb_dof());
1811  gmm::resize(BT2, nbc * d, mf_u2.nb_dof());
1812  compute_contact_matrices(mf_u1, mf_u2, cnpl, gap, &BN1, &BN2, &BT1, &BT2);
1813  }
1814  }
1815 
1816  // computation of alpha vector.
1817  scalar_type l = scalar_type(0);
1818  for (size_type i = 0; i <= size_type(two_variables ? 1 : 0); i++) {
1819  const mesh_fem &mf_u = i ? mf_u2 : mf_u1;
1820  base_node Pmin, Pmax;
1821  mf_u.linked_mesh().bounding_box(Pmin, Pmax);
1822  for (size_type j = 0; j < Pmin.size(); ++j)
1823  l = std::max(l, gmm::abs(Pmax[j] - Pmin[j]));
1824  }
1825  gmm::resize(alpha, nbc);
1826  size_type mult_id = 0;
1827  for (size_type it = 0; it < rg1.size() && it < rg2.size(); ++it) {
1828  for (size_type swap = 0; swap <= 1; ++swap) {
1829  if (swap ? slave2 : slave1) {
1830  size_type rg = swap ? rg2[it] : rg1[it];
1831  const mesh_fem &mf_u = swap ? mf_u2 : mf_u1;
1832  const mesh_im &mim = swap ? mim2 : mim1;
1833  CONTACT_B_MATRIX MM(mf_u.nb_dof(), mf_u.nb_dof());
1834  asm_mass_matrix(MM, mim, mf_u, rg);
1835  size_type qdim = mf_u.get_qdim();
1836  dal::bit_vector rg_dofs = mf_u.basic_dof_on_region(rg);
1837  for (dal::bv_visitor id(rg_dofs); !id.finished(); ++id)
1838  if (id % qdim == 0) alpha[mult_id++] = MM(id, id) / l;
1839  }
1840  }
1841  }
1842  }
1843 
1844  const model_real_plain_vector dummy_wt;
1845  basic_asm_real_tangent_terms
1846  (u1, u2, lambda_n, lambda_t, dummy_wt, dummy_wt, matl, vecl, version);
1847  }
1848  }
1849 
1850  Coulomb_friction_brick_nonmatching_meshes
1851  (int aug_version, bool contact_only_, bool two_variables_,
1852  const std::vector<size_type> &rg1_, const std::vector<size_type> &rg2_,
1853  bool slave1_=true, bool slave2_=false)
1854  : Coulomb_friction_brick(aug_version, contact_only_, two_variables_),
1855  rg1(rg1_), rg2(rg2_), slave1(slave1_), slave2(slave2_) {}
1856 
1857  };
1858 
1859 
1860  //=========================================================================
1861  // Add a frictionless contact condition between two faces of one or two
1862  // elastic bodies.
1863  //=========================================================================
1864 
1865  size_type add_nodal_contact_between_nonmatching_meshes_brick
1866  (model &md, const mesh_im &mim1, const mesh_im &mim2,
1867  const std::string &varname_u1, const std::string &varname_u2,
1868  std::string &multname_n, const std::string &dataname_r,
1869  const std::vector<size_type> &rg1, const std::vector<size_type> &rg2,
1870  bool slave1, bool slave2, int aug_version) {
1871 
1872  bool two_variables = (varname_u1.compare(varname_u2) != 0);
1873 
1874  pbrick pbr = std::make_shared<Coulomb_friction_brick_nonmatching_meshes>
1875  (aug_version, true, two_variables, rg1, rg2, slave1, slave2);
1876 
1877  // Calculate multipliers size
1878  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(varname_u1);
1879  const mesh_fem &mf_u2 = md.mesh_fem_of_variable(varname_u2);
1880  size_type nbc = 0;
1881  for (size_type it = 0; it < rg1.size() && it < rg2.size(); ++it) {
1882  for (size_type swap = 0; swap <= 1; ++swap) {
1883  if (swap ? slave2 : slave1) {
1884  const mesh_fem &mf = swap ? mf_u2 : mf_u1;
1885  size_type rg = swap ? rg2[it] : rg1[it];
1886  dal::bit_vector rg_dofs = mf.basic_dof_on_region(rg);
1887  nbc += rg_dofs.card() / mf.get_qdim();
1888  }
1889  }
1890  }
1891 
1892  if (multname_n.size() == 0)
1893  multname_n = md.new_name("contact_multiplier");
1894  else
1895  GMM_ASSERT1(multname_n.compare(md.new_name(multname_n)) == 0,
1896  "The given name for the multiplier is already reserved in the model");
1897  md.add_fixed_size_variable(multname_n, nbc);
1898 
1899  model::termlist tl;
1900  tl.push_back(model::term_description(varname_u1, varname_u1, false));
1901  if (two_variables) {
1902  tl.push_back(model::term_description(varname_u2, varname_u2, false));
1903  }
1904  tl.push_back(model::term_description(varname_u1, multname_n, false));
1905  tl.push_back(model::term_description(multname_n, varname_u1, false));
1906  if (two_variables) {
1907  tl.push_back(model::term_description(varname_u2, multname_n, false));
1908  tl.push_back(model::term_description(multname_n, varname_u2, false));
1909  }
1910  tl.push_back(model::term_description(multname_n, multname_n, false));
1911 
1912  // Variables (order: varname_u, multname_n)
1913  model::varnamelist vl;
1914  vl.push_back(varname_u1);
1915  if (two_variables) vl.push_back(varname_u2);
1916  vl.push_back(multname_n);
1917 
1918  // Parameters (order: r, ...)
1919  model::varnamelist dl;
1920  dl.push_back(dataname_r);
1921 
1922  model::mimlist ml;
1923  ml.push_back(&mim1);
1924  ml.push_back(&mim2);
1925 
1926  return md.add_brick(pbr, vl, dl, tl, ml, size_type(-1));
1927  }
1928 
1929 
1930  //=========================================================================
1931  // Add a contact with friction condition between two faces of one or two
1932  // elastic bodies.
1933  //=========================================================================
1934 
1935  size_type add_nodal_contact_between_nonmatching_meshes_brick
1936  (model &md, const mesh_im &mim1, const mesh_im &mim2,
1937  const std::string &varname_u1, const std::string &varname_u2,
1938  std::string &multname_n, std::string &multname_t,
1939  const std::string &dataname_r, const std::string &dataname_friction_coeff,
1940  const std::vector<size_type> &rg1, const std::vector<size_type> &rg2,
1941  bool slave1, bool slave2, int aug_version) {
1942 
1943  bool two_variables = (varname_u1.compare(varname_u2) != 0);
1944 
1945  pbrick pbr = std::make_shared<Coulomb_friction_brick_nonmatching_meshes>
1946  (aug_version, false, two_variables, rg1, rg2, slave1, slave2);
1947 
1948  // Calculate multipliers size
1949  const mesh_fem &mf_u1 = md.mesh_fem_of_variable(varname_u1);
1950  const mesh_fem &mf_u2 = md.mesh_fem_of_variable(varname_u2);
1951  size_type nbc = 0;
1952  for (size_type it = 0; it < rg1.size() && it < rg2.size(); ++it) {
1953  for (size_type swap = 0; swap <= 1; ++swap) {
1954  if (swap ? slave2 : slave1) {
1955  const mesh_fem &mf = swap ? mf_u2 : mf_u1;
1956  size_type rg = swap ? rg2[it] : rg1[it];
1957  dal::bit_vector rg_dofs = mf.basic_dof_on_region(rg);
1958  nbc += rg_dofs.card() / mf.get_qdim();
1959  }
1960  }
1961  }
1962 
1963  if (multname_n.size() == 0)
1964  multname_n = md.new_name("contact_normal_multiplier");
1965  else
1966  GMM_ASSERT1(multname_n.compare(md.new_name(multname_n)) == 0,
1967  "The given name for the multiplier is already reserved in the model");
1968  md.add_fixed_size_variable(multname_n, nbc);
1969  if (multname_t.size() == 0)
1970  multname_t = md.new_name("contact_tangent_multiplier");
1971  else
1972  GMM_ASSERT1(multname_t.compare(md.new_name(multname_t)) == 0,
1973  "The given name for the multiplier is already reserved in the model");
1974  md.add_fixed_size_variable(multname_t, nbc * (mf_u1.get_qdim() - 1) ); // ??
1975 
1976  model::termlist tl;
1977  tl.push_back(model::term_description(varname_u1, varname_u1, false));
1978  if (two_variables) {
1979  tl.push_back(model::term_description(varname_u2, varname_u2, false));
1980  }
1981 
1982  tl.push_back(model::term_description(varname_u1, multname_n, false));
1983  tl.push_back(model::term_description(multname_n, varname_u1, false));
1984  if (two_variables) {
1985  tl.push_back(model::term_description(varname_u2, multname_n, false));
1986  tl.push_back(model::term_description(multname_n, varname_u2, false));
1987  }
1988  tl.push_back(model::term_description(multname_n, multname_n, false));
1989 
1990  tl.push_back(model::term_description(varname_u1, multname_t, false));
1991  tl.push_back(model::term_description(multname_t, varname_u1, false));
1992  if (two_variables) {
1993  tl.push_back(model::term_description(varname_u2, multname_t, false));
1994  tl.push_back(model::term_description(multname_t, varname_u2, false));
1995  }
1996  tl.push_back(model::term_description(multname_t, multname_t, false));
1997  tl.push_back(model::term_description(multname_t, multname_n,
1998  (aug_version == 4)));
1999 
2000  // Variables (order: varname_u, multname_n, multname_t)
2001  model::varnamelist vl;
2002  vl.push_back(varname_u1);
2003  if (two_variables) vl.push_back(varname_u2);
2004  vl.push_back(multname_n);
2005  vl.push_back(multname_t);
2006 
2007  // Parameters (order: r, friction_coeff)
2008  model::varnamelist dl;
2009  dl.push_back(dataname_r);
2010  dl.push_back(dataname_friction_coeff);
2011 
2012  model::mimlist ml;
2013  ml.push_back(&mim1);
2014  ml.push_back(&mim2);
2015 
2016  return md.add_brick(pbr, vl, dl, tl, ml, size_type(-1));
2017  }
2018 
2019 } /* end of namespace getfem. */
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:312
bgeot_kdtree.h
Simple implementation of a KD-tree.
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:64
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
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::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::add_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:1866
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:175
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
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::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::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
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:726
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot::kdtree
Balanced tree over a set of points.
Definition: bgeot_kdtree.h:102
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
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::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:743
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:209
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:1031
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::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:244
getfem_contact_and_friction_nodal.h
Unilateral contact and Coulomb friction condition brick.
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:2891
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:1649
getfem::model::add_brick
size_type add_brick(pbrick pbr, const varnamelist &varnames, const varnamelist &datanames, const termlist &terms, const mimlist &mims, size_type region)
Add a brick to the model.
Definition: getfem_models.cc:1033
getfem::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::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:85
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
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
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::virtual_brick
The virtual brick has to be derived to describe real model bricks.
Definition: getfem_models.h:1440
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:61
bgeot::index_node_pair
store a point and the associated index for the kdtree.
Definition: bgeot_kdtree.h:59
getfem_contact_and_friction_common.h
Comomon tools for unilateral contact and Coulomb friction bricks.
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
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
getfem::pbrick
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
Definition: getfem_models.h:49
gmm::mult_add
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1781
getfem_assembling.h
Miscelleanous assembly routines for common terms. Use the low-level generic assembly....
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:120