GetFEM  5.4.2
getfem_level_set_contact.cc
1 /*===========================================================================
2 
3  Copyright (C) 2012-2020 Andriy Andreykiv
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 
25 #include <algorithm>
29 #include <math.h>
30 
31 level_set_contact::contact_body::contact_body(model& _md, std::string _var_name):
32  var_name(_var_name),
33  is_deformed(false),
34  own_mesh(const_cast<mesh&>(_md.mesh_fem_of_variable(_var_name).
35  linked_mesh())),
36  own_mesh_fem(_md.mesh_fem_of_variable(_var_name)),
37  md(_md)
38 {}
39 
40 
41 
43  getfem::model& _md, const std::string& _var_name,
44  getfem::mesh_im* _pmim_ls) : contact_body(_md,_var_name),
45  ls_name("ls_on_"+_var_name),
46  ls_mesh_fem(_md.mesh_fem_of_variable(_var_name)),
47  pmim(_pmim_ls)
48 
49 {
50  ls_mesh_fem.set_qdim(size_type(1));
51  model_real_plain_vector LS(ls_mesh_fem.nb_dof());
52  boundary_level_set_field(get_mesh(),ls_mesh_fem,
53  *pmim,LS);
54  md.add_initialized_fem_data(ls_name,ls_mesh_fem,LS);
55 }
56 
57 
59  getfem::model& _md, std::string _var_name, std::string _ls_name):
60 contact_body(_md,_var_name),
61  ls_name(_ls_name),
62  pmim(0)
63 { }
64 
66 {
67  for(size_type i=0;i<ls_values().size();i++)
68  ls_values()[i]+=off;
69 }
70 
72  model& _md,
73  const std::string& _var_name,
74  size_type _mult_order,
75  size_type _mult_mim_order) :
76 
77  contact_body(_md,_var_name),
78  mult_mim_order(_mult_mim_order),
79  mult_int_method(""),
80  mult_mf_order(_mult_order),
81  integration(PER_ELEMENT),
82  regularized_tollerance(0),
83  small_weight_multiplier(0),
84  max_contact_angle(45)
85 {
86  //store existing elements in VOLUME_ELEMENTS region
87  //before boundary elements are created
88  VOLUME_ELEMENTS = getfem::mesh_region::free_region_id(get_mesh());
89  get_mesh().region(VOLUME_ELEMENTS).add(get_mesh().convex_index());
90 
91  //create boundary elements (current mesh_fem should be automatically
92  // extended with these elements)
93  BOUNDARY_ELEMENTS = getfem::mesh_region::free_region_id(get_mesh());
94  get_mesh().region(BOUNDARY_ELEMENTS).clear();
95 
96  masters.push_back(this);
97 }
98 
100  model& _md,
101  const std::string& _var_name,
102  size_type _mult_order,
103  const std::string& _mult_mim_method,
104  contact_integration _integration,
105  scalar_type _regularized_tollerance,
106  scalar_type _small_weight_multiplier,
107  scalar_type _max_contact_angle):
108 
109  contact_body(_md,_var_name),
110  mult_mim_order(size_type(-1)),
111  mult_int_method(_mult_mim_method),
112  mult_mf_order(_mult_order),
113  integration(_integration),
114  regularized_tollerance(_regularized_tollerance),
115  small_weight_multiplier(_small_weight_multiplier),
116  max_contact_angle(_max_contact_angle)
117 {
118  //store existing elements in VOLUME_ELEMENTS region
119  //before boundary elements are created
120  VOLUME_ELEMENTS = getfem::mesh_region::free_region_id(get_mesh());
121  get_mesh().region(VOLUME_ELEMENTS).add(get_mesh().convex_index());
122 
123  //create boundary elements (current mesh_fem should be automatically
124  // extended with these elements)
125  BOUNDARY_ELEMENTS = getfem::mesh_region::free_region_id(get_mesh());
126  get_mesh().region(BOUNDARY_ELEMENTS).clear();
127  masters.push_back(this);
128 }
129 
130 
133  const std::string& slave_var_name) const
134 {
135  auto it = contact_table.find(slave_var_name);
136  if (it!=contact_table.end()) return *(it->second);
137  GMM_ASSERT1(false,"did not find info on slave contact body, \
138  defined on variable "+slave_var_name);
139 }
140 
143  const std::string& slave_var_name)
144 {
145  auto it = contact_table.find(slave_var_name);
146  if (it!=contact_table.end()) return *(it->second);
147  GMM_ASSERT1(false,"did not find info on slave contact body, \
148  defined on variable "+slave_var_name);
149 }
150 
151 
152 level_set_contact::face_type level_set_contact::master_contact_body::
154 {
155  auto it = border_faces.find(i);
156  if (it!=border_faces.end()) return it->second;
157  GMM_ASSERT1(false,"did not find a face, corresponding to element "<<i);
158 }
159 
160 
162  add_slave(slave_contact_body& scb, size_type assumed_contact_region)
163 {
164  //check input
165  GMM_ASSERT1(&md==&scb.get_model(),
166  "Model objects of master and slave are not the same");
167  if (assumed_contact_region!=size_type(-1))
168  GMM_ASSERT1(get_mesh().region(assumed_contact_region).is_boundary(),
169  "Assumed_contact_region must be on the boundary");
170 
171  //add surface elements where contact will be computed
172  size_type assumed_contact_elems = getfem::mesh_region::free_region_id(get_mesh());
173  getfem::mesh_region& contact_elems = get_mesh().region(assumed_contact_elems);
174  getfem::mesh_region& boundary_elems = get_mesh().region(BOUNDARY_ELEMENTS);
175  std::shared_ptr<getfem::mr_visitor> i;
176  getfem::mesh_region outer_faces;
177  outer_faces.clear();
178  getfem::outer_faces_of_mesh(get_mesh(), outer_faces);
179 
180  if (assumed_contact_region==size_type(-1)){ //all faces will be searched for contact
181  i = std::make_shared<getfem::mr_visitor>(outer_faces);
182  }
183  else // only specified faces will be searched
184  {
185  getfem::mesh_region& assumed_region = get_mesh().
186  region(assumed_contact_region);
187  i = std::make_shared<getfem::mr_visitor>(assumed_region);
188  }
189 
190  for (; !i->finished(); ++(*i)){
191  getfem::size_type new_elem =
192  get_mesh().add_convex(
194  get_mesh().trans_of_convex(i->cv())),
195  get_mesh().ind_points_of_face_of_convex(
196  i->cv(), i->f()).begin());
197 
198 
199  border_faces[new_elem] = face_type(*i);
200  contact_elems.add(new_elem);
201  boundary_elems.add(new_elem);
202  }
203 
204  GMM_ASSERT1(get_mesh().region(BOUNDARY_ELEMENTS).index().card()!=0,
205  "No boundary elements added !!!");
206 
207  GMM_ASSERT1(get_mesh().region(assumed_contact_elems).index().card()!=0,
208  "No contact elements added !!!");
209 
210  //creating Lagrange multiplier
211  std::string mult_name = md.new_name("mult_on_"+get_var_name()+
212  "_and_"+scb.get_var_name());
213  const mesh_fem &mf_mult =
214  getfem::classical_mesh_fem(get_mesh(),
215  bgeot::dim_type(mult_mf_order), bgeot::dim_type(1));
216  md.add_multiplier(mult_name,mf_mult,get_var_name());
217 
218  //adding variable to store level set, projected from the slave
219  const mesh_fem& mf_ls =
220  getfem::classical_mesh_fem(get_mesh(),bgeot::dim_type(mult_mf_order+1));
221  plain_vector LS(mf_ls.nb_dof());
222  md.add_initialized_fem_data("ls_on"+get_var_name()+
223  "_from_"+scb.get_var_name(),mf_ls,LS);
224 
225  //register contact pair
226  contact_table[scb.get_var_name()] =
227  std::make_shared<contact_pair_info>(*this,scb,mult_name,assumed_contact_elems);
228 
229 }
230 
231 std::vector<level_set_contact::master_contact_body*>
232  level_set_contact::master_contact_body::masters;
233 
234 
236 {
237  if (masters.size()==0) GMM_WARNING3("Running contact detection, while no \
238  contact bodies are registered");
239  bool contact_surfaces_changed = false;
240  for(size_type i=0;i<masters.size();i++)
241  if (masters[i]->master_contact_changed())
242  contact_surfaces_changed=true;
243 
244  return contact_surfaces_changed;
245 }
246 
248 {
249  if (masters.size()==0) GMM_WARNING3("Clearing contact lists, while no \
250  contact bodies are registered");
251  for(size_type i=0;i<masters.size();i++)
252  masters[i]->clear_contact_history();
253 }
254 
256 {
257  std::map<std::string, std::shared_ptr<contact_pair_info> >::
258  iterator it = contact_table.begin();
259  for(;it!=contact_table.end();it++)
260  it->second->clear_contact_history();
261 }
262 
264 {
265  bool contact_surfaces_changed = false;
266  std::map<std::string, std::shared_ptr<contact_pair_info> >::
267  iterator it = contact_table.begin();
268  for(;it!=contact_table.end();it++)
269  if (it->second->contact_changed())
270  contact_surfaces_changed=true;
271 
272  return contact_surfaces_changed;
273 }
274 
275 std::shared_ptr<getfem::mesh_im> level_set_contact::master_contact_body::
277 {
278 
279  std::shared_ptr<getfem::mesh_im> pmim_contact;
280 
281  pmim_contact = std::make_shared<mesh_im>(get_mesh());
282  if (mult_mim_order!=size_type(-1)){
283  pmim_contact->set_integration_method
284  (get_mesh().region(region_id).index(),
285  bgeot::dim_type(contact_mim_order()));
286  }
287  else
288  {
289  pmim_contact->set_integration_method
290  (get_mesh().region(region_id).index(), contact_int_method());
291  }
292 
293  return pmim_contact;
294 }
295 
296 
297 level_set_contact::contact_pair_update::contact_pair_update(
298  master_contact_body& _mcb,
299  slave_contact_body& _scb,
300  update_depth ud):
301 mcb(_mcb), scb(_scb)
302 
303 {
304 
305  GMM_ASSERT1(!mcb.is_mesh_deformed(),"Trying to deform \
306  already deformed Master Contact Body");
307  GMM_ASSERT1(!scb.is_mesh_deformed(),"Trying to deform \
308  already deformed Slave Contact Body");
309 
310  const model_real_plain_vector&
311  Umaster=mcb.get_model().real_variable(mcb.get_var_name());
312  // size_type dof_check = Umaster.size();
313  // size_type node_check = mcb.get_mesh().nb_points();
314  def_master = std::make_shared<getfem::temporary_mesh_deformator<>>
315  (mcb.get_mesh_fem(),Umaster);
316  mcb.is_deformed=true;
317  if (&mcb.get_mesh()!=&scb.get_mesh()){
318  // not deforming the slave if the master and the slave are the same
319  const model_real_plain_vector&
320  Uslave=scb.get_model().real_variable(scb.get_var_name());
321  def_slave = std::make_shared<getfem::temporary_mesh_deformator<>>
322  (scb.get_mesh_fem(),Uslave);
323  scb.is_deformed=true;
324  }
325  if (ud == FULL_UPDATE) mcb.update_for_slave(scb.get_var_name());
326 }
327 
328 level_set_contact::contact_pair_update::~contact_pair_update(){
329  mcb.is_deformed=false;
330  scb.is_deformed=false;
331 }
332 
333 
334 
335 level_set_contact::contact_pair_info::contact_pair_info(
336  master_contact_body& underformed_mcb,
337  slave_contact_body& underformed_scb,
338  const std::string& _mult_name,
339  size_type _GIVEN_CONTACT_REGION) :
340 
341  master_cb(underformed_mcb),
342  slave_cb(underformed_scb),
343  mult_name(_mult_name),
344  GIVEN_CONTACT_REGION(_GIVEN_CONTACT_REGION),
345 
346  ACTIVE_CONTACT_REGION(getfem::mesh_region::free_region_id(master_cb.get_mesh())),
347  members_are_computed(false),
348  init_cont_detect_done(false)
349 
350 {
351  //input check (if mult_name is incorrect, exception will be generated)
352  // const mesh_fem& mf_mult=
353  // master_cb.get_model().mesh_fem_of_variable(mult_name);
354  GMM_ASSERT1(master_cb.get_mesh().
355  region(GIVEN_CONTACT_REGION).index().card()!=0,
356  "provided contact region for contact_pair_info class is empty!!!");
357 }
358 
360 {
361  old_contact_elm_list.clear();
362  pre_old_ct_list.clear();
363 }
364 
366 {
367  //deform master and slave meshes
369  temp_mesh_deformation(master_cb,slave_cb,DEFORM_MESHES_ONLY);
370 
371  // create mf on the boundary of the master (copy from the master)
372  mesh_fem mf_scalar(master_cb.get_mesh());
373  for(size_type i=0;i<mf_scalar.linked_mesh().nb_convex();i++)
374  mf_scalar.set_finite_element(i,master_cb.get_mesh_fem().fem_of_element(i));
375 
376  mf_scalar.set_qdim(1);
377  getfem::partial_mesh_fem mf_boundary(mf_scalar);
378  mf_boundary.adapt(mf_scalar.dof_on_region(GIVEN_CONTACT_REGION));
379 
380  // interpolate level set from the slave to the master
381  model_real_plain_vector LS_on_contour(mf_boundary.nb_dof());
382  getfem::interpolation(slave_cb.get_ls_mesh_fem(), mf_boundary,
383  slave_cb.ls_values(), LS_on_contour);
384  model_real_plain_vector LS(mf_scalar.nb_dof());
385  mf_boundary.extend_vector(LS_on_contour,LS);
386  gmm::copy(LS,master_cb.get_model().set_real_variable
387  ("ls_on"+master_cb.get_var_name()+"_from_"+slave_cb.get_var_name()));
388 
389  // interpolate the gradient of the level set onto the master surfaces
390  // (this is to obtain the normal direction of the level set)
391  mesh_fem mf_gradient_ls(slave_cb.get_mesh());
392  mf_gradient_ls.set_classical_discontinuous_finite_element(bgeot::dim_type(master_cb.mult_mf_order));
393  mesh_fem mf_gradient_ls_vect(mf_gradient_ls);
394  mf_gradient_ls_vect.set_qdim(slave_cb.get_mesh().dim());
395  plain_vector GradLS(mf_gradient_ls.nb_dof()*slave_cb.get_mesh().dim());
396  getfem::compute_gradient(slave_cb.get_ls_mesh_fem(), mf_gradient_ls, slave_cb.ls_values(), GradLS);
397  getfem::partial_mesh_fem mf_boundary_vect(master_cb.get_mesh_fem());
398  mf_boundary_vect.adapt(master_cb.get_mesh_fem().dof_on_region(GIVEN_CONTACT_REGION));
399  plain_vector GradLS_boundary(mf_boundary.nb_dof()*slave_cb.get_mesh().dim());
400  getfem::interpolation(mf_gradient_ls_vect,mf_boundary_vect,GradLS,GradLS_boundary);
401 
402  size_type dim = slave_cb.get_mesh().dim();
403  const scalar_type TINY = 1e-15;
404  for(size_type i=0;i<mf_boundary.nb_dof();i++){ //normalizing the projected ls field
405  bgeot::base_node ls_node(dim);
406  for(size_type j=0;j<dim;j++) ls_node[j]=GradLS_boundary[dim*i+j];
407  ls_node/= (gmm::vect_norm2(ls_node)+TINY);
408  for(size_type j=0;j<dim;j++) GradLS_boundary[dim*i+j]=ls_node[j];
409  }
410  plain_vector normLS_master(master_cb.get_mesh_fem().nb_dof());
411  mf_boundary_vect.extend_vector(GradLS_boundary,normLS_master);
412 
413 
414  // extend Lagrange Multiplier onto the whole boundary of the master
415  const mesh_fem& mf_mult =
416  master_cb.get_model().mesh_fem_of_variable(mult_name);
417  const model_real_plain_vector& lambda =
418  master_cb.get_model().real_variable(mult_name);
419  model_real_plain_vector lambda_full(mf_mult.nb_basic_dof());
420  if (lambda.size()>0) mf_mult.extend_vector(lambda,lambda_full);
421  // update contact region
422  dal::bit_vector cc = master_cb.get_mesh().
423  region(GIVEN_CONTACT_REGION).index();
424  master_cb.get_mesh().region(ACTIVE_CONTACT_REGION).clear();
425  master_cb.get_mesh().region(ACTIVE_CONTACT_REGION).add(cc);
427  for (o << cc; o != bgeot::size_type(-1); o << cc) {
428  getfem::mesh_fem::ind_dof_ct dof_ls = mf_scalar.ind_basic_dof_of_element(o);
429  getfem::mesh_fem::ind_dof_ct dof_lm = mf_mult.ind_basic_dof_of_element(o);
430 
431  //measure the angle between ls countour and the master face
432  face_type face = master_cb.ext_face_of_elem(o);
433  bgeot::base_node unit_face_normal =
434  master_cb.get_mesh().normal_of_face_of_convex(face.cv,face.f);
435  unit_face_normal/=gmm::vect_norm2(unit_face_normal);
436  scalar_type cosine_alpha = 0;
437  for (size_type j = 0; j < dof_ls.size(); j++){
438  bgeot::base_node ls_grad_node(dim);
439  for(size_type k=0; k<dim; k++)
440  ls_grad_node[k]=normLS_master[dim*dof_ls[j]+k];
441  cosine_alpha += gmm::vect_sp(ls_grad_node,unit_face_normal);
442  }
443  cosine_alpha /= scalar_type(dof_ls.size());
444  scalar_type alpha = acos(cosine_alpha)*360/(2*M_PI); // now this is average angle
445  // between master surface and ls zero contour
446 
447  scalar_type LS_extreeme = LS[dof_ls[0]];
448  if (master_cb.integration==master_contact_body::PER_ELEMENT)
449  for (size_type j = 0; j < dof_ls.size(); j++)
450  LS_extreeme=std::min(LS[dof_ls[j]],LS_extreeme);
451  else
452  for (size_type j = 0; j < dof_ls.size(); j++)
453  LS_extreeme=std::max(LS[dof_ls[j]],LS_extreeme);
454 
455  scalar_type LM_sum = 0;
456  for (size_type j = 0; j < dof_lm.size(); j++)
457  LM_sum+=lambda_full[dof_lm[j]];
458 
459  const scalar_type TINY_2 = 1e-9;
460 
461  if (LS_extreeme+LM_sum < TINY_2 || alpha > master_cb.max_contact_angle)
462  master_cb.get_mesh().region(ACTIVE_CONTACT_REGION).sup(o);
463  }
464 
465  // check whether contact areas have changed
466  bool contact_surface_changed;
467  const dal::bit_vector& current_contact_elm_list =
468  master_cb.get_mesh().region(ACTIVE_CONTACT_REGION).index();
469  GMM_TRACE2("Current contact elements: "<< current_contact_elm_list);
470  GMM_TRACE2("Old contact elements: "<< old_contact_elm_list);
471  GMM_TRACE2("Pre-old contact elements: "<< pre_old_ct_list);
472 
473  if (current_contact_elm_list == old_contact_elm_list &&
474  current_contact_elm_list.card() == old_contact_elm_list.card()) {
475  contact_surface_changed = false;
476  GMM_TRACE2(" the contact area has not changed");
477  } else {
478  if (current_contact_elm_list == pre_old_ct_list &&
479  current_contact_elm_list.card() == pre_old_ct_list.card()) {
480  contact_surface_changed = false;
481  GMM_TRACE2(" the contact area has changed, but cycling, \
482  so exiting active set search");
483  } else {
484  contact_surface_changed = true;
485  GMM_TRACE2(" the contact area has changed");
486  pre_old_ct_list = old_contact_elm_list;
487  old_contact_elm_list = current_contact_elm_list;
488  }
489  }
490 
491  init_cont_detect_done = true;
492  force_update();
493 
494 
495  //building integration method
496  pmim_contact = master_cb.build_mesh_im_on_boundary(ACTIVE_CONTACT_REGION);
497  n_integrated_elems = pmim_contact->convex_index().card();
498  GMM_ASSERT1(n_integrated_elems==current_contact_elm_list.card(),
499  "Failure in integration method: The number of integrated elements "
500  "does not correspond to the number of contact elements");
501 
502  return contact_surface_changed;
503 
504 }
505 
506 
507 
509 {
510  GMM_ASSERT1(master_cb.is_mesh_deformed(),"Master mesh is not deformed, \
511  cannot calucalte contact info");
512 
513  GMM_ASSERT1(slave_cb.is_mesh_deformed(),"Slave mesh is not deformed, \
514  cannot calucalte contact info");
515 
516  GMM_ASSERT1(master_cb.get_mesh().region(ACTIVE_CONTACT_REGION).index().card()>0,
517  "Internal error: Contact area is empty");
518 
519  //pinterpolated_fem for level set
520  pinterpolated_fem = std::make_shared<mesh_fem>(master_cb.get_mesh());
521 
522  pinterpolated_fem_U = std::shared_ptr<mesh_fem>();
523 
524  if (ifem_srf.get()) getfem::del_interpolated_fem(ifem_srf);
525 
526  ifem_srf = getfem::new_interpolated_fem(slave_cb.get_ls_mesh_fem(),
527  *pmim_contact, 0, dal::bit_vector(), false);
528  pinterpolated_fem->set_finite_element(
529  master_cb.get_mesh().region(ACTIVE_CONTACT_REGION).index(),ifem_srf);
530  pinterpolated_fem->set_qdim(1);
531 
532 
533  //pinterpolated_fem_U
534  pinterpolated_fem_U = std::make_shared<mesh_fem>(master_cb.get_mesh());
535  pinterpolated_fem_U->set_finite_element(master_cb.get_mesh().
536  region(ACTIVE_CONTACT_REGION).index(),ifem_srf);
537  pinterpolated_fem_U->set_qdim(master_cb.get_mesh().dim());
538 
539  //slave_ls_dofs
540  std::vector<size_type> index(pinterpolated_fem->nb_dof());
541  dal::bit_vector cc =
542  master_cb.get_mesh().region(ACTIVE_CONTACT_REGION).index();
543  for (dal::bv_visitor icv(cc); !icv.finished(); ++icv){
544  for (size_type j = 0; j < pinterpolated_fem->nb_basic_dof_of_element(icv);
545  ++j)
546  {index[pinterpolated_fem->ind_basic_dof_of_element(icv)[j]]
547  = ifem_srf->index_of_global_dof(icv, j);}
548  }
549 
550  slave_ls_dofs = std::make_shared<gmm::unsorted_sub_index>(index);
551 
552  //slave_U_dofs
553  std::vector<size_type> indexU(pinterpolated_fem_U->nb_dof());
554  size_type dim = pinterpolated_fem_U->get_qdim();
555  for(size_type d=0;d<dim;d++)
556  for(size_type i=0;i<pinterpolated_fem->nb_dof();i++)
557  indexU[dim*i+d] = dim*index[i]+d;
558  slave_U_dofs = std::make_shared<gmm::unsorted_sub_index>(indexU);
559 
560  members_are_computed=true;
561 }
562 
563 
564 
566  model& md,
567  master_contact_body& mcb,
568  slave_contact_body& scb,
569  size_type rg) {
570  //level set contact class
571  getfem::pbrick pbr = std::make_shared<level_set_contact_brick>(md,mcb,scb,rg);
572 
573  //term description
574  const std::string& name_Um = mcb.get_var_name();
575  const std::string& name_Us = scb.get_var_name();
576  const std::string& name_LM = mcb.get_pair_info(name_Us).get_mult_name();
577  model::termlist terms;
578  terms.push_back(model::term_description(name_Um,name_Um,false));
579  terms.push_back(model::term_description(name_Us,name_Us,false));
580  terms.push_back(model::term_description(name_LM,name_LM,false));
581  terms.push_back(model::term_description(name_Um,name_Us,true));
582  terms.push_back(model::term_description(name_Um,name_LM,true));
583  terms.push_back(model::term_description(name_Us,name_LM,true));
584 
585  //variables
586  model::varnamelist variables;
587  variables.push_back(name_Um);
588  variables.push_back(name_Us);
589  variables.push_back(name_LM);
590 
591  //empty data and integration method lists
592  //(we don't have any properties or initial data,
593  //while integration methods are created per iteration)
594  model::varnamelist datalist;
595  model::mimlist mimlist;
596 
597  //register the brick with the model and return its number
598  return md.add_brick(pbr,variables,datalist,terms,mimlist,rg);
599 }
600 
601 
602 level_set_contact::level_set_contact_brick::
603  level_set_contact_brick(
604  model& _md,
605  master_contact_body& _mcb,
606  slave_contact_body& _scb,
607  size_type rg) :
608 md(_md),mcb(_mcb),scb(_scb), given_contact_id(rg)
609 {
610  GMM_ASSERT1(&md == &mcb.get_model(),
611  "Master body is defined on a different model then the input");
612 
613  //register master/slave pair
614  mcb.add_slave(scb,given_contact_id);
615 
616  //Reduce computation to own MPI region
617  contact_region_id = mcb.get_pair_info(scb.get_var_name()).contact_region();
618  getfem::mesh_region& contact_region_ = mcb.get_mesh().region(contact_region_id);
619  mcb.get_mesh().intersect_with_mpi_region(contact_region_);
620  contact_region_id = contact_region_.id(); //probably not needed, but still
621 
622 
623  set_flags("Level set contact brick",
624  false /* is linear*/,
625  true /* is symmetric */,
626  false /* is coercive */,
627  true /* is real */,
628  false /* is complex */);
629 
630 }
631 
633  const model &mdd, size_type /* ib */,
634  const model::varnamelist &vl,
635  const model::varnamelist &/* dl */,
636  const model::mimlist &/* mims */,
637  model::real_matlist &matl,
638  model::real_veclist &vecl,
639  model::real_veclist &,
640  size_type region,
641  build_version version) const {
642  //input check
643  GMM_ASSERT1(vl.size() == 3,
644  "Level set contact brick needs three variables");
645  GMM_ASSERT1(matl.size() == 6,
646  "Level set contact brick needs six matrices");
647  GMM_ASSERT1(vecl.size() == 6,
648  "Level set contact brick assembles size RHSs");
649  GMM_ASSERT1(region==given_contact_id,
650  "Assumed contact region has changed!!! \
651  This implementation does not handle this \
652  for efficiency reasons!!");
653 
654  if (version & model::BUILD_MATRIX )
655  for(size_type i=0;i<matl.size();i++) gmm::clear(matl[i]);
656  if (version & model::BUILD_RHS )
657  for(size_type i=0;i<vecl.size();i++) gmm::clear(vecl[i]);
658 
659  const getfem::mesh_region& active_contact_region =
660  mcb.get_mesh().region(contact_region_id);
661  if (active_contact_region.index().card()==0) return; //no contact -> no contact assembly
662 
663  //deform the meshes, update contact info
664  contact_pair_update cp_update(mcb,scb,FULL_UPDATE);
665 
666  //extract DOF vectors
667  const plain_vector &LM = mdd.real_variable(vl[2]);
668 
669  //Assemble Tangent Matrix
670  if (version & model::BUILD_MATRIX ) {
671  GMM_TRACE2("Level set contact brick stiffness matrix assembly on "
672  << mcb.get_pair_info(scb.get_var_name()).num_of_integr_elems()
673  << " elements");
674  asm_level_set_contact_tangent_matrix(matl,mcb,scb,LM,active_contact_region);
675  }
676 
677  //Assemble RHS
678  if (version & model::BUILD_RHS ) {
679  GMM_TRACE2("Level set contact brick RHS assembly on "
680  << mcb.get_pair_info(scb.get_var_name()).num_of_integr_elems()
681  << " elements");
682  asm_level_set_contact_rhs(vecl,mcb,scb,LM,active_contact_region);
683  for(size_type i=0;i<vecl.size();i++)
684  gmm::scale(vecl[i], scalar_type(-1));
685  }
686 }
687 
688 
689 void level_set_contact::NormalTerm::compute(
691  bgeot::base_tensor &t)
692 {
693  size_type cv = ctx.convex_num();
694  size_type cv_volume = mcb.ext_face_of_elem(cv).cv;
695  bgeot::short_type f_volume = mcb.ext_face_of_elem(cv).f;
696  bgeot::base_node un = mcb.get_mesh().normal_of_face_of_convex
697  (cv_volume, bgeot::short_type(f_volume), ctx.xref());
698  un /= gmm::vect_norm2(un);
699 
700  if (version == 1) {
701  for (size_type i = 0; i < dim; i++) t[i] = un[i];
702  } else {
703  for (size_type i = 0; i < dim; i++)
704  for (size_type j = 0; j < dim; j++)
705  if (i == j) t(i, j) = 1.0 - un[i] * un[j];
706  else t(i, j) =-un[i] * un[j];
707  }
708 }
709 
710 
711 level_set_contact::HFunction::HFunction(
712  const mesh_fem &lsmf_,
713  const plain_vector &LS_U_,
714  scalar_type epsilon,
715  scalar_type small_h_):
716 
717 lsmf(lsmf_),
718  LS_U(LS_U_),
719  m_Epsilon(epsilon),
720  small_h(small_h_),
721  sizes_(1)
722 {sizes_[0]=1;}
723 
724 const bgeot::multi_index& level_set_contact::HFunction::
725  sizes(size_type) const { return sizes_;}
726 
727 void level_set_contact::HFunction::
728 prepare(getfem::fem_interpolation_context& /*ctx*/, size_type /*nl_part*/) {}
729 
730 void level_set_contact::HFunction::compute(getfem::fem_interpolation_context& ctx,
731  bgeot::base_tensor &t)
732 {
733  size_type cv = ctx.convex_num();
734  plain_vector U;
735  slice_vector_on_basic_dof_of_element(lsmf, LS_U, cv, U);
736  // plain_vector U(lsmf.nb_basic_dof_of_element(cv));
737  // gmm::copy(gmm::sub_vector(LS_U,gmm::sub_index(lsmf.ind_basic_dof_of_element(cv))),U);
738  plain_vector ls_interpolated(1);
739  ctx.pf()->interpolation(ctx,U,ls_interpolated,1);
740  t[0] = hRegularized(ls_interpolated[0],m_Epsilon,small_h);
741 }
742 
743 bgeot::scalar_type level_set_contact::HFunction::
744  hRegularized(scalar_type f, scalar_type epsilon, scalar_type small_h_)
745 {
746  if (f>epsilon) return 1.0;
747  if (f<(-epsilon)) return small_h_;
748  return 0.5+0.125*(9.0*f/(epsilon)-5.0*pow(f/(epsilon),scalar_type(3)));
749 }
750 
751 
752 level_set_contact::Unity::Unity(const mesh_fem &mf_):mf(mf_),sizes_(1)
753 {sizes_[0]=1;}
754 const bgeot::multi_index& level_set_contact::Unity::sizes(size_type) const {return sizes_;}
755 void level_set_contact::Unity::
756 prepare(getfem::fem_interpolation_context& /*ctx*/, size_type /*nl_part*/) {}
757 void level_set_contact::Unity::
758 compute(getfem::fem_interpolation_context& /*ctx*/, bgeot::base_tensor &t){t[0]=1.0;}
759 
760 
761 void level_set_contact::
762  solve_with_contact(
763  SOLVE_FUNCTION sf,
764  getfem::model& md,
765  gmm::iteration& it_newton,
766  gmm::iteration& it_staggered,
767  const std::string& lsolver_name,
768  getfem::abstract_newton_line_search &ls)
769 {
770  bool active_set_converged = false;
771  it_staggered.set_iteration(0);
773 
774  do {
775  active_set_converged = !master_contact_body::any_contact_change();
776  it_newton.set_iteration(0);
777  getfem::rmodel_plsolver_type plsolver=getfem::select_linear_solver<sparse_matrix,plain_vector>(md,lsolver_name);
778  (*sf)(md,it_newton,plsolver,ls);
779  GMM_TRACE2("Newton converged? - "<<it_newton.converged());
780  GMM_TRACE2("active set converged? - "<<active_set_converged);
781  GMM_ASSERT1(it_newton.converged(),"Newton method did not converge");
782  it_staggered++;
783 
784  } while(!active_set_converged &&
785  !it_staggered.finished(it_staggered.get_resmax()+1.0));
786 
787  if (active_set_converged && it_newton.converged()) it_staggered.enforce_converged();
788 }
789 
792  bgeot::pgeometric_trans pelem_trans)
793 {
794  std::string name = bgeot::name_of_geometric_trans(pelem_trans);
795  std::stringstream fname;
796  fname<<name.substr(0,5);
797  GMM_ASSERT1((fname.str()=="GT_QK" || fname.str()=="GT_PK"),
798  "Cannot handle other transformations but QK or PK,\
799  Sorry, to be implemented" );
800  std::stringstream str1(name.substr(6,1));
801  size_type dim;
802  str1>>dim;
803 
804  std::istringstream str2(name.substr(8,1));
805  size_type order;
806  str2>>order;
807 
808  fname<<"("<<dim-1<<","<<order<<")";
809  return bgeot::geometric_trans_descriptor(fname.str());
810 }
getfem::partial_mesh_fem::adapt
void adapt(const dal::bit_vector &kept_dof, const dal::bit_vector &rejected_elt=dal::bit_vector())
build the mesh_fem keeping only the dof of the original mesh_fem which are listed in kept_dof.
Definition: getfem_partial_mesh_fem.cc:41
getfem::partial_mesh_fem
a subclass of mesh_fem which allows to eliminate a number of dof of the original mesh_fem.
Definition: getfem_partial_mesh_fem.h:55
level_set_contact::face_trans_of_elem
bgeot::pgeometric_trans face_trans_of_elem(bgeot::pgeometric_trans pelem_trans)
Determines geometric transformation on the face of the element based on the geometric transformation ...
Definition: getfem_level_set_contact.cc:791
level_set_contact::slave_contact_body
Contact body that will be projected on the boundary of the master.
Definition: getfem_level_set_contact.h:144
level_set_contact::contact_pair_info
Prepares the final information needed to pass to the contact brick for every contact pair to assemble...
Definition: getfem_level_set_contact.h:178
level_set_contact::contact_pair_update
temporary object that updates contact pair, deformes meshes and undeformes when it selfdestructs
Definition: getfem_level_set_contact.h:437
bgeot::name_of_geometric_trans
std::string name_of_geometric_trans(pgeometric_trans p)
Get the string name of a geometric transformation.
Definition: bgeot_geometric_trans.cc:1173
level_set_contact::master_contact_body::get_pair_info
const contact_pair_info & get_pair_info(const std::string &slave_var_name) const
access to a structure that contains all the info about contact pair between this master and a slave,...
Definition: getfem_level_set_contact.cc:132
level_set_contact::master_contact_body::master_contact_body
master_contact_body(model &_md, const std::string &_var_name, size_type _mult_order, size_type _mult_mim_order)
create master contact body with a model, name where masters displacements are defined,...
Definition: getfem_level_set_contact.cc:71
getfem::mesh::region
const mesh_region region(size_type id) const
Return the region of index 'id'.
Definition: getfem_mesh.h:421
level_set_contact::master_contact_body::master_contact_changed
bool master_contact_changed(void)
contact detection for all slaves
Definition: getfem_level_set_contact.cc:263
level_set_contact::master_contact_body::clear_contact_history
void clear_contact_history(void)
clearing previous contact elem lists
Definition: getfem_level_set_contact.cc:255
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::compute_gradient
void compute_gradient(const mesh_fem &mf, const mesh_fem &mf_target, const VECT1 &UU, VECT2 &VV)
Compute the gradient of a field on a getfem::mesh_fem.
Definition: getfem_derivatives.h:62
getfem::mesh_fem::set_finite_element
void set_finite_element(size_type cv, pfem pf)
Set the finite element method of a convex.
Definition: getfem_mesh_fem.cc:127
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::partial_mesh_fem::nb_dof
size_type nb_dof(void) const
Return the total number of degrees of freedom.
Definition: getfem_partial_mesh_fem.h:133
level_set_contact::slave_contact_body::slave_contact_body
slave_contact_body(model &_md, const std::string &_var_name, mesh_im *_pmim)
default constructor.
Definition: getfem_level_set_contact.cc:42
getfem::mesh_fem::dof_on_region
dal::bit_vector dof_on_region(const mesh_region &b) const
Get a list of dof lying on a given mesh_region.
Definition: getfem_mesh_fem.cc:114
getfem_mesh_level_set.h
Keep informations about a mesh crossed by level-sets.
getfem::model::add_initialized_fem_data
void add_initialized_fem_data(const std::string &name, const mesh_fem &mf, const VECT &v)
Add an initialized fixed size data to the model, assumed to be a vector field if the size of the vect...
Definition: getfem_models.h:835
bgeot::geometric_trans_descriptor
pgeometric_trans geometric_trans_descriptor(std::string name)
Get the geometric transformation from its string name.
Definition: bgeot_geometric_trans.cc:1163
level_set_contact::master_contact_body::add_slave
void add_slave(slave_contact_body &scb, size_type slave_contact_region=-1)
associate a slave contact body with this master.
Definition: getfem_level_set_contact.cc:162
bgeot::geotrans_interpolation_context::xref
const base_node & xref() const
coordinates of the current point, in the reference convex.
Definition: bgeot_geometric_trans.cc:291
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::fem_interpolation_context::convex_num
size_type convex_num() const
get the current convex number
Definition: getfem_fem.cc:52
getfem::mesh_region::free_region_id
static size_type free_region_id(const getfem::mesh &m)
Extract the next region number that does not yet exists in the mesh.
Definition: getfem_mesh_region.cc:480
getfem::interpolation
void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target, const VECTU &U, VECTV &V, int extrapolation=0, double EPS=1E-10, mesh_region rg_source=mesh_region::all_convexes(), mesh_region rg_target=mesh_region::all_convexes())
interpolation/extrapolation of (mf_source, U) on mf_target.
Definition: getfem_interpolation.h:693
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
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
level_set_contact::level_set_contact_brick::asm_real_tangent_terms
virtual void asm_real_tangent_terms(const model &md, size_type, const model::varnamelist &vl, const model::varnamelist &dl, const model::mimlist &mims, model::real_matlist &matl, model::real_veclist &vecl, model::real_veclist &, size_type region, build_version version) const
Assembly of bricks real tangent terms.
Definition: getfem_level_set_contact.cc:632
level_set_contact::master_contact_body::any_contact_change
static bool any_contact_change()
contact detection for all masters/slave couples
Definition: getfem_level_set_contact.cc:235
getfem::mesh_fem::nb_basic_dof
virtual size_type nb_basic_dof() const
Return the total number of basic degrees of freedom (before the optional reduction).
Definition: getfem_mesh_fem.h:556
getfem_interpolated_fem.h
FEM which interpolates a mesh_fem on a different mesh.
getfem::del_interpolated_fem
void del_interpolated_fem(const pfem &pf)
release an interpolated fem
Definition: getfem_interpolated_fem.h:187
level_set_contact::contact_pair_info::clear_contact_history
void clear_contact_history()
clearing contact element lists
Definition: getfem_level_set_contact.cc:359
getfem::mesh_region::index
const dal::bit_vector & index() const
Index of the region convexes, or the convexes from the partition on the current thread.
Definition: getfem_mesh_region.cc:237
level_set_contact::contact_pair_info::contact_changed
bool contact_changed()
Actual master/slave contact detection.
Definition: getfem_level_set_contact.cc:365
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
bgeot::small_vector< scalar_type >
getfem_derivatives.h
Compute the gradient of a field on a getfem::mesh_fem.
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
level_set_contact::master_contact_body::build_mesh_im_on_boundary
std::shared_ptr< mesh_im > build_mesh_im_on_boundary(size_type region)
return a pointer to mesh_im used for contact surface calculations
Definition: getfem_level_set_contact.cc:276
level_set_contact::contact_body
base class for the master and the slave contact bodies.
Definition: getfem_level_set_contact.h:119
level_set_contact::contact_pair_info::update
void update(void) const
updating contact information (mesh_fem's, mesh_im's) with the deformation.
Definition: getfem_level_set_contact.cc:508
getfem::new_interpolated_fem
pfem new_interpolated_fem(const mesh_fem &mef, const mesh_im &mim, pinterpolated_func pif=0, dal::bit_vector blocked_dof=dal::bit_vector(), bool store_val=true)
create a new interpolated FEM.
Definition: getfem_interpolated_fem.cc:394
getfem_level_set_contact.h
Non frictional level set based large sliding contact; for details see: A. Andreykiv et al....
level_set_contact::slave_contact_body::offset_level_set
void offset_level_set(scalar_type off)
adds a fixed value "off" to the level set field
Definition: getfem_level_set_contact.cc:65
getfem_level_set.h
Define level-sets.
level_set_contact::boundary_level_set_field
void boundary_level_set_field(const getfem::mesh &_mesh, const getfem::mesh_fem &mf, const getfem::mesh_im &mim, VECT &LS)
build a level set function on mesh with zero on the boundary.
Definition: getfem_level_set_contact.h:72
getfem::slice_vector_on_basic_dof_of_element
void slice_vector_on_basic_dof_of_element(const mesh_fem &mf, const VEC1 &vec, size_type cv, VEC2 &coeff, size_type qmult1=size_type(-1), size_type qmult2=size_type(-1))
Given a mesh_fem.
Definition: getfem_mesh_fem.h:659
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::mesh_fem::ind_basic_dof_of_element
virtual ind_dof_ct ind_basic_dof_of_element(size_type cv) const
Give an array of the dof numbers a of convex.
Definition: getfem_mesh_fem.h:450
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem::fem_interpolation_context::pf
const pfem pf() const
get the current FEM descriptor
Definition: getfem_fem.h:782
level_set_contact::master_contact_body::ext_face_of_elem
face_type ext_face_of_elem(size_type cv) const
gives a face, corresponding to newly created boundary element
Definition: getfem_level_set_contact.cc:153
getfem_mesh_im_level_set.h
a subclass of mesh_im which is conformal to a number of level sets.
getfem::outer_faces_of_mesh
void APIDECL outer_faces_of_mesh(const mesh &m, const dal::bit_vector &cvlst, convex_face_ct &flist)
returns a list of "exterior" faces of a mesh (i.e.
Definition: getfem_mesh.cc:822
level_set_contact::master_contact_body
Master contact body which surface will be used to project contact stresses and stiffness terms.
Definition: getfem_level_set_contact.h:287
getfem::pbrick
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
Definition: getfem_models.h:49
getfem::mesh_fem::set_classical_discontinuous_finite_element
void set_classical_discontinuous_finite_element(size_type cv, dim_type fem_degree, scalar_type alpha=0, bool complete=false)
Similar to set_classical_finite_element, but uses discontinuous lagrange elements.
Definition: getfem_mesh_fem.cc:200
getfem::mesh_fem::nb_dof
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
Definition: getfem_mesh_fem.h:562
getfem::mesh_fem::set_qdim
virtual void set_qdim(dim_type q)
Change the Q dimension.
Definition: getfem_mesh_fem.h:316
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:2959
getfem::classical_mesh_fem
const mesh_fem & classical_mesh_fem(const mesh &mesh, dim_type degree, dim_type qdim=1, bool complete=false)
Gives the descriptor of a classical finite element method of degree K on mesh.
Definition: getfem_mesh_fem.cc:862
level_set_contact::master_contact_body::clear_all_contact_history
static void clear_all_contact_history()
should be used in the beginning of a step to clean data structures that store previous contact elemen...
Definition: getfem_level_set_contact.cc:247
level_set_contact::add_level_set_normal_contact_brick
size_type add_level_set_normal_contact_brick(model &md, master_contact_body &mcb, slave_contact_body &scb, size_type rg=-1)
adding level set based normal contact brick to the model.
Definition: getfem_level_set_contact.cc:565