libMesh
Public Member Functions | Private Attributes | List of all members
LinearElasticityWithContact Class Reference

This class encapsulate all functionality required for assembling and solving a linear elastic model with contact. More...

#include <linear_elasticity_with_contact.h>

Inheritance diagram for LinearElasticityWithContact:
[legend]

Public Member Functions

 LinearElasticityWithContact (NonlinearImplicitSystem &sys_in, Real contact_penalty_in)
 Constructor. More...
 
void set_contact_penalty (Real contact_penalty_in)
 Update the penalty parameter. More...
 
Real get_contact_penalty () const
 Get the penalty parameter. More...
 
Real kronecker_delta (unsigned int i, unsigned int j)
 Kronecker delta function. More...
 
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. More...
 
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. More...
 
void initialize_contact_load_paths ()
 Set up the load paths on the contact surfaces. More...
 
void add_contact_edge_elements ()
 Add edge elements into the mesh based on the contact load paths. More...
 
virtual void residual_and_jacobian (const NumericVector< Number > &soln, NumericVector< Number > *residual, SparseMatrix< Number > *jacobian, NonlinearImplicitSystem &)
 Evaluate the Jacobian of the nonlinear system. More...
 
void compute_stresses ()
 Compute the Cauchy stress for the current solution. More...
 
std::pair< Real, Real > update_lambdas ()
 Update the lambda parameters in the augmented Lagrangian method. More...
 
std::pair< Real, Real > get_least_and_max_gap_function ()
 

Private Attributes

NonlinearImplicitSystem_sys
 Keep a reference to the NonlinearImplicitSystem. More...
 
Real _contact_penalty
 Penalize overlapping elements. More...
 
std::map< dof_id_type, Real > _lambda_plus_penalty_values
 Store the intermediate values of lambda plus penalty. More...
 
std::map< dof_id_type, Real > _lambdas
 Augmented Lagrangian values at each contact node. More...
 
std::map< dof_id_type, dof_id_type > _contact_node_map
 This provides a map between contact nodes. More...
 

Detailed Description

This class encapsulate all functionality required for assembling and solving a linear elastic model with contact.

Definition at line 31 of file linear_elasticity_with_contact.h.

Constructor & Destructor Documentation

◆ LinearElasticityWithContact()

LinearElasticityWithContact::LinearElasticityWithContact ( NonlinearImplicitSystem sys_in,
Real  contact_penalty_in 
)

Constructor.

Definition at line 35 of file linear_elasticity_with_contact.C.

36  :
37  _sys(sys_in),
38  _contact_penalty(contact_penalty_in)
39 {
40 }
Real _contact_penalty
Penalize overlapping elements.
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.

Member Function Documentation

◆ add_contact_edge_elements()

void LinearElasticityWithContact::add_contact_edge_elements ( )

Add edge elements into the mesh based on the contact load paths.

This ensure proper parallel communication of data, e.g. if a node on one side of the contact surface has a constraint on it, then adding contact elements into the mesh ensures that the constraint will be enforced properly.

Definition at line 209 of file linear_elasticity_with_contact.C.

References _contact_node_map, _sys, libMesh::MeshBase::add_elem(), libMesh::EDGE2, libMesh::System::get_mesh(), mesh, libMesh::MeshBase::node_ref(), libMesh::MeshBase::prepare_for_use(), libMesh::Elem::set_node(), and libMesh::Elem::subdomain_id().

210 {
211  MeshBase & mesh = _sys.get_mesh();
212 
213  for (const auto & [master_node_id, slave_node_id] : _contact_node_map)
214  {
215  Node & master_node = mesh.node_ref(master_node_id);
216  Node & slave_node = mesh.node_ref(slave_node_id);
217 
218  Elem * connector_elem = mesh.add_elem(Elem::build(EDGE2));
219  connector_elem->set_node(0) = &master_node;
220  connector_elem->set_node(1) = &slave_node;
221 
222  connector_elem->subdomain_id() = 10;
223  }
224 
226 }
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2381
A Node is like a Point, but with more information.
Definition: node.h:52
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:710
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
const MeshBase & get_mesh() const
Definition: system.h:2277
This is the MeshBase class.
Definition: mesh_base.h:74
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
subdomain_id_type subdomain_id() const
Definition: elem.h:2391
std::map< dof_id_type, dof_id_type > _contact_node_map
This provides a map between contact nodes.
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:575

◆ compute_stresses()

void LinearElasticityWithContact::compute_stresses ( )

Compute the Cauchy stress for the current solution.

Definition at line 472 of file linear_elasticity_with_contact.C.

References _sys, libMesh::TypeTensor< T >::add_scaled(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), libMesh::EDGE2, elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::System::n_vars(), libMesh::EquationSystems::parameters, libMesh::Utility::pow(), libMesh::Real, libMesh::System::solution, std::sqrt(), libMesh::System::update(), libMesh::System::variable_number(), libMesh::DofMap::variable_type(), and libMesh::TypeTensor< T >::zero().

473 {
475  const Real young_modulus = es.parameters.get<Real>("young_modulus");
476  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
477 
478  const MeshBase & mesh = _sys.get_mesh();
479  const unsigned int dim = mesh.mesh_dimension();
480 
481  unsigned int displacement_vars[] = {
482  _sys.variable_number ("u"),
483  _sys.variable_number ("v"),
484  _sys.variable_number ("w")};
485  const unsigned int u_var = _sys.variable_number ("u");
486 
487  const DofMap & dof_map = _sys.get_dof_map();
488  FEType fe_type = dof_map.variable_type(u_var);
489  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
490  QGauss qrule (dim, fe_type.default_quadrature_order());
491  fe->attach_quadrature_rule (&qrule);
492 
493  const std::vector<Real> & JxW = fe->get_JxW();
494  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
495 
496  // Also, get a reference to the ExplicitSystem
497  ExplicitSystem & stress_system = es.get_system<ExplicitSystem>("StressSystem");
498  const DofMap & stress_dof_map = stress_system.get_dof_map();
499  unsigned int sigma_vars[] = {
500  stress_system.variable_number ("sigma_00"),
501  stress_system.variable_number ("sigma_01"),
502  stress_system.variable_number ("sigma_02"),
503  stress_system.variable_number ("sigma_11"),
504  stress_system.variable_number ("sigma_12"),
505  stress_system.variable_number ("sigma_22")};
506  unsigned int vonMises_var = stress_system.variable_number ("vonMises");
507 
508  // Storage for the stress dof indices on each element
509  std::vector<std::vector<dof_id_type>> dof_indices_var(_sys.n_vars());
510  std::vector<dof_id_type> stress_dof_indices_var;
511 
512  // To store the stress tensor on each element
513  TensorValue<Number> elem_avg_stress_tensor;
514 
515  for (const auto & elem : mesh.active_local_element_ptr_range())
516  {
517  if( elem->type() == EDGE2 )
518  {
519  // We do not compute stress on the contact connector elements.
520  continue;
521  }
522 
523  for (unsigned int var=0; var<3; var++)
524  dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
525 
526  const unsigned int n_var_dofs = dof_indices_var[0].size();
527 
528  fe->reinit (elem);
529 
530  // clear the stress tensor
531  elem_avg_stress_tensor.zero();
532 
533  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
534  {
535  // Row is variable u1, u2, or u3, column is x, y, or z
536  TensorValue<Number> grad_u;
537  for (unsigned int var_i=0; var_i<3; var_i++)
538  for (unsigned int var_j=0; var_j<3; var_j++)
539  for (unsigned int j=0; j<n_var_dofs; j++)
540  grad_u(var_i, var_j) +=
541  dphi[j][qp](var_j) * _sys.current_solution(dof_indices_var[var_i][j]);
542 
543  TensorValue<Number> stress_tensor;
544  for (unsigned int i=0; i<3; i++)
545  for (unsigned int j=0; j<3; j++)
546  for (unsigned int k=0; k<3; k++)
547  for (unsigned int l=0; l<3; l++)
548  stress_tensor(i,j) +=
549  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * grad_u(k,l);
550 
551  // We want to plot the average stress on each element, hence
552  // we integrate stress_tensor
553  elem_avg_stress_tensor.add_scaled(stress_tensor, JxW[qp]);
554  }
555 
556  // Get the average stress per element by dividing by volume
557  elem_avg_stress_tensor /= elem->volume();
558 
559  // load elem_sigma data into stress_system
560  unsigned int stress_var_index = 0;
561  for (unsigned int i=0; i<3; i++)
562  for (unsigned int j=i; j<3; j++)
563  {
564  stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[stress_var_index]);
565 
566  // We are using CONSTANT MONOMIAL basis functions, hence we only need to get
567  // one dof index per variable
568  dof_id_type dof_index = stress_dof_indices_var[0];
569 
570  if ((stress_system.solution->first_local_index() <= dof_index) &&
571  (dof_index < stress_system.solution->last_local_index()))
572  stress_system.solution->set(dof_index, elem_avg_stress_tensor(i,j));
573 
574  stress_var_index++;
575  }
576 
577  // Also, the von Mises stress
578  Number vonMises_value = std::sqrt(0.5*(pow(elem_avg_stress_tensor(0,0) - elem_avg_stress_tensor(1,1), 2.) +
579  pow(elem_avg_stress_tensor(1,1) - elem_avg_stress_tensor(2,2), 2.) +
580  pow(elem_avg_stress_tensor(2,2) - elem_avg_stress_tensor(0,0), 2.) +
581  6.*(pow(elem_avg_stress_tensor(0,1), 2.) +
582  pow(elem_avg_stress_tensor(1,2), 2.) +
583  pow(elem_avg_stress_tensor(2,0), 2.))));
584 
585  stress_dof_map.dof_indices (elem, stress_dof_indices_var, vonMises_var);
586  dof_id_type dof_index = stress_dof_indices_var[0];
587 
588  if ((stress_system.solution->first_local_index() <= dof_index) &&
589  (dof_index < stress_system.solution->last_local_index()))
590  stress_system.solution->set(dof_index, vonMises_value);
591  }
592 
593  // Should call close and update when we set vector entries directly
594  stress_system.solution->close();
595  stress_system.update();
596 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
This is the EquationSystems class.
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:1992
unsigned int dim
void zero()
Set all entries of the tensor to 0.
Definition: type_tensor.h:1338
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2144
const EquationSystems & get_equation_systems() const
Definition: system.h:730
MeshBase & mesh
Order default_quadrature_order() const
Definition: fe_type.h:357
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:53
const T_sys & get_system(std::string_view name) const
const MeshBase & get_mesh() const
Definition: system.h:2277
This is the MeshBase class.
Definition: mesh_base.h:74
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:157
unsigned int variable_number(std::string_view var) const
Definition: system.C:1557
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:169
T pow(const T &x)
Definition: utility.h:328
const T & get(std::string_view) const
Definition: parameters.h:426
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1573
void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
Definition: type_tensor.h:851
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.
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:493
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:324
Parameters parameters
Data structure holding arbitrary parameters.
unsigned int n_vars() const
Definition: system.h:2349
const DofMap & get_dof_map() const
Definition: system.h:2293
Manages consistently variables, degrees of freedom, and coefficient vectors for explicit systems...
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
uint8_t dof_id_type
Definition: id_types.h:67

◆ elasticity_tensor()

Real LinearElasticityWithContact::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.

Definition at line 58 of file linear_elasticity_with_contact.C.

References kronecker_delta(), and libMesh::Real.

Referenced by compute_stresses(), and residual_and_jacobian().

64 {
65  // Define the Lame constants
66  const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
67  const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));
68 
69  return lambda_1 * kronecker_delta(i, j) * kronecker_delta(k, l) +
70  lambda_2 * (kronecker_delta(i, k) * kronecker_delta(j, l) + kronecker_delta(i, l) * kronecker_delta(j, k));
71 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real kronecker_delta(unsigned int i, unsigned int j)
Kronecker delta function.

◆ get_contact_penalty()

Real LinearElasticityWithContact::get_contact_penalty ( ) const

Get the penalty parameter.

Definition at line 47 of file linear_elasticity_with_contact.C.

References _contact_penalty.

48 {
49  return _contact_penalty;
50 }
Real _contact_penalty
Penalize overlapping elements.

◆ get_least_and_max_gap_function()

std::pair< Real, Real > LinearElasticityWithContact::get_least_and_max_gap_function ( )
Returns
the least and max gap function values for the current solution.

Definition at line 624 of file linear_elasticity_with_contact.C.

References _contact_node_map, _sys, libMesh::MeshBase::clone(), libMesh::System::get_mesh(), move_mesh(), libMesh::Real, and libMesh::System::solution.

625 {
626  std::unique_ptr<MeshBase> mesh_clone = _sys.get_mesh().clone();
627  move_mesh(*mesh_clone, *_sys.solution);
628 
629  Real least_value = std::numeric_limits<Real>::max();
630  Real max_value = std::numeric_limits<Real>::min();
631 
632  for (const auto & [lower_point_id, upper_point_id] : _contact_node_map)
633  {
634  Point upper_to_lower;
635  {
636  Point lower_node_moved = mesh_clone->point(lower_point_id);
637  Point upper_node_moved = mesh_clone->point(upper_point_id);
638 
639  upper_to_lower = (lower_node_moved - upper_node_moved);
640  }
641 
642  // Set the contact force direction. Usually this would be calculated
643  // separately on each master node, but here we just set it to (0, 0, 1)
644  // everywhere.
645  Point contact_force_direction(0., 0., 1.);
646 
647  // gap_function > 0. means that contact has been detected
648  // gap_function < 0. means that we have a gap
649  // This sign convention matches Simo & Laursen (1992).
650  Real gap_function = upper_to_lower * contact_force_direction;
651 
652  if (gap_function < least_value)
653  least_value = gap_function;
654 
655  if (gap_function > max_value)
656  max_value = gap_function;
657  }
658 
659  return std::make_pair(least_value, max_value);
660 }
const MeshBase & get_mesh() const
Definition: system.h:2277
virtual std::unique_ptr< MeshBase > clone() const =0
Virtual "copy constructor".
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.
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1573
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< dof_id_type, dof_id_type > _contact_node_map
This provides a map between contact nodes.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ initialize_contact_load_paths()

void LinearElasticityWithContact::initialize_contact_load_paths ( )

Set up the load paths on the contact surfaces.

Definition at line 144 of file linear_elasticity_with_contact.C.

References _contact_node_map, _lambdas, _sys, distance(), libMesh::MeshBase::get_boundary_info(), libMesh::System::get_mesh(), libMesh::BoundaryInfo::has_boundary_id(), libMesh::libmesh_assert(), mesh, libMesh::TensorTools::norm(), libMesh::MeshBase::point(), and libMesh::Real.

145 {
146  const MeshBase & mesh = _sys.get_mesh();
147 
148  std::vector<dof_id_type> nodes_on_lower_surface;
149  std::vector<dof_id_type> nodes_on_upper_surface;
150 
151  _lambdas.clear();
152  for (const auto & elem : mesh.active_element_ptr_range())
153  for (auto side : elem->side_index_range())
154  if (elem->neighbor_ptr(side) == nullptr)
155  {
156  bool on_lower_contact_surface =
157  mesh.get_boundary_info().has_boundary_id (elem, side, CONTACT_BOUNDARY_LOWER);
158 
159  bool on_upper_contact_surface =
160  mesh.get_boundary_info().has_boundary_id (elem, side, CONTACT_BOUNDARY_UPPER);
161 
162  libmesh_error_msg_if(on_lower_contact_surface && on_upper_contact_surface,
163  "Should not be on both surfaces at the same time");
164 
165  if (on_lower_contact_surface || on_upper_contact_surface)
166  {
167  for (auto node_index : elem->node_index_range())
168  if (elem->is_node_on_side(node_index, side))
169  {
170  if (on_lower_contact_surface)
171  nodes_on_lower_surface.push_back(elem->node_id(node_index));
172  else
173  {
174  _lambdas[elem->node_id(node_index)] = 0.;
175  nodes_on_upper_surface.push_back(elem->node_id(node_index));
176  }
177  }
178  }
179  } // end if neighbor(side_) != nullptr
180 
181  // In this example, we expect the number of upper and lower nodes to match
182  libmesh_assert(nodes_on_lower_surface.size() == nodes_on_upper_surface.size());
183 
184  // Do an N^2 search to match the contact nodes
185  _contact_node_map.clear();
186  for (std::size_t i=0; i<nodes_on_lower_surface.size(); i++)
187  {
188  dof_id_type lower_node_id = nodes_on_lower_surface[i];
189  Point p_lower = mesh.point(lower_node_id);
190 
191  Real min_distance = std::numeric_limits<Real>::max();
192 
193  for (std::size_t j=0; j<nodes_on_upper_surface.size(); j++)
194  {
195  dof_id_type upper_node_id = nodes_on_upper_surface[j];
196  Point p_upper = mesh.point(upper_node_id);
197 
198  Real distance = (p_upper - p_lower).norm();
199 
200  if (distance < min_distance)
201  {
202  _contact_node_map[lower_node_id] = upper_node_id;
203  min_distance = distance;
204  }
205  }
206  }
207 }
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
MeshBase & mesh
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:159
const MeshBase & get_mesh() const
Definition: system.h:2277
Real distance(const Point &p)
This is the MeshBase class.
Definition: mesh_base.h:74
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
std::map< dof_id_type, Real > _lambdas
Augmented Lagrangian values at each contact node.
libmesh_assert(ctx)
auto norm(const T &a) -> decltype(std::abs(a))
Definition: tensor_tools.h:74
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< dof_id_type, dof_id_type > _contact_node_map
This provides a map between contact nodes.
virtual const Point & point(const dof_id_type i) const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
uint8_t dof_id_type
Definition: id_types.h:67

◆ kronecker_delta()

Real LinearElasticityWithContact::kronecker_delta ( unsigned int  i,
unsigned int  j 
)

Kronecker delta function.

Definition at line 52 of file linear_elasticity_with_contact.C.

Referenced by elasticity_tensor().

54 {
55  return i == j ? 1. : 0.;
56 }

◆ move_mesh()

void LinearElasticityWithContact::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.

Definition at line 73 of file linear_elasticity_with_contact.C.

References _sys, libMesh::ParallelObject::comm(), libMesh::DofMap::dof_indices(), libMesh::MeshBase::elem_ptr(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::DofObject::id(), libMesh::NumericVector< T >::init(), libMesh::NumericVector< T >::localize(), libMesh::SERIAL, libMesh::NumericVector< T >::size(), value, libMesh::System::variable_number(), and libMesh::DofMap::variable_type().

Referenced by get_least_and_max_gap_function(), and residual_and_jacobian().

75 {
76  // Maintain a set of node ids that we've encountered.
77  std::unordered_set<dof_id_type> encountered_node_ids;
78 
79  // Localize input_solution so that we have the data to move all
80  // elements (not just elements local to this processor).
81  std::unique_ptr<NumericVector<Number>> localized_input_solution =
82  NumericVector<Number>::build(input_solution.comm());
83 
84  localized_input_solution->init (input_solution.size(), false, SERIAL);
85  input_solution.localize(*localized_input_solution);
86 
87  for (const auto & elem : input_mesh.active_element_ptr_range())
88  {
89  Elem * orig_elem = _sys.get_mesh().elem_ptr(elem->id());
90 
91  for (auto node_id : elem->node_index_range())
92  {
93  Node & node = elem->node_ref(node_id);
94 
95  if (encountered_node_ids.find(node.id()) != encountered_node_ids.end())
96  continue;
97 
98  encountered_node_ids.insert(node.id());
99 
100  std::vector<std::string> uvw_names(3);
101  uvw_names[0] = "u";
102  uvw_names[1] = "v";
103  uvw_names[2] = "w";
104 
105  {
106  const Point master_point = elem->master_point(node_id);
107 
108  Point uvw;
109  for (std::size_t index=0; index<uvw_names.size(); index++)
110  {
111  const unsigned int var = _sys.variable_number(uvw_names[index]);
112  const FEType & fe_type = _sys.get_dof_map().variable_type(var);
113 
114  FEComputeData data (_sys.get_equation_systems(), master_point);
115 
116  FEInterface::compute_data(elem->dim(),
117  fe_type,
118  elem,
119  data);
120 
121  std::vector<dof_id_type> dof_indices_var;
122  _sys.get_dof_map().dof_indices (orig_elem, dof_indices_var, var);
123 
124  for (std::size_t i=0; i<dof_indices_var.size(); i++)
125  {
126  Number value = (*localized_input_solution)(dof_indices_var[i]) * data.shape[i];
127 
128 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
129  // We explicitly store the real part in uvw
130  uvw(index) += value.real();
131 #else
132  uvw(index) += value;
133 #endif
134  }
135  }
136 
137  // Update the node's location
138  node += uvw;
139  }
140  }
141  }
142 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
A Node is like a Point, but with more information.
Definition: node.h:52
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:1992
class FEComputeData hides arbitrary data to be passed to and from children of FEBase through the FEIn...
virtual numeric_index_type size() const =0
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2144
const EquationSystems & get_equation_systems() const
Definition: system.h:730
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: vector_fe_ex5.C:43
const Parallel::Communicator & comm() const
const MeshBase & get_mesh() const
Definition: system.h:2277
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
Change the dimension of the vector to n.
unsigned int variable_number(std::string_view var) const
Definition: system.C:1557
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
dof_id_type id() const
Definition: dof_object.h:823
virtual const Elem * elem_ptr(const dof_id_type i) const =0
static const bool value
Definition: xdr_io.C:54
const DofMap & get_dof_map() const
Definition: system.h:2293
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.

◆ residual_and_jacobian()

void LinearElasticityWithContact::residual_and_jacobian ( const NumericVector< Number > &  soln,
NumericVector< Number > *  residual,
SparseMatrix< Number > *  jacobian,
NonlinearImplicitSystem  
)
virtual

Evaluate the Jacobian of the nonlinear system.

Definition at line 228 of file linear_elasticity_with_contact.C.

References _contact_node_map, _contact_penalty, _lambda_plus_penalty_values, _lambdas, _sys, libMesh::SparseMatrix< T >::add(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::MeshBase::clone(), libMesh::ParallelObject::comm(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), libMesh::DofObject::dof_number(), libMesh::EDGE2, elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), mesh, libMesh::MeshBase::mesh_dimension(), move_mesh(), libMesh::MeshBase::node_ref(), libMesh::System::number(), libMesh::EquationSystems::parameters, TIMPI::Communicator::rank(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::System::variable_number(), libMesh::DofMap::variable_type(), libMesh::NumericVector< T >::zero(), and libMesh::SparseMatrix< T >::zero().

232 {
234  const Real young_modulus = es.parameters.get<Real>("young_modulus");
235  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
236 
237  const MeshBase & mesh = _sys.get_mesh();
238  const unsigned int dim = mesh.mesh_dimension();
239 
240  const unsigned int u_var = _sys.variable_number ("u");
241 
242  DofMap & dof_map = _sys.get_dof_map();
243 
244  FEType fe_type = dof_map.variable_type(u_var);
245  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
246  QGauss qrule (dim, fe_type.default_quadrature_order());
247  fe->attach_quadrature_rule (&qrule);
248 
249  std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
250  QGauss qface (dim-1, fe_type.default_quadrature_order());
251  fe_face->attach_quadrature_rule (&qface);
252 
253  std::unique_ptr<FEBase> fe_neighbor_face (FEBase::build(dim, fe_type));
254  fe_neighbor_face->attach_quadrature_rule (&qface);
255 
256  const std::vector<Real> & JxW = fe->get_JxW();
257  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
258 
259  if (jacobian)
260  jacobian->zero();
261 
262  if (residual)
263  residual->zero();
264 
265  // Do jacobian and residual assembly, including contact forces
267 
268  DenseSubVector<Number> Re_var[3] =
272 
274  DenseSubMatrix<Number> Ke_var[3][3] =
275  {
279  };
280 
281  std::vector<dof_id_type> dof_indices;
282  std::vector<std::vector<dof_id_type>> dof_indices_var(3);
283 
284  for (const auto & elem : mesh.active_local_element_ptr_range())
285  {
286  if( elem->type() == EDGE2 )
287  {
288  // We do not do any assembly on the contact connector elements.
289  // The contact force assembly is handled in a separate loop.
290  continue;
291  }
292 
293  dof_map.dof_indices (elem, dof_indices);
294  for (unsigned int var=0; var<3; var++)
295  dof_map.dof_indices (elem, dof_indices_var[var], var);
296 
297  const unsigned int n_dofs = dof_indices.size();
298  const unsigned int n_var_dofs = dof_indices_var[0].size();
299 
300  fe->reinit (elem);
301 
302  Re.resize (n_dofs);
303  for (unsigned int var=0; var<3; var++)
304  Re_var[var].reposition (var*n_var_dofs, n_var_dofs);
305 
306  Ke.resize (n_dofs, n_dofs);
307  for (unsigned int var_i=0; var_i<3; var_i++)
308  for (unsigned int var_j=0; var_j<3; var_j++)
309  Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
310 
311  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
312  {
313  // Row is variable u, v, or w column is x, y, or z
314  TensorValue<Number> grad_u;
315  for (unsigned int var_i=0; var_i<3; var_i++)
316  for (unsigned int var_j=0; var_j<3; var_j++)
317  for (unsigned int j=0; j<n_var_dofs; j++)
318  grad_u(var_i, var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
319 
320  // - C_ijkl u_k,l v_i,j
321  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
322  for (unsigned int i=0; i<3; i++)
323  for (unsigned int j=0; j<3; j++)
324  for (unsigned int k=0; k<3; k++)
325  for (unsigned int l=0; l<3; l++)
326  Re_var[i](dof_i) -= JxW[qp] *
327  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * grad_u(k,l) * dphi[dof_i][qp](j);
328 
329  // assemble \int_Omega C_ijkl u_k,l v_i,j \dx
330  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
331  for (unsigned int dof_j=0; dof_j<n_var_dofs; dof_j++)
332  for (unsigned int i=0; i<3; i++)
333  for (unsigned int j=0; j<3; j++)
334  for (unsigned int k=0; k<3; k++)
335  for (unsigned int l=0; l<3; l++)
336  Ke_var[i][k](dof_i, dof_j) -= JxW[qp] *
337  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * dphi[dof_j][qp](l) * dphi[dof_i][qp](j);
338  }
339 
340  dof_map.constrain_element_matrix_and_vector (Ke, Re, dof_indices);
341 
342  if (jacobian)
343  jacobian->add_matrix (Ke, dof_indices);
344 
345  if (residual)
346  residual->add_vector (Re, dof_indices);
347  }
348 
349  // Move a copy of the mesh based on the solution. This is used to get
350  // the contact forces. We could compute the contact forces based on the
351  // current displacement solution, which would not require a clone of the
352  // mesh. Avoiding the mesh clone would be important for production-scale
353  // contact solves, but for the sake of this example, using the clone is
354  // simple and fast enough.
355  std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
356  move_mesh(*mesh_clone, soln);
357 
358  // Add contributions due to contact penalty forces. Only need to do this on
359  // one processor.
361 
362  for (const auto & [lower_point_id, upper_point_id] : _contact_node_map)
363  {
364  Point upper_to_lower;
365  {
366  Point lower_node_moved = mesh_clone->point(lower_point_id);
367  Point upper_node_moved = mesh_clone->point(upper_point_id);
368 
369  upper_to_lower = (lower_node_moved - upper_node_moved);
370  }
371 
372  // Set the contact force direction. Usually this would be calculated
373  // separately on each master node, but here we just set it to (0, 0, 1)
374  // everywhere.
375  Point contact_force_direction(0., 0., 1.);
376 
377  // gap_function > 0. means that contact has been detected
378  // gap_function < 0. means that we have a gap
379  // This sign convention matches Simo & Laursen (1992).
380  Real gap_function = upper_to_lower * contact_force_direction;
381 
382  // We use the sign of lambda_plus_penalty to determine whether or
383  // not we need to impose a contact force.
384  Real lambda_plus_penalty =
385  (_lambdas[upper_point_id] + gap_function * _contact_penalty);
386 
387  if (lambda_plus_penalty < 0.)
388  lambda_plus_penalty = 0.;
389 
390  // Store lambda_plus_penalty, we'll need to use it later to update _lambdas
391  _lambda_plus_penalty_values[upper_point_id] = lambda_plus_penalty;
392 
393  const Node & lower_node = mesh.node_ref(lower_point_id);
394  const Node & upper_node = mesh.node_ref(upper_point_id);
395 
396  std::vector<dof_id_type> dof_indices_on_lower_node(3);
397  std::vector<dof_id_type> dof_indices_on_upper_node(3);
398  DenseVector<Number> lower_contact_force_vec(3);
399  DenseVector<Number> upper_contact_force_vec(3);
400 
401  for (unsigned int var=0; var<3; var++)
402  {
403  dof_indices_on_lower_node[var] = lower_node.dof_number(_sys.number(), var, 0);
404  lower_contact_force_vec(var) = -lambda_plus_penalty * contact_force_direction(var);
405 
406  dof_indices_on_upper_node[var] = upper_node.dof_number(_sys.number(), var, 0);
407  upper_contact_force_vec(var) = lambda_plus_penalty * contact_force_direction(var);
408  }
409 
410  if (lambda_plus_penalty > 0.)
411  {
412  if (residual && (_sys.comm().rank() == 0))
413  {
414  residual->add_vector (lower_contact_force_vec, dof_indices_on_lower_node);
415  residual->add_vector (upper_contact_force_vec, dof_indices_on_upper_node);
416  }
417 
418  // Add the Jacobian terms due to the contact forces. The lambda contribution
419  // is not relevant here because it doesn't depend on the solution.
420  if (jacobian && (_sys.comm().rank() == 0))
421  for (unsigned int var=0; var<3; var++)
422  for (unsigned int j=0; j<3; j++)
423  {
424  jacobian->add(dof_indices_on_lower_node[var],
425  dof_indices_on_upper_node[j],
426  _contact_penalty * contact_force_direction(j) * contact_force_direction(var));
427 
428  jacobian->add(dof_indices_on_lower_node[var],
429  dof_indices_on_lower_node[j],
430  -_contact_penalty * contact_force_direction(j) * contact_force_direction(var));
431 
432  jacobian->add(dof_indices_on_upper_node[var],
433  dof_indices_on_lower_node[j],
434  _contact_penalty * contact_force_direction(j) * contact_force_direction(var));
435 
436  jacobian->add(dof_indices_on_upper_node[var],
437  dof_indices_on_upper_node[j],
438  -_contact_penalty * contact_force_direction(j) * contact_force_direction(var));
439  }
440  }
441  else
442  {
443  // We add zeros to the matrix even when lambda_plus_penalty = 0.
444  // We do this because some linear algebra libraries (e.g. PETSc)
445  // will condense out any unused entries from the sparsity pattern,
446  // so adding these zeros in ensures that these entries are not
447  // condensed out.
448  if (jacobian && (_sys.comm().rank() == 0))
449  for (unsigned int var=0; var<3; var++)
450  for (unsigned int j=0; j<3; j++)
451  {
452  jacobian->add(dof_indices_on_lower_node[var],
453  dof_indices_on_upper_node[j],
454  0.);
455 
456  jacobian->add(dof_indices_on_lower_node[var],
457  dof_indices_on_lower_node[j],
458  0.);
459 
460  jacobian->add(dof_indices_on_upper_node[var],
461  dof_indices_on_lower_node[j],
462  0.);
463 
464  jacobian->add(dof_indices_on_upper_node[var],
465  dof_indices_on_upper_node[j],
466  0.);
467  }
468  }
469  }
470 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
This is the EquationSystems class.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1025
A Node is like a Point, but with more information.
Definition: node.h:52
std::map< dof_id_type, Real > _lambda_plus_penalty_values
Store the intermediate values of lambda plus penalty.
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:2254
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:1992
unsigned int dim
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:374
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]...
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2144
const EquationSystems & get_equation_systems() const
Definition: system.h:730
MeshBase & mesh
processor_id_type rank() const
Real _contact_penalty
Penalize overlapping elements.
const Parallel::Communicator & comm() const
Order default_quadrature_order() const
Definition: fe_type.h:357
const MeshBase & get_mesh() const
Definition: system.h:2277
virtual std::unique_ptr< MeshBase > clone() const =0
Virtual "copy constructor".
This is the MeshBase class.
Definition: mesh_base.h:74
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.
unsigned int variable_number(std::string_view var) const
Definition: system.C:1557
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
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:169
unsigned int number() const
Definition: system.h:2269
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 T & get(std::string_view) const
Definition: parameters.h:426
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
virtual void zero()=0
Set all entries to 0.
std::map< dof_id_type, Real > _lambdas
Augmented Lagrangian values at each contact node.
Defines a dense submatrix for use in Finite Element-type computations.
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.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< dof_id_type, dof_id_type > _contact_node_map
This provides a map between contact nodes.
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:324
Parameters parameters
Data structure holding arbitrary parameters.
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:575
const DofMap & get_dof_map() const
Definition: system.h:2293
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.

◆ set_contact_penalty()

void LinearElasticityWithContact::set_contact_penalty ( Real  contact_penalty_in)

Update the penalty parameter.

Definition at line 42 of file linear_elasticity_with_contact.C.

References _contact_penalty.

43 {
44  _contact_penalty = contact_penalty_in;
45 }
Real _contact_penalty
Penalize overlapping elements.

◆ update_lambdas()

std::pair< Real, Real > LinearElasticityWithContact::update_lambdas ( )

Update the lambda parameters in the augmented Lagrangian method.

Returns
the largest change in the lambdas, and the largest lambda value.

Definition at line 598 of file linear_elasticity_with_contact.C.

References _lambda_plus_penalty_values, _lambdas, std::abs(), and libMesh::Real.

599 {
600  Real max_delta_lambda = 0.;
601  Real max_new_lambda = 0.;
602 
603  for (auto & [upper_node_id, lambda_val] : _lambdas)
604  {
605  auto new_lambda_it = _lambda_plus_penalty_values.find(upper_node_id);
606  libmesh_error_msg_if(new_lambda_it == _lambda_plus_penalty_values.end(), "New lambda value not found");
607 
608  Real new_lambda = new_lambda_it->second;
609  Real old_lambda = lambda_val;
610 
611  lambda_val = new_lambda;
612 
613  Real delta_lambda = std::abs(new_lambda-old_lambda);
614  if (delta_lambda > max_delta_lambda)
615  max_delta_lambda = delta_lambda;
616 
617  if (std::abs(new_lambda) > max_new_lambda)
618  max_new_lambda = std::abs(new_lambda);
619  }
620 
621  return std::make_pair(max_delta_lambda, max_new_lambda);
622 }
std::map< dof_id_type, Real > _lambda_plus_penalty_values
Store the intermediate values of lambda plus penalty.
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:57
std::map< dof_id_type, Real > _lambdas
Augmented Lagrangian values at each contact node.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

Member Data Documentation

◆ _contact_node_map

std::map<dof_id_type, dof_id_type> LinearElasticityWithContact::_contact_node_map
private

This provides a map between contact nodes.

Definition at line 61 of file linear_elasticity_with_contact.h.

Referenced by add_contact_edge_elements(), get_least_and_max_gap_function(), initialize_contact_load_paths(), and residual_and_jacobian().

◆ _contact_penalty

Real LinearElasticityWithContact::_contact_penalty
private

Penalize overlapping elements.

Definition at line 44 of file linear_elasticity_with_contact.h.

Referenced by get_contact_penalty(), residual_and_jacobian(), and set_contact_penalty().

◆ _lambda_plus_penalty_values

std::map<dof_id_type, Real> LinearElasticityWithContact::_lambda_plus_penalty_values
private

Store the intermediate values of lambda plus penalty.

The dof IDs refer to the nodes on the upper contact surface.

Definition at line 50 of file linear_elasticity_with_contact.h.

Referenced by residual_and_jacobian(), and update_lambdas().

◆ _lambdas

std::map<dof_id_type, Real> LinearElasticityWithContact::_lambdas
private

Augmented Lagrangian values at each contact node.

The dof IDs refer to the nodes on the upper contact surface.

Definition at line 56 of file linear_elasticity_with_contact.h.

Referenced by initialize_contact_load_paths(), residual_and_jacobian(), and update_lambdas().

◆ _sys

NonlinearImplicitSystem& LinearElasticityWithContact::_sys
private

The documentation for this class was generated from the following files: