GetFEM  5.4.2
getfem_level_set_contact.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2012-2020 Andriy Andreykiv
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /** @file getfem_level_set_contact.h
33  @author "Andriy Andreykiv" <andriy.andreykiv@gmail.com>
34  @date July, 2012.
35  @brief Non frictional level set based large sliding contact;
36  for details see:
37  A. Andreykiv et al. A level set based large sliding contact
38  algorithm for an easy analysis of implant positioning
39  2012 International Journal for Numerical Methods in Engineering,
40  89, pp. 1317-1336
41  2D and 3D Examples of the usage: test_contact.cpp
42  */
43 
44 
45 #pragma once
46 
47 #ifndef GETFEM_LEVEL_SET_CONTACT_H__
48 #define GETFEM_LEVEL_SET_CONTACT_H__
49 
51 #include <getfem/getfem_models.h>
54 #include <gmm/gmm_except.h>
55 
56 namespace level_set_contact {
57 
58  using getfem::mesh_fem;
59  using getfem::mesh_im;
60  using getfem::mesh;
61  using getfem::model;
62  using getfem::size_type;
63  using getfem::scalar_type;
64  using getfem::model_real_plain_vector;
65  typedef getfem::model_real_plain_vector plain_vector;
66  typedef getfem::model_real_sparse_matrix sparse_matrix;
67 
68 
69  /**build a level set function on mesh with zero on the boundary.
70  Solves Laplace equation with zero Dirichlet on the boundary.
71  Used to create simple level sets for contact testing*/
72  template<class VECT> void boundary_level_set_field(
73  const getfem::mesh& _mesh,
74  const getfem::mesh_fem& mf,
75  const getfem::mesh_im& mim,
76  VECT& LS)
77  {
78  getfem::mesh& mesh = const_cast<getfem::mesh&>(_mesh);
79  //model and vars
80  getfem::model md;
81  md.add_fem_variable("LS",mf);
82  getfem::model_real_plain_vector RHS(mf.nb_dof());
83 
84  //calculating the size of the LS RHS based on the size of the geometry
85  getfem::base_node Pmin(mesh.dim()),Pmax(mesh.dim()),range(mesh.dim());
86  mesh.bounding_box(Pmin,Pmax);
87  gmm::add(Pmax,gmm::scaled(Pmin,-1.0),range);
88  getfem::scalar_type mesh_size = *(std::max_element(range.begin(),range.end()));
89  gmm::fill(RHS,mesh_size*5.0);
90  md.add_initialized_fem_data("RHS",mf,RHS);
91 
92  //border region
93  getfem::mesh_region border_faces;
94  getfem::outer_faces_of_mesh(mesh, border_faces);
96  for (getfem::mr_visitor i(border_faces); !i.finished(); ++i)
97  mesh.region(BORDER).add(i.cv(),i.f());
98 
99  //describing the PDE problem
100  getfem::add_Laplacian_brick(md,mim,"LS");
102  getfem::add_source_term_brick(md,mim,"LS","RHS");
103 
104  //solving
105  gmm::iteration iter;
106  GMM_TRACE2("building scalar level set with laplace equation..");
107  getfem::standard_solve(md,iter);
108 
109  //extracting the result
110  gmm::copy(md.real_variable("LS"),LS);
111 
112  //so, now the mesh is as it was, hence const is still valid
113  mesh.sup_region(BORDER);
114  GMM_TRACE2("..done");
115  }
116 
117 
118  /**base class for the master and the slave contact bodies.*/
120 
121  const std::string var_name;
122  bool is_deformed;
123  friend class contact_pair_update;
124 
125  protected:
126  mesh& own_mesh;
127  const mesh_fem& own_mesh_fem;
128  model& md;
129 
130  public:
131 
132  contact_body(model& _md, std::string _var_name);
133  inline std::string get_var_name() const {return var_name;}
134  inline mesh& get_mesh() {return own_mesh;}
135  inline const mesh& get_mesh()const {return own_mesh;}
136  inline const mesh_fem& get_mesh_fem() const {return own_mesh_fem;}
137  inline const model& get_model() const {return md;}
138  inline bool is_mesh_deformed() const {return is_deformed;}
139  };
140 
141 
142  /** Contact body that will be projected on the boundary
143  of the master. */
145 
146  std::string ls_name;
147  mesh_fem ls_mesh_fem;
148  mesh_im* pmim;
149 
150  public:
151 
152  /**default constructor. Level set field will have zero value
153  right on the boundary of the contact body*/
154  slave_contact_body(model& _md, const std::string& _var_name,
155  mesh_im* _pmim);
156 
157  /**Level set field is provided via the model variable name*/
158  slave_contact_body(model& _md, std::string _var_name,
159  std::string _ls_name);
160  inline std::string get_ls_name() const {return ls_name;}
161  inline const plain_vector& ls_values() const
162  {return md.real_variable(ls_name);}
163  inline plain_vector& ls_values()
164  {return md.set_real_variable(ls_name);}
165  inline const mesh_fem& get_ls_mesh_fem() const {return md.mesh_fem_of_variable(ls_name);}
166  template<class VECTOR> void set_level_set(const VECTOR& ls)
167  {gmm::copy(ls,md.set_real_variable(ls_name));}
168 
169  /**adds a fixed value "off" to the level set field */
170  void offset_level_set(scalar_type off);
171  };
172 
173 
174  class master_contact_body;
175 
176  /**Prepares the final information needed to pass to the contact
177  brick for every contact pair to assemble tangent terms*/
179 
180  //obtain on construction
181  master_contact_body& master_cb;
182  slave_contact_body& slave_cb;
183  const std::string mult_name;
184  const size_type GIVEN_CONTACT_REGION;
185 
186  //to be built
187  dal::bit_vector old_contact_elm_list;
188  dal::bit_vector pre_old_ct_list;
189  size_type ACTIVE_CONTACT_REGION;
190  mutable std::shared_ptr<mesh_im> pmim_contact;
191  mutable getfem::pfem ifem_srf;
192  mutable std::shared_ptr<mesh_fem> pinterpolated_fem;
193  mutable std::shared_ptr<mesh_fem> pinterpolated_fem_U;
194  mutable std::shared_ptr<gmm::unsorted_sub_index> slave_ls_dofs;
195  mutable std::shared_ptr<gmm::unsorted_sub_index> slave_U_dofs;
196  mutable size_type n_integrated_elems;
197 
198  // state of the object
199  mutable bool members_are_computed;
200  mutable bool init_cont_detect_done;
201  public:
202 
203  // accessors
204  inline const mesh_fem& slave_scalar_fem() const {
205  if (dependecies_changed()) update();
206  return *pinterpolated_fem;
207  }
208  inline const mesh_fem& slave_vector_fem() const {
209  if (dependecies_changed()) update();
210  return *pinterpolated_fem_U;
211  }
212  inline const gmm::unsorted_sub_index& slave_scalar_dofs() const {
213  if (dependecies_changed()) update();
214  return *slave_ls_dofs;
215  }
216  inline const gmm::unsorted_sub_index& slave_vector_dofs() const {
217  if (dependecies_changed()) update();
218  return *slave_U_dofs;
219  }
220  inline const mesh_im& contact_mesh_im() const {
221  if (dependecies_changed()) update();
222  return *pmim_contact;
223  }
224 
225  inline size_type contact_region() const
226  {return ACTIVE_CONTACT_REGION;}
227 
228  inline const std::string& get_mult_name() const
229  {return mult_name;}
230 
231  inline size_type num_of_integr_elems() const {return n_integrated_elems;}
232  // update
233  inline bool dependecies_changed() const
234  {return !members_are_computed;}
235  inline void force_update() const
236  {members_are_computed=false;}
237 
238  /** Actual master/slave contact detection. Level set field is projected on the
239  boundary of the master and only the elements which nodes satisfy
240  level_set + Multiplier > 0
241  become contact elements*/
242  bool contact_changed();
243 
244  /**clearing contact element lists*/
245  void clear_contact_history();
246 
247  /** updating contact information (mesh_fem's, mesh_im's)
248  with the deformation. Contact detection is not performed*/
249  void update(void) const;
250 
251  contact_pair_info(master_contact_body& underformed_mcb,
252  slave_contact_body& underformed_scb, const std::string& _mult_name,
253  size_type _GIVEN_CONTACT_REGION);
254 
255  private:
256  /**prohibiting copying*/
258  contact_pair_info& operator=(const contact_pair_info&);
259 
260 
261  };
262 
263  struct face_type{
265  face_type(size_type _cv=0, bgeot::short_type _f=0):cv(_cv),f(_f){}
266  face_type(const getfem::mr_visitor& i): cv(i.cv()),f(i.f()){}
267  };
268 
269  /**Determines geometric transformation on the face of the element
270  based on the geometric transformation of the element itself. Works
271  only for PK and QK elements*/
273 
274 
275  /** Master contact body which surface will be used to project contact
276  stresses and stiffness terms. It contains and manages the slaves and
277  knows other masters.
278  Master contact body must be created with mesh_fem that allows automatic
279  addition of mesh_fem description on new elements (use mesh_fem::set_auto_add or
280  use set_classical_finite_element). This feature is used when new boundary
281  elements are created from faces. At the same time the mesh_im object that
282  is used to add for instance some structural bricks on the volume (elastostatic,
283  nonlinear_elastostatic, updated_lagrangian) should be either created before
284  master contact body, or set on master_contact_body::volume_region() if it's
285  created after. This is to avoid integration of the volume integrals on the
286  boundary elemenents of lower dimension. */
288 
289 
290  const size_type mult_mim_order;
291  const std::string mult_int_method;
292  size_type BOUNDARY_ELEMENTS, VOLUME_ELEMENTS;
293  std::vector<size_type> face_to_belem_ind;
294  static std::vector<master_contact_body*> masters;
295  std::map<std::string, std::shared_ptr<contact_pair_info> > contact_table;
296  std::map<size_type,face_type> border_faces;
297 
298  protected:
299 
300  /**contact detection for all slaves*/
301  bool master_contact_changed(void);
302 
303  /** clearing previous contact elem lists*/
304  void clear_contact_history(void);
305 
306  public:
307 
308  enum contact_integration{PER_ELEMENT=1,REGULARIZED_LEVEL_SET=2};
309 
310  /**approximation order for Lagrange multiplier on the contact surface*/
312 
313  /**integration approach for contact elements that are partially
314  crossed by level sets:
315  PER_ELEMENT - a whole element is incuded into contact (default)
316  REGULARIZED_LEVEL_SET - Gauss points where projected value of the level
317  set is < zero are set to zero or small value
318  (with gradual transition)*/
319  const contact_integration integration;
320 
321  /**width of transition for a regularazied Heaviside function in
322  case of REGULARIZED_LEVEL_SET*/
323  const scalar_type regularized_tollerance;
324 
325  /**in case of REGULARIZED_LEVEL_SET this value
326  scales weight of Gauss points that have negative level
327  set value*/
328  const scalar_type small_weight_multiplier;
329 
330  /**if the angle (in degrees) between contact element and
331  level set contour exceed this value, this element is not included in
332  contact algorithm*/
333  const scalar_type max_contact_angle;
334 
335 
336  /** create master contact body with a model,
337  name where masters displacements are defined, order for
338  Lagrange multiplier, order for the integration method*/
340  const std::string& _var_name,
341  size_type _mult_order, size_type _mult_mim_order);
342 
343  /**the same as above, but specifically provide itegration
344  * method on the contact surface (_mult_int_method), additionally,
345  * specify if surface contact elements have to be cut by the level set.
346  * The later ensures that contact surface is strictly a domain
347  * that overlaps with the slave, hence this allows smooth growth of the contact
348  * surface. The level set cutting is done using regularized Heaviside function
349  */
351  const std::string& _var_name,
352  size_type _mult_order,
353  const std::string& _mult_int_method,
354  contact_integration _integration = PER_ELEMENT,
355  scalar_type _regularized_tollerance = 1e-6,
356  scalar_type _small_weight_multiplier = 0.001,
357  scalar_type _max_contact_angle = 45);
358 
359  /** associate a slave contact body with this master. \
360  specify a region of expected contact interaction. \
361  (takes the whole master boundary if not specified)*/
362  void add_slave(slave_contact_body& scb,
363  size_type slave_contact_region = -1);
364 
365  /** order of integration of boundary contact terms*/
366  inline size_type contact_mim_order() const
367  {
368  GMM_ASSERT1(mult_mim_order!=size_type(-1),
369  "master body was not created with "
370  "order of integration for contact area");
371  return mult_mim_order;
372  }
373 
374  /** integration method on the contact surface,
375  * use it when the master is created with a specific
376  * integration method and not the approx_order*/
377  inline getfem::pintegration_method contact_int_method() const
378  {
379  GMM_ASSERT1(mult_mim_order==size_type(-1),
380  "master body was not created with integration "
381  "method for contact area");
382  return getfem::int_method_descriptor(mult_int_method);
383  }
384 
385  /** region of all volume elements without the boundary*/
386  inline size_type volume_region() const
387  {return VOLUME_ELEMENTS;}
388 
389  /**boundary elements, added after creation of
390  the master contact body */
391  inline size_type boundary_region() const
392  {return BOUNDARY_ELEMENTS;}
393 
394  /**access to a structure that contains all the info
395  about contact pair between this master and a slave, defined
396  on @param slave_var_name*/
397  const contact_pair_info& get_pair_info(const std::string& slave_var_name) const;
398 
399  /**the same as above, but non-const*/
400  contact_pair_info& get_pair_info(const std::string& slave_var_name);
401 
402 
403  /**contact detection for all masters/slave couples
404  @return true if any of the contact areas where changed
405  (which requires new Newton method run)*/
406  static bool any_contact_change();
407 
408  /** should be used in the beginning of a step
409  to clean data structures that store previous
410  contact element lists (used to verify if contact surface
411  is converged to one list)
412  */
413  static void clear_all_contact_history();
414 
415  inline void update_for_slave(std::string slave_var_name)
416  { contact_table[slave_var_name]->update(); }
417 
418  /** return a pointer to mesh_im used for contact surface calculations
419  */
420  std::shared_ptr<mesh_im> build_mesh_im_on_boundary(size_type region);
421 
422  /**gives a face, corresponding to newly created
423  boundary element @param cv*/
424  face_type ext_face_of_elem(size_type cv) const;
425 
426  private:
427  /**prohibiting copying*/
429  master_contact_body& operator=(const master_contact_body&);
430 
431  };
432 
433  enum update_depth{DEFORM_MESHES_ONLY,FULL_UPDATE};
434 
435  /**temporary object that updates contact pair,
436  deformes meshes and undeformes when it selfdestructs*/
438  std::shared_ptr<getfem::temporary_mesh_deformator<> > def_master;
439  std::shared_ptr<getfem::temporary_mesh_deformator<> > def_slave;
440  master_contact_body& mcb;
441  slave_contact_body& scb;
442  public:
444  slave_contact_body& _scb,
445  update_depth ud = FULL_UPDATE);
446 
448  };
449 
450 
451  /** adding level set based normal contact brick to the model.
452  The contact is etablished between the
453  @param mcb - master contact body and
454  @param scb - slave contact body, defined on
455  @param md - model object
456  @param rg - optional assumed contact region
457  helping to narrow down contact search
458  Note, this contact algorithm is note stabilized, hence,
459  master contact body mesh should be coarser than slave's mesh.
460  Otherwise this contact constraint will violate inf-sub condition
461  and the solver will fail (or diverge, if it's iterative)
462  */
464  master_contact_body& mcb,
465  slave_contact_body& scb,
466  size_type rg = -1);
467 
468 
469  /** assembles normal contact terms on the boundary of
470  two contact bodies (master/slave)*/
472 
473  model& md;
474  master_contact_body& mcb;
475  slave_contact_body& scb;
476 
477  /**id of the region of faces where contact has to be checked*/
478  size_type given_contact_id;
479 
480  /**id of the region of boundary elements,
481  corresponding to the above faces*/
482  size_type contact_region_id;
483 
484  /**actual region object, with id = contact_region_id*/
485  getfem::mesh_region contact_region;
486 
487  public:
488  virtual void asm_real_tangent_terms(
489  const model &md, size_type /* ib */,
490  const model::varnamelist &vl,
491  const model::varnamelist &dl,
492  const model::mimlist &mims,
493  model::real_matlist &matl,
494  model::real_veclist &vecl,
495  model::real_veclist &,
496  size_type region,
497  build_version version) const;
498 
500  model& _md,
501  master_contact_body& _mcb,
502  slave_contact_body& scb,
503  size_type rg = -1);
504  };
505 
506 
507  /** A term, used in level set contact assemblies that
508  builds a surface projection matrix R = N^t X N
509  (where N is normal vector to the boundary)
510  */
512  { private:
513  const master_contact_body& mcb;
514  bgeot::multi_index sizes_;
515  bgeot::size_type version;
516  bgeot::size_type dim;
517 
518  public:
519 
520  NormalTerm(const master_contact_body& _mcb, size_type version_ = 1) :
521  mcb(_mcb),
522  sizes_(version_),
523  version(version_),
524  dim(_mcb.get_mesh().dim()) {
525 
526  GMM_ASSERT1(dim==2 || dim==3, "NormalTerm: wrong space dimension ");
527  GMM_ASSERT1(version==1 || version==2,"NormalTerm:: wrong version ");
528 
529  if (version == 1)
530  if (dim == 2)
531  sizes_[0] = 2;
532  else
533  sizes_[0] = 3;
534  else
535  if (dim == 2) {
536  sizes_[0] = 2;
537  sizes_[1] = 2;
538  }
539  else {
540  sizes_[0] = 3;
541  sizes_[1] = 3;
542  }
543  }
544  const bgeot::multi_index &sizes(size_type) const {return sizes_;};
545  void compute(getfem::fem_interpolation_context& ctx, bgeot::base_tensor &t);
546  void prepare(getfem::fem_interpolation_context& /* ctx */, size_type /* nl_part */) {}
547 
548  };
549 
550  /** Regularized Heaviside function.
551  Can be used instead of mesh_im_level_set in assemblies.
552  It's more stable, as it never fails in comparison to Delauney method
553  (used inside mesh_im_level_set), but less accurate, as it has a
554  transition zone from 1 to 0 of epsilon width.
555  The idea is taken from one of the articles of Ted Belytschko on XFem*/
557  {
558  private:
559  const mesh_fem &lsmf;
560  const plain_vector &LS_U;
561  scalar_type m_Epsilon;
562  scalar_type small_h;
563  bgeot::multi_index sizes_;
564 
565 
566  public:
567  HFunction(
568  const mesh_fem &lsmf_,
569  const plain_vector &LS_U_,
570  scalar_type epsilon=1e-9,
571  scalar_type small_h_=0);
572  const bgeot::multi_index &sizes(size_type) const;
573  void prepare(getfem::fem_interpolation_context& ctx, size_type nl_part);
574  void compute(getfem::fem_interpolation_context& ctx, bgeot::base_tensor &t);
575  scalar_type hRegularized(scalar_type x, scalar_type epsion, scalar_type small);
576  };
577 
578  //A dummy nonlinear term, does nothing
579  class Unity : public getfem::nonlinear_elem_term
580  {
581  private:
582  const mesh_fem &mf;
583  bgeot::multi_index sizes_;
584 
585  public:
586  Unity(const mesh_fem &mf_);
587  const bgeot::multi_index &sizes(size_type) const;
588  void prepare(getfem::fem_interpolation_context& ctx, size_type nl_part);
589  void compute(getfem::fem_interpolation_context& ctx, bgeot::base_tensor &t);
590  };
591 
592 
593 
594  template<typename MAT, typename VECT>
595  void asm_level_set_contact_tangent_matrix(
596  std::vector<MAT>& matl,
597  const master_contact_body& mcb,
598  const slave_contact_body& scb,
599  const VECT& LM,
600  const getfem::mesh_region& contact_region)
601  {
602  //extract matrix references
603  MAT& Kmm = matl[0];
604  MAT& Kss = matl[1];
605  //MAT& Kll = matl[2] remains zero
606  MAT& Kms = matl[3];
607  MAT& Kml = matl[4];
608  MAT& Ksl = matl[5];
609 
610  const std::string& mult_name =
611  mcb.get_pair_info(scb.get_var_name()).get_mult_name();
612  const std::string ls_name = "ls_on"+mcb.get_var_name()+"_from_"+scb.get_var_name();
613 
614  //extract mfs, and mims
615  const mesh_fem& mf_U_line = mcb.get_mesh_fem();
616  const mesh_fem& mf_lambda = mcb.get_model().mesh_fem_of_variable(mult_name);
617  const mesh_fem& mf_interpolate =
618  mcb.get_pair_info(scb.get_var_name()).slave_scalar_fem();
619  const mesh_fem& mf_U_interpolate =
620  mcb.get_pair_info(scb.get_var_name()).slave_vector_fem();
621  const mesh_fem& mf_master_ls = mcb.get_model().mesh_fem_of_variable(ls_name);
622  const mesh_im& mim_line =
623  mcb.get_pair_info(scb.get_var_name()).contact_mesh_im();
624 
625  //build temp vectors for interpolated fems
626  plain_vector LS_small(mf_interpolate.nb_dof());
627  gmm::copy(gmm::sub_vector(scb.ls_values(),
628  mcb.get_pair_info(scb.get_var_name()).slave_scalar_dofs()),LS_small);
629 
630  //nonlinear term to compute normal vector and R matrix
631  NormalTerm R_matrix(mcb,2);
632 
633  //nonlinear term that describes regularized integration or dummy (unity) multiplier
634  std::shared_ptr<getfem::nonlinear_elem_term> integration;
635  if (mcb.integration==master_contact_body::REGULARIZED_LEVEL_SET) {
636  integration = std::make_shared<HFunction>
637  (mf_master_ls,mcb.get_model().real_variable(ls_name),
638  mcb.regularized_tollerance,mcb.small_weight_multiplier);
639  } else { integration = std::make_shared<Unity>(mf_master_ls); }
640 
641 
642  //temp matrices due to different DOF indeces of the slave
643  sparse_matrix Kms_small(mf_U_interpolate.nb_dof(), mf_U_line.nb_dof());
644  sparse_matrix Kss_small(mf_U_interpolate.nb_dof(),mf_U_interpolate.nb_dof());
645  sparse_matrix Ksl_small(mf_U_interpolate.nb_dof(),mf_lambda.nb_dof());
646 
647  //assembly
648  getfem::generic_assembly assem_boundary;
649 
650  assem_boundary.set(
651  "F=data$1(#3);"
652  "L=data$2(#1);"
653  "Kmm1 = comp(Base(#1).Grad(#3).vBase(#2).NonLin$1(#2).vGrad(#2).NonLin$2(#5))(i,j,k,:,k,m,n,:,n,m,1).L(i).F(j);"
654  "Kmm2 = comp(Base(#1).NonLin$1(#2).vGrad(#2).Grad(#3).vBase(#2).NonLin$2(#5))(i,m,n,:,n,m,j,k,:,k,1).L(i).F(j);"
655  "Kmm3 = comp(Base(#1).Base(#3).NonLin$1(#2).vGrad(#2).NonLin$1(#2).vGrad(#2).NonLin$2(#5))(i,j,k,l,:,l,k,m,n,:,n,m,1).L(i).F(j);"
656  "Kmm4 = (-1.0)*comp(Base(#1).Base(#3).NonLin$1(#2).vGrad(#2).vGrad(#2).NonLin$2(#5))(i,j,m,n,:,n,l,:,l,m,1).L(i).F(j);"
657  "M$1(#2,#2)+= sym(Kmm1+Kmm2+Kmm3+Kmm4);"
658  "Ksm1=(-1.0)*comp(Base(#1).Grad(#3).vBase(#4).NonLin$1(#2).vGrad(#2).NonLin$2(#5))(i,j,k,:,k,m,n,:,n,m,1).L(i).F(j);"
659  "Ksm2=(-1.0)*comp(Base(#1).Grad(#3).vGrad(#4).vBase(#2).NonLin$2(#5))(i,j,m,:,m,n,:,n,1).L(i).F(j);"
660  "M$2(#4,#2)+= Ksm1+Ksm2;"
661  "Kml1=comp(Base(#3).NonLin$1(#2).vGrad(#2).Base(#1).NonLin$2(#5))(i,m,n,:,n,m,:,1).F(i);"
662  "Kml2=comp(Grad(#3).vBase(#2).Base(#1).NonLin$2(#5))(i,j,:,j,:,1).F(i);"
663  "M$3(#2,#1)+= Kml1+Kml2;"
664  "Kss_part = comp(Base(#1).Grad(#3).vGrad(#4).vBase(#4).NonLin$2(#5))(i,j,m,:,m,n,:,n,1).L(i).F(j);"
665  "M$4(#4,#4)+=sym(Kss_part{1,2}+Kss_part{2,1});"
666  "M$5(#4,#1)+=(-1.0)*comp(Grad(#3).vBase(#4).Base(#1).NonLin$2(#5))(i,k,:,k,:,1).F(i);"
667  ); /* Here we don't compute matrices that contain Hessian of
668  the level set function, as Getfem does not compute Hessian
669  for interpolated_fem class that we use for level set function */
670  assem_boundary.push_mi(mim_line); //mim on the contact surface
671  assem_boundary.push_mf(mf_lambda); //mf 1 Lambda
672  assem_boundary.push_mf(mf_U_line); //mf 2 Umaster
673  assem_boundary.push_mf(mf_interpolate); //mf 3 LSslave
674  assem_boundary.push_mf(mf_U_interpolate);//mf 4 Uslave
675  assem_boundary.push_mf(mf_master_ls); //mf 5 ls_on_master
676  assem_boundary.push_nonlinear_term(&R_matrix); //matrix of the normal products
677  assem_boundary.push_nonlinear_term(integration.get()); //term to limit integration domain
678  assem_boundary.push_data(LS_small); // data Level set on interpolated
679  assem_boundary.push_data(LM); // data Lagrange mult values
680  assem_boundary.push_mat(Kmm); // result mat 1
681  assem_boundary.push_mat(Kms_small); // .. mat 2
682  assem_boundary.push_mat(Kml); // .. mat 3
683  assem_boundary.push_mat(Kss_small); // .. mat 4
684  assem_boundary.push_mat(Ksl_small); // .. mat 5
685  assem_boundary.assembly(contact_region);
686 
687  //transfering from interpolated mesh_fem into full slave mesh_fem mat's
688  const gmm::sub_interval& Um_dof = gmm::sub_interval(0,mf_U_line.nb_dof());
689  const gmm::unsorted_sub_index& Us_dof =
690  mcb.get_pair_info(scb.get_var_name()).slave_vector_dofs();
691  const gmm::sub_interval& LM_dof = gmm::sub_interval(0,mf_lambda.nb_dof());
692  gmm::copy(gmm::transposed(Kms_small),gmm::sub_matrix(Kms,Um_dof,Us_dof));
693  gmm::copy(Kss_small,gmm::sub_matrix(Kss,Us_dof,Us_dof));
694  gmm::copy(Ksl_small,gmm::sub_matrix(Ksl,Us_dof,LM_dof));
695 
696  }
697 
698  template<typename VECT0,typename VECT1>
699  void asm_level_set_contact_rhs(
700  std::vector<VECT0>& vecl,
701  const master_contact_body& mcb,
702  const slave_contact_body& scb,
703  const VECT1& LM,
704  const getfem::mesh_region& contact_region)
705  {
706  //extract vector references
707  VECT0& RHS_Um = vecl[0];
708  VECT0& RHS_Us = vecl[1];
709  VECT0& RHS_LM = vecl[2];
710  // vecl[3, 4 and 5] remain zero
711 
712 
713  const std::string& mult_name =
714  mcb.get_pair_info(scb.get_var_name()).get_mult_name();
715  const std::string ls_name = "ls_on"+mcb.get_var_name()+"_from_"+scb.get_var_name();
716 
717  //extract mfs, and mims
718  const mesh_fem& mf_U_line = mcb.get_mesh_fem();
719  const mesh_fem& mf_lambda =
720  mcb.get_model().mesh_fem_of_variable(mult_name);
721  const mesh_fem& mf_interpolate =
722  mcb.get_pair_info(scb.get_var_name()).slave_scalar_fem();
723  const mesh_fem& mf_U_interpolate =
724  mcb.get_pair_info(scb.get_var_name()).slave_vector_fem();
725  const mesh_fem& mf_master_ls = mcb.get_model().mesh_fem_of_variable(ls_name);
726  const mesh_im& mim_line =
727  mcb.get_pair_info(scb.get_var_name()).contact_mesh_im();
728 
729  //build temp vectors for interpolated fems
730  plain_vector LS_small(mf_interpolate.nb_dof());
731  gmm::copy(gmm::sub_vector(scb.ls_values(),
732  mcb.get_pair_info(scb.get_var_name()).slave_scalar_dofs()),LS_small);
733 
734  //nonlinear term to compute normal vector and R matrix
735  NormalTerm R_matrix(mcb,2);
736 
737  //nonlinear term that describes regularized integration or dummy (unity) multiplier
738  std::shared_ptr<getfem::nonlinear_elem_term> integration;
739  if (mcb.integration==master_contact_body::REGULARIZED_LEVEL_SET) {
740  integration = std::make_shared<HFunction>
741  (mf_master_ls,mcb.get_model().real_variable(ls_name),
742  mcb.regularized_tollerance,mcb.small_weight_multiplier);
743  } else { integration = std::make_shared<Unity>(mf_master_ls); }
744 
745  // temp RHS vector due to diff DOF indeces for mesh_fem object of the slave
746  plain_vector RHS_Us_small(mf_U_interpolate.nb_dof());
747 
748  getfem::generic_assembly assem_boundary;
749  assem_boundary.set(
750  "F=data$1(#3);"
751  "L=data$2(#1);"
752  "RHS_L_Us_1=comp(Base(#1).Base(#3).NonLin$1(#2).vGrad(#2).NonLin$2(#5))(i,j,m,n,:,n,m,1).L(i).F(j);"
753  "RHS_L_Us_2=comp(Base(#1).Grad(#3).vBase(#2).NonLin$2(#5))(i,j,k,:,k,1).L(i).F(j);"
754  "V$1(#2)+=RHS_L_Us_1+RHS_L_Us_2;"
755  "V$2(#4)+=(-1.0)*comp(Base(#1).Grad(#3).vBase(#4).NonLin$2(#5))(i,j,k,:,k,1).L(i).F(j);"
756  "V$3(#1)+=comp(Base(#1).Base(#3).NonLin$2(#5))(:,i,1).F(i);"
757  );
758  assem_boundary.push_mi(mim_line); //mim on the contact surface
759  assem_boundary.push_mf(mf_lambda); //mf 1 Lambda
760  assem_boundary.push_mf(mf_U_line); //mf 2 Umaster
761  assem_boundary.push_mf(mf_interpolate); //mf 3 LSslave
762  assem_boundary.push_mf(mf_U_interpolate);//mf 4 Uslave
763  assem_boundary.push_mf(mf_master_ls); //mf 5 ls_on_master
764  assem_boundary.push_nonlinear_term(&R_matrix); //matrix of the normal products
765  assem_boundary.push_nonlinear_term(integration.get()); //term to limit integration domain
766  assem_boundary.push_data(LS_small); // data Level set on interpolated
767  assem_boundary.push_data(LM); // data Lagrange mult values
768  assem_boundary.push_vec(RHS_Um); // result vec 1
769  assem_boundary.push_vec(RHS_Us_small); // .. vec 2
770  assem_boundary.push_vec(RHS_LM); // .. vec 3
771  assem_boundary.assembly(contact_region);
772 
773  //transfering from interpolated mesh_fem into full slave mesh_fem RHS
774  const gmm::unsorted_sub_index& Us_dof =
775  mcb.get_pair_info(scb.get_var_name()).slave_vector_dofs();
776  gmm::copy(RHS_Us_small, gmm::sub_vector(RHS_Us,Us_dof));
777 
778  }
779 
780 
781 
782  typedef void(*SOLVE_FUNCTION)(
783  getfem::model &md,
784  gmm::iteration &iter,
785  getfem::rmodel_plsolver_type solver,
786  getfem::abstract_newton_line_search &ls);
787 
788  /** Solves a model that has contact in it.
789  Function checks wheather the contact area has converged
790  @param sf - a pointer to a newton solver function,
791  can be, for instance, getfem::standard_solve
792  @param it_newton - iteration object for newton method
793  @param it_staggered - iteration object for staggered calculation
794  between conact detection and newton method (only max
795  num. of iterations should be provided)
796  @param lsolver - solver for a linear system
797  @param ls - reference to line search method*/
798 
799  void solve_with_contact(
800  SOLVE_FUNCTION sf,
801  getfem::model& md,
802  gmm::iteration& it_newton,
803  gmm::iteration& it_staggered,
804  const std::string& lsolver,
805  getfem::abstract_newton_line_search &ls);
806 
807 } //end of the namespace level_set_contact
808 
809 #endif //GETFEM_LEVEL_SET_CONTACT_H__
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
getfem_model_solvers.h
Standard solvers for model bricks.
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
level_set_contact::level_set_contact_brick
assembles normal contact terms on the boundary of two contact bodies (master/slave)
Definition: getfem_level_set_contact.h:471
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::NormalTerm
A term, used in level set contact assemblies that builds a surface projection matrix R = N^t X N (whe...
Definition: getfem_level_set_contact.h:511
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::contact_int_method
getfem::pintegration_method contact_int_method() const
integration method on the contact surface, use it when the master is created with a specific integrat...
Definition: getfem_level_set_contact.h:377
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::int_method_descriptor
pintegration_method int_method_descriptor(std::string name, bool throw_if_not_found=true)
Get an integration method from its name .
Definition: getfem_integration.cc:1130
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
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::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
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
getfem::generic_assembly::push_mf
void push_mf(const mesh_fem &mf_)
Add a new mesh_fem.
Definition: getfem_assembling_tensors.h:650
getfem::generic_assembly::push_mat
void push_mat(const MAT &m)
Add a new output matrix (fake const version..)
Definition: getfem_assembling_tensors.h:671
getfem::mesh::sup_region
void sup_region(size_type b)
Remove the region of index b.
Definition: getfem_mesh.h:447
gmm_except.h
Definition of basic exceptions.
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem_deformable_mesh.h
A class adaptor to deform a mesh.
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::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::generic_assembly::push_data
void push_data(const VEC &d)
Add a new data (dense array)
Definition: getfem_assembling_tensors.h:659
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
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
getfem::generic_assembly::assembly
void assembly(const mesh_region &region=mesh_region::all_convexes())
do the assembly on the specified region (boundary or set of convexes)
Definition: getfem_assembling_tensors.cc:1814
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_models.h
Model representation in Getfem.
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
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::model::add_fem_variable
void add_fem_variable(const std::string &name, const mesh_fem &mf, size_type niter=1)
Add a variable being the dofs of a finite element method to the model.
Definition: getfem_models.cc:834
level_set_contact::contact_pair_info::contact_changed
bool contact_changed()
Actual master/slave contact detection.
Definition: getfem_level_set_contact.cc:365
level_set_contact::master_contact_body::integration
const contact_integration integration
integration approach for contact elements that are partially crossed by level sets: PER_ELEMENT - a w...
Definition: getfem_level_set_contact.h:319
getfem::generic_assembly
Generic assembly of vectors, matrices.
Definition: getfem_assembling_tensors.h:607
level_set_contact::master_contact_body::small_weight_multiplier
const scalar_type small_weight_multiplier
in case of REGULARIZED_LEVEL_SET this value scales weight of Gauss points that have negative level se...
Definition: getfem_level_set_contact.h:328
getfem::mesh::bounding_box
void bounding_box(base_node &Pmin, base_node &Pmax) const
Return the bounding box [Pmin - Pmax] of the mesh.
Definition: getfem_mesh.cc:261
level_set_contact::master_contact_body::regularized_tollerance
const scalar_type regularized_tollerance
width of transition for a regularazied Heaviside function in case of REGULARIZED_LEVEL_SET
Definition: getfem_level_set_contact.h:323
level_set_contact::master_contact_body::contact_mim_order
size_type contact_mim_order() const
order of integration of boundary contact terms
Definition: getfem_level_set_contact.h:366
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
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
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
getfem::generic_assembly::push_vec
void push_vec(VEC &v)
Add a new output vector.
Definition: getfem_assembling_tensors.h:663
getfem::generic_assembly::push_mi
void push_mi(const mesh_im &im_)
Add a new mesh_im.
Definition: getfem_assembling_tensors.h:652
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::add_Dirichlet_condition_with_penalization
size_type APIDECL add_Dirichlet_condition_with_penalization(model &md, const mesh_im &mim, const std::string &varname, scalar_type penalization_coeff, size_type region, const std::string &dataname=std::string(), const mesh_fem *mf_mult=0)
Add a Dirichlet condition on the variable varname and the mesh region region.
Definition: getfem_models.cc:4710
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
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::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
getfem::nonlinear_elem_term
abstract class for integration of non-linear terms into the mat_elem computations the nonlinear term ...
Definition: getfem_mat_elem_type.h:67
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
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
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::add_source_term_brick
size_type APIDECL add_source_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region=size_type(-1), const std::string &directdataname=std::string())
Add a source term on the variable varname.
Definition: getfem_models.cc:4127
level_set_contact::solve_with_contact
void solve_with_contact(SOLVE_FUNCTION sf, getfem::model &md, gmm::iteration &it_newton, gmm::iteration &it_staggered, const std::string &lsolver, getfem::abstract_newton_line_search &ls)
Solves a model that has contact in it.
Definition: getfem_level_set_contact.cc:762
level_set_contact::master_contact_body::volume_region
size_type volume_region() const
region of all volume elements without the boundary
Definition: getfem_level_set_contact.h:386
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::generic_assembly::push_nonlinear_term
void push_nonlinear_term(pnonlinear_elem_term net)
Add a new non-linear term.
Definition: getfem_assembling_tensors.h:655
getfem::standard_solve
void standard_solve(model &md, gmm::iteration &iter, rmodel_plsolver_type lsolver, abstract_newton_line_search &ls)
A default solver for the model brick system.
Definition: getfem_model_solvers.cc:410
getfem::add_Laplacian_brick
size_type APIDECL add_Laplacian_brick(model &md, const mesh_im &mim, const std::string &varname, size_type region=size_type(-1))
Add a Laplacian term on the variable varname (in fact with a minus : :math:-\text{div}(\nabla u)).
Definition: getfem_models.cc:3921
level_set_contact::HFunction
Regularized Heaviside function.
Definition: getfem_level_set_contact.h:556
level_set_contact::master_contact_body::max_contact_angle
const scalar_type max_contact_angle
if the angle (in degrees) between contact element and level set contour exceed this value,...
Definition: getfem_level_set_contact.h:333
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::model::set_real_variable
model_real_plain_vector & set_real_variable(const std::string &name, size_type niter) const
Gives the write access to the vector value of a variable.
Definition: getfem_models.cc:3013
level_set_contact::master_contact_body::boundary_region
size_type boundary_region() const
boundary elements, added after creation of the master contact body
Definition: getfem_level_set_contact.h:391
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_assembling.h
Miscelleanous assembly routines for common terms. Use the low-level generic assembly....
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
level_set_contact::master_contact_body::mult_mf_order
const size_type mult_mf_order
approximation order for Lagrange multiplier on the contact surface
Definition: getfem_level_set_contact.h:311