libMesh
linear_elasticity_with_contact.C
Go to the documentation of this file.
1 // Local includes
3 
4 // C++ include files that we need
5 #include <iostream>
6 #include <algorithm>
7 #include <cmath>
8 #include <unordered_set>
9 
10 // Various include files needed for the mesh & solver functionality.
11 #include "libmesh/libmesh.h"
12 #include "libmesh/mesh.h"
13 #include "libmesh/equation_systems.h"
14 #include "libmesh/fe.h"
15 #include "libmesh/quadrature_gauss.h"
16 #include "libmesh/dof_map.h"
17 #include "libmesh/sparse_matrix.h"
18 #include "libmesh/numeric_vector.h"
19 #include "libmesh/dense_matrix.h"
20 #include "libmesh/dense_vector.h"
21 #include "libmesh/elem.h"
22 #include "libmesh/fe_interface.h"
23 #include "libmesh/fe_compute_data.h"
24 #include "libmesh/petsc_matrix.h"
25 #include "libmesh/edge_edge2.h"
26 
27 // The nonlinear solver and system we will be using
28 #include "libmesh/nonlinear_solver.h"
29 #include "libmesh/nonlinear_implicit_system.h"
30 
31 using namespace libMesh;
32 
34  Real contact_penalty_in) :
35  _sys(sys_in),
36  _contact_penalty(contact_penalty_in)
37 {
38 }
39 
41 {
42  _contact_penalty = contact_penalty_in;
43 }
44 
46 {
47  return _contact_penalty;
48 }
49 
51  unsigned int j)
52 {
53  return i == j ? 1. : 0.;
54 }
55 
57  Real poisson_ratio,
58  unsigned int i,
59  unsigned int j,
60  unsigned int k,
61  unsigned int l)
62 {
63  // Define the Lame constants
64  const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
65  const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));
66 
67  return lambda_1 * kronecker_delta(i, j) * kronecker_delta(k, l) +
68  lambda_2 * (kronecker_delta(i, k) * kronecker_delta(j, l) + kronecker_delta(i, l) * kronecker_delta(j, k));
69 }
70 
72  const NumericVector<Number> & input_solution)
73 {
74  // Maintain a set of node ids that we've encountered.
75  std::unordered_set<dof_id_type> encountered_node_ids;
76 
77  // Localize input_solution so that we have the data to move all
78  // elements (not just elements local to this processor).
79  UniquePtr<NumericVector<Number>> localized_input_solution =
80  NumericVector<Number>::build(input_solution.comm());
81 
82  localized_input_solution->init (input_solution.size(), false, SERIAL);
83  input_solution.localize(*localized_input_solution);
84 
85  for (const auto & elem : input_mesh.active_element_ptr_range())
86  {
87  Elem * orig_elem = _sys.get_mesh().elem_ptr(elem->id());
88 
89  for (auto node_id : elem->node_index_range())
90  {
91  Node & node = elem->node_ref(node_id);
92 
93  if (encountered_node_ids.find(node.id()) != encountered_node_ids.end())
94  continue;
95 
96  encountered_node_ids.insert(node.id());
97 
98  std::vector<std::string> uvw_names(3);
99  uvw_names[0] = "u";
100  uvw_names[1] = "v";
101  uvw_names[2] = "w";
102 
103  {
104  const Point master_point = elem->master_point(node_id);
105 
106  Point uvw;
107  for (std::size_t index=0; index<uvw_names.size(); index++)
108  {
109  const unsigned int var = _sys.variable_number(uvw_names[index]);
110  const FEType & fe_type = _sys.get_dof_map().variable_type(var);
111 
112  FEComputeData data (_sys.get_equation_systems(), master_point);
113 
114  FEInterface::compute_data(elem->dim(),
115  fe_type,
116  elem,
117  data);
118 
119  std::vector<dof_id_type> dof_indices_var;
120  _sys.get_dof_map().dof_indices (orig_elem, dof_indices_var, var);
121 
122  for (std::size_t i=0; i<dof_indices_var.size(); i++)
123  {
124  Number value = (*localized_input_solution)(dof_indices_var[i]) * data.shape[i];
125 
126 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
127  // We explicitly store the real part in uvw
128  uvw(index) += value.real();
129 #else
130  uvw(index) += value;
131 #endif
132  }
133  }
134 
135  // Update the node's location
136  node += uvw;
137  }
138  }
139  }
140 }
141 
143 {
144  const MeshBase & mesh = _sys.get_mesh();
145 
146  std::vector<dof_id_type> nodes_on_lower_surface;
147  std::vector<dof_id_type> nodes_on_upper_surface;
148 
149  _lambdas.clear();
150  for (const auto & elem : mesh.active_element_ptr_range())
151  for (auto side : elem->side_index_range())
152  if (elem->neighbor_ptr(side) == libmesh_nullptr)
153  {
154  bool on_lower_contact_surface =
155  mesh.get_boundary_info().has_boundary_id (elem, side, CONTACT_BOUNDARY_LOWER);
156 
157  bool on_upper_contact_surface =
158  mesh.get_boundary_info().has_boundary_id (elem, side, CONTACT_BOUNDARY_UPPER);
159 
160  if (on_lower_contact_surface && on_upper_contact_surface)
161  libmesh_error_msg("Should not be on both surfaces at the same time");
162 
163  if (on_lower_contact_surface || on_upper_contact_surface)
164  {
165  for (auto node_index : elem->node_index_range())
166  if (elem->is_node_on_side(node_index, side))
167  {
168  if (on_lower_contact_surface)
169  nodes_on_lower_surface.push_back(elem->node_id(node_index));
170  else
171  {
172  _lambdas[elem->node_id(node_index)] = 0.;
173  nodes_on_upper_surface.push_back(elem->node_id(node_index));
174  }
175  }
176  }
177  } // end if neighbor(side_) != libmesh_nullptr
178 
179  // In this example, we expect the number of upper and lower nodes to match
180  libmesh_assert(nodes_on_lower_surface.size() == nodes_on_upper_surface.size());
181 
182  // Do an N^2 search to match the contact nodes
183  _contact_node_map.clear();
184  for (std::size_t i=0; i<nodes_on_lower_surface.size(); i++)
185  {
186  dof_id_type lower_node_id = nodes_on_lower_surface[i];
187  Point p_lower = mesh.point(lower_node_id);
188 
189  Real min_distance = std::numeric_limits<Real>::max();
190 
191  for (std::size_t j=0; j<nodes_on_upper_surface.size(); j++)
192  {
193  dof_id_type upper_node_id = nodes_on_upper_surface[j];
194  Point p_upper = mesh.point(upper_node_id);
195 
196  Real distance = (p_upper - p_lower).norm();
197 
198  if (distance < min_distance)
199  {
200  _contact_node_map[lower_node_id] = upper_node_id;
201  min_distance = distance;
202  }
203  }
204  }
205 }
206 
208 {
209  MeshBase & mesh = _sys.get_mesh();
210 
211  std::map<dof_id_type, dof_id_type>::iterator it = _contact_node_map.begin();
212  std::map<dof_id_type, dof_id_type>::iterator it_end = _contact_node_map.end();
213  for( ; it != it_end ; ++it )
214  {
215  dof_id_type master_node_id = it->first;
216  dof_id_type slave_node_id = it->second;
217 
218  Node & master_node = mesh.node_ref(master_node_id);
219  Node & slave_node = mesh.node_ref(slave_node_id);
220 
221  Elem * connector_elem = mesh.add_elem (new Edge2);
222  connector_elem->set_node(0) = &master_node;
223  connector_elem->set_node(1) = &slave_node;
224 
225  connector_elem->subdomain_id() = 10;
226  }
227 
228  mesh.prepare_for_use();
229 }
230 
232  NumericVector<Number> * residual,
233  SparseMatrix<Number> * jacobian,
234  NonlinearImplicitSystem & /*sys*/)
235 {
237  const Real young_modulus = es.parameters.get<Real>("young_modulus");
238  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
239 
240  const MeshBase & mesh = _sys.get_mesh();
241  const unsigned int dim = mesh.mesh_dimension();
242 
243  const unsigned int u_var = _sys.variable_number ("u");
244 
245  DofMap & dof_map = _sys.get_dof_map();
246 
247  FEType fe_type = dof_map.variable_type(u_var);
248  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
249  QGauss qrule (dim, fe_type.default_quadrature_order());
250  fe->attach_quadrature_rule (&qrule);
251 
252  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
253  QGauss qface (dim-1, fe_type.default_quadrature_order());
254  fe_face->attach_quadrature_rule (&qface);
255 
256  UniquePtr<FEBase> fe_neighbor_face (FEBase::build(dim, fe_type));
257  fe_neighbor_face->attach_quadrature_rule (&qface);
258 
259  const std::vector<Real> & JxW = fe->get_JxW();
260  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
261 
262  if (jacobian)
263  jacobian->zero();
264 
265  if (residual)
266  residual->zero();
267 
268  // Do jacobian and residual assembly, including contact forces
270 
271  DenseSubVector<Number> Re_var[3] =
275 
277  DenseSubMatrix<Number> Ke_var[3][3] =
278  {
282  };
283 
284  std::vector<dof_id_type> dof_indices;
285  std::vector<std::vector<dof_id_type>> dof_indices_var(3);
286 
287  for (const auto & elem : mesh.active_local_element_ptr_range())
288  {
289  if( elem->type() == EDGE2 )
290  {
291  // We do not do any assembly on the contact connector elements.
292  // The contact force assembly is handled in a separate loop.
293  continue;
294  }
295 
296  dof_map.dof_indices (elem, dof_indices);
297  for (unsigned int var=0; var<3; var++)
298  dof_map.dof_indices (elem, dof_indices_var[var], var);
299 
300  const unsigned int n_dofs = dof_indices.size();
301  const unsigned int n_var_dofs = dof_indices_var[0].size();
302 
303  fe->reinit (elem);
304 
305  Re.resize (n_dofs);
306  for (unsigned int var=0; var<3; var++)
307  Re_var[var].reposition (var*n_var_dofs, n_var_dofs);
308 
309  Ke.resize (n_dofs, n_dofs);
310  for (unsigned int var_i=0; var_i<3; var_i++)
311  for (unsigned int var_j=0; var_j<3; var_j++)
312  Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
313 
314  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
315  {
316  // Row is variable u, v, or w column is x, y, or z
317  DenseMatrix<Number> grad_u(3, 3);
318  for (unsigned int var_i=0; var_i<3; var_i++)
319  for (unsigned int var_j=0; var_j<3; var_j++)
320  for (unsigned int j=0; j<n_var_dofs; j++)
321  grad_u(var_i, var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
322 
323  // - C_ijkl u_k,l v_i,j
324  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
325  for (unsigned int i=0; i<3; i++)
326  for (unsigned int j=0; j<3; j++)
327  for (unsigned int k=0; k<3; k++)
328  for (unsigned int l=0; l<3; l++)
329  Re_var[i](dof_i) -= JxW[qp] *
330  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * grad_u(k,l) * dphi[dof_i][qp](j);
331 
332  // assemble \int_Omega C_ijkl u_k,l v_i,j \dx
333  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
334  for (unsigned int dof_j=0; dof_j<n_var_dofs; dof_j++)
335  for (unsigned int i=0; i<3; i++)
336  for (unsigned int j=0; j<3; j++)
337  for (unsigned int k=0; k<3; k++)
338  for (unsigned int l=0; l<3; l++)
339  Ke_var[i][k](dof_i, dof_j) -= JxW[qp] *
340  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * dphi[dof_j][qp](l) * dphi[dof_i][qp](j);
341  }
342 
343  dof_map.constrain_element_matrix_and_vector (Ke, Re, dof_indices);
344 
345  if (jacobian)
346  jacobian->add_matrix (Ke, dof_indices);
347 
348  if (residual)
349  residual->add_vector (Re, dof_indices);
350  }
351 
352  // Move a copy of the mesh based on the solution. This is used to get
353  // the contact forces. We could compute the contact forces based on the
354  // current displacement solution, which would not require a clone of the
355  // mesh. Avoiding the mesh clone would be important for production-scale
356  // contact solves, but for the sake of this example, using the clone is
357  // simple and fast enough.
358  UniquePtr<MeshBase> mesh_clone = mesh.clone();
359  move_mesh(*mesh_clone, soln);
360 
361  // Add contributions due to contact penalty forces. Only need to do this on
362  // one processor.
364 
365  std::map<dof_id_type, dof_id_type>::iterator it = _contact_node_map.begin();
366  std::map<dof_id_type, dof_id_type>::iterator it_end = _contact_node_map.end();
367 
368  for ( ; it != it_end; ++it)
369  {
370  dof_id_type lower_point_id = it->first;
371  dof_id_type upper_point_id = it->second;
372 
373  Point upper_to_lower;
374  {
375  Point lower_node_moved = mesh_clone->point(lower_point_id);
376  Point upper_node_moved = mesh_clone->point(upper_point_id);
377 
378  upper_to_lower = (lower_node_moved - upper_node_moved);
379  }
380 
381  // Set the contact force direction. Usually this would be calculated
382  // separately on each master node, but here we just set it to (0, 0, 1)
383  // everywhere.
384  Point contact_force_direction(0., 0., 1.);
385 
386  // gap_function > 0. means that contact has been detected
387  // gap_function < 0. means that we have a gap
388  // This sign convention matches Simo & Laursen (1992).
389  Real gap_function = upper_to_lower * contact_force_direction;
390 
391  // We use the sign of lambda_plus_penalty to determine whether or
392  // not we need to impose a contact force.
393  Real lambda_plus_penalty =
394  (_lambdas[upper_point_id] + gap_function * _contact_penalty);
395 
396  if (lambda_plus_penalty < 0.)
397  lambda_plus_penalty = 0.;
398 
399  // Store lambda_plus_penalty, we'll need to use it later to update _lambdas
400  _lambda_plus_penalty_values[upper_point_id] = lambda_plus_penalty;
401 
402  const Node & lower_node = mesh.node_ref(lower_point_id);
403  const Node & upper_node = mesh.node_ref(upper_point_id);
404 
405  std::vector<dof_id_type> dof_indices_on_lower_node(3);
406  std::vector<dof_id_type> dof_indices_on_upper_node(3);
407  DenseVector<Number> lower_contact_force_vec(3);
408  DenseVector<Number> upper_contact_force_vec(3);
409 
410  for (unsigned int var=0; var<3; var++)
411  {
412  dof_indices_on_lower_node[var] = lower_node.dof_number(_sys.number(), var, 0);
413  lower_contact_force_vec(var) = -lambda_plus_penalty * contact_force_direction(var);
414 
415  dof_indices_on_upper_node[var] = upper_node.dof_number(_sys.number(), var, 0);
416  upper_contact_force_vec(var) = lambda_plus_penalty * contact_force_direction(var);
417  }
418 
419  if (lambda_plus_penalty > 0.)
420  {
421  if (residual && (_sys.comm().rank() == 0))
422  {
423  residual->add_vector (lower_contact_force_vec, dof_indices_on_lower_node);
424  residual->add_vector (upper_contact_force_vec, dof_indices_on_upper_node);
425  }
426 
427  // Add the Jacobian terms due to the contact forces. The lambda contribution
428  // is not relevant here because it doesn't depend on the solution.
429  if (jacobian && (_sys.comm().rank() == 0))
430  for (unsigned int var=0; var<3; var++)
431  for (unsigned int j=0; j<3; j++)
432  {
433  jacobian->add(dof_indices_on_lower_node[var],
434  dof_indices_on_upper_node[j],
435  _contact_penalty * contact_force_direction(j) * contact_force_direction(var));
436 
437  jacobian->add(dof_indices_on_lower_node[var],
438  dof_indices_on_lower_node[j],
439  -_contact_penalty * contact_force_direction(j) * contact_force_direction(var));
440 
441  jacobian->add(dof_indices_on_upper_node[var],
442  dof_indices_on_lower_node[j],
443  _contact_penalty * contact_force_direction(j) * contact_force_direction(var));
444 
445  jacobian->add(dof_indices_on_upper_node[var],
446  dof_indices_on_upper_node[j],
447  -_contact_penalty * contact_force_direction(j) * contact_force_direction(var));
448  }
449  }
450  else
451  {
452  // We add zeros to the matrix even when lambda_plus_penalty = 0.
453  // We do this because some linear algebra libraries (e.g. PETSc)
454  // will condense out any unused entries from the sparsity pattern,
455  // so adding these zeros in ensures that these entries are not
456  // condensed out.
457  if (jacobian && (_sys.comm().rank() == 0))
458  for (unsigned int var=0; var<3; var++)
459  for (unsigned int j=0; j<3; j++)
460  {
461  jacobian->add(dof_indices_on_lower_node[var],
462  dof_indices_on_upper_node[j],
463  0.);
464 
465  jacobian->add(dof_indices_on_lower_node[var],
466  dof_indices_on_lower_node[j],
467  0.);
468 
469  jacobian->add(dof_indices_on_upper_node[var],
470  dof_indices_on_lower_node[j],
471  0.);
472 
473  jacobian->add(dof_indices_on_upper_node[var],
474  dof_indices_on_upper_node[j],
475  0.);
476  }
477  }
478  }
479 }
480 
482 {
484  const Real young_modulus = es.parameters.get<Real>("young_modulus");
485  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
486 
487  const MeshBase & mesh = _sys.get_mesh();
488  const unsigned int dim = mesh.mesh_dimension();
489 
490  unsigned int displacement_vars[3];
491  displacement_vars[0] = _sys.variable_number ("u");
492  displacement_vars[1] = _sys.variable_number ("v");
493  displacement_vars[2] = _sys.variable_number ("w");
494  const unsigned int u_var = _sys.variable_number ("u");
495 
496  const DofMap & dof_map = _sys.get_dof_map();
497  FEType fe_type = dof_map.variable_type(u_var);
498  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
499  QGauss qrule (dim, fe_type.default_quadrature_order());
500  fe->attach_quadrature_rule (&qrule);
501 
502  const std::vector<Real> & JxW = fe->get_JxW();
503  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
504 
505  // Also, get a reference to the ExplicitSystem
506  ExplicitSystem & stress_system = es.get_system<ExplicitSystem>("StressSystem");
507  const DofMap & stress_dof_map = stress_system.get_dof_map();
508  unsigned int sigma_vars[6];
509  sigma_vars[0] = stress_system.variable_number ("sigma_00");
510  sigma_vars[1] = stress_system.variable_number ("sigma_01");
511  sigma_vars[2] = stress_system.variable_number ("sigma_02");
512  sigma_vars[3] = stress_system.variable_number ("sigma_11");
513  sigma_vars[4] = stress_system.variable_number ("sigma_12");
514  sigma_vars[5] = stress_system.variable_number ("sigma_22");
515  unsigned int vonMises_var = stress_system.variable_number ("vonMises");
516 
517  // Storage for the stress dof indices on each element
518  std::vector<std::vector<dof_id_type>> dof_indices_var(_sys.n_vars());
519  std::vector<dof_id_type> stress_dof_indices_var;
520 
521  // To store the stress tensor on each element
522  DenseMatrix<Number> elem_avg_stress_tensor(3, 3);
523 
524  for (const auto & elem : mesh.active_local_element_ptr_range())
525  {
526  if( elem->type() == EDGE2 )
527  {
528  // We do not compute stress on the contact connector elements.
529  continue;
530  }
531 
532  for (unsigned int var=0; var<3; var++)
533  dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
534 
535  const unsigned int n_var_dofs = dof_indices_var[0].size();
536 
537  fe->reinit (elem);
538 
539  // clear the stress tensor
540  elem_avg_stress_tensor.resize(3, 3);
541 
542  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
543  {
544  // Row is variable u1, u2, or u3, column is x, y, or z
545  DenseMatrix<Number> grad_u(3, 3);
546  for (unsigned int var_i=0; var_i<3; var_i++)
547  for (unsigned int var_j=0; var_j<3; var_j++)
548  for (unsigned int j=0; j<n_var_dofs; j++)
549  grad_u(var_i, var_j) +=
550  dphi[j][qp](var_j) * _sys.current_solution(dof_indices_var[var_i][j]);
551 
552  DenseMatrix<Number> stress_tensor(3, 3);
553  for (unsigned int i=0; i<3; i++)
554  for (unsigned int j=0; j<3; j++)
555  for (unsigned int k=0; k<3; k++)
556  for (unsigned int l=0; l<3; l++)
557  stress_tensor(i,j) +=
558  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * grad_u(k,l);
559 
560  // We want to plot the average stress on each element, hence
561  // we integrate stress_tensor
562  elem_avg_stress_tensor.add(JxW[qp], stress_tensor);
563  }
564 
565  // Get the average stress per element by dividing by volume
566  elem_avg_stress_tensor.scale(1./elem->volume());
567 
568  // load elem_sigma data into stress_system
569  unsigned int stress_var_index = 0;
570  for (unsigned int i=0; i<3; i++)
571  for (unsigned int j=i; j<3; j++)
572  {
573  stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[stress_var_index]);
574 
575  // We are using CONSTANT MONOMIAL basis functions, hence we only need to get
576  // one dof index per variable
577  dof_id_type dof_index = stress_dof_indices_var[0];
578 
579  if ((stress_system.solution->first_local_index() <= dof_index) &&
580  (dof_index < stress_system.solution->last_local_index()))
581  stress_system.solution->set(dof_index, elem_avg_stress_tensor(i,j));
582 
583  stress_var_index++;
584  }
585 
586  // Also, the von Mises stress
587  Number vonMises_value = std::sqrt(0.5*(pow(elem_avg_stress_tensor(0,0) - elem_avg_stress_tensor(1,1), 2.) +
588  pow(elem_avg_stress_tensor(1,1) - elem_avg_stress_tensor(2,2), 2.) +
589  pow(elem_avg_stress_tensor(2,2) - elem_avg_stress_tensor(0,0), 2.) +
590  6.*(pow(elem_avg_stress_tensor(0,1), 2.) +
591  pow(elem_avg_stress_tensor(1,2), 2.) +
592  pow(elem_avg_stress_tensor(2,0), 2.))));
593 
594  stress_dof_map.dof_indices (elem, stress_dof_indices_var, vonMises_var);
595  dof_id_type dof_index = stress_dof_indices_var[0];
596 
597  if ((stress_system.solution->first_local_index() <= dof_index) &&
598  (dof_index < stress_system.solution->last_local_index()))
599  stress_system.solution->set(dof_index, vonMises_value);
600  }
601 
602  // Should call close and update when we set vector entries directly
603  stress_system.solution->close();
604  stress_system.update();
605 }
606 
608 {
609  Real max_delta_lambda = 0.;
610  Real max_new_lambda = 0.;
611 
612  std::map<dof_id_type, Real>::iterator it = _lambdas.begin();
613  std::map<dof_id_type, Real>::iterator it_end = _lambdas.end();
614  for ( ; it != it_end; ++it)
615  {
616  dof_id_type upper_node_id = it->first;
617 
618  std::map<dof_id_type, Real>::iterator new_lambda_it = _lambda_plus_penalty_values.find(upper_node_id);
619  if (new_lambda_it == _lambda_plus_penalty_values.end())
620  libmesh_error_msg("New lambda value not found");
621 
622  Real new_lambda = new_lambda_it->second;
623  Real old_lambda = it->second;
624 
625  it->second = new_lambda;
626 
627  Real delta_lambda = std::abs(new_lambda-old_lambda);
628  if (delta_lambda > max_delta_lambda)
629  max_delta_lambda = delta_lambda;
630 
631  if (std::abs(new_lambda) > max_new_lambda)
632  max_new_lambda = std::abs(new_lambda);
633  }
634 
635  return std::make_pair(max_delta_lambda, max_new_lambda);
636 }
637 
639 {
640  UniquePtr<MeshBase> mesh_clone = _sys.get_mesh().clone();
641  move_mesh(*mesh_clone, *_sys.solution);
642 
643  Real least_value = std::numeric_limits<Real>::max();
644  Real max_value = std::numeric_limits<Real>::min();
645 
646  std::map<dof_id_type, dof_id_type>::iterator it = _contact_node_map.begin();
647  std::map<dof_id_type, dof_id_type>::iterator it_end = _contact_node_map.end();
648  for ( ; it != it_end; ++it)
649  {
650  dof_id_type lower_point_id = it->first;
651  dof_id_type upper_point_id = it->second;
652 
653  Point upper_to_lower;
654  {
655  Point lower_node_moved = mesh_clone->point(lower_point_id);
656  Point upper_node_moved = mesh_clone->point(upper_point_id);
657 
658  upper_to_lower = (lower_node_moved - upper_node_moved);
659  }
660 
661  // Set the contact force direction. Usually this would be calculated
662  // separately on each master node, but here we just set it to (0, 0, 1)
663  // everywhere.
664  Point contact_force_direction(0., 0., 1.);
665 
666  // gap_function > 0. means that contact has been detected
667  // gap_function < 0. means that we have a gap
668  // This sign convention matches Simo & Laursen (1992).
669  Real gap_function = upper_to_lower * contact_force_direction;
670 
671  if (gap_function < least_value)
672  least_value = gap_function;
673 
674  if (gap_function > max_value)
675  max_value = gap_function;
676  }
677 
678  return std::make_pair(least_value, max_value);
679 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
double abs(double a)
This is the EquationSystems class.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
A Node is like a Point, but with more information.
Definition: node.h:52
virtual numeric_index_type size() const =0
virtual const Point & point(const dof_id_type i) const =0
Real get_contact_penalty() const
Get the penalty parameter.
std::map< dof_id_type, Real > _lambda_plus_penalty_values
Store the intermediate values of lambda plus penalty.
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int dim
class FEComputeData hides arbitrary data to be passed to and from children of FEBase through the FEIn...
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:420
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
Real _contact_penalty
Penalize overlapping elements.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
const T & get(const std::string &) const
Definition: parameters.h:430
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
The libMesh namespace provides an interface to certain functionality in the library.
void initialize_contact_load_paths()
Set up the load paths on the contact surfaces.
long double max(long double a, double b)
Real distance(const Point &p)
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void zero()=0
Set all entries to zero.
void move_mesh(MeshBase &input_mesh, const NumericVector< Number > &input_solution)
Move the mesh nodes of input_mesh based on the displacement field in input_solution.
libmesh_assert(j)
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
virtual void zero()=0
Set all entries to 0.
void add_contact_edge_elements()
Add edge elements into the mesh based on the contact load paths.
std::map< dof_id_type, Real > _lambdas
Augmented Lagrangian values at each contact node.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
Order default_quadrature_order() const
Definition: fe_type.h:332
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseMatrix< T3 > &mat)
Adds factor times mat to this matrix.
Definition: dense_matrix.h:883
Defines a dense submatrix for use in Finite Element-type computations.
virtual void residual_and_jacobian(const NumericVector< Number > &soln, NumericVector< Number > *residual, SparseMatrix< Number > *jacobian, NonlinearImplicitSystem &)
Evaluate the Jacobian of the nonlinear system.
This class provides a specific system class.
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:1806
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:174
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
std::pair< Real, Real > update_lambdas()
Update the lambda parameters in the augmented Lagrangian method.
double pow(double a, int b)
Real elasticity_tensor(Real young_modulus, Real poisson_ratio, unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
void set_contact_penalty(Real contact_penalty_in)
Update the penalty parameter.
void scale(const T factor)
Multiplies every element in the matrix by factor.
Definition: dense_matrix.h:851
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
unsigned int number() const
Definition: system.h:2006
std::map< dof_id_type, dof_id_type > _contact_node_map
This provides a map between contact nodes.
const Parallel::Communicator & comm() const
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
static const bool value
Definition: xdr_io.C:108
void compute_stresses()
Compute the Cauchy stress for the current solution.
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
The Edge2 is an element in 1D composed of 2 nodes.
Definition: edge_edge2.h:43
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
LinearElasticityWithContact(NonlinearImplicitSystem &sys_in, Real contact_penalty_in)
Constructor.
unsigned int rank() const
Definition: parallel.h:724
Real kronecker_delta(unsigned int i, unsigned int j)
Kronecker delta function.
const T_sys & get_system(const std::string &name) const
unsigned int n_vars() const
Definition: system.h:2086
IterBase * data
Ideally this private member data should have protected access.
dof_id_type id() const
Definition: dof_object.h:632
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual UniquePtr< MeshBase > clone() const =0
Virtual "copy constructor".
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
std::pair< Real, Real > get_least_and_max_gap_function()
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
The ExplicitSystem provides only "right hand side" storage, which should be sufficient for solving mo...
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917