1 #ifndef LINEAR_ELASTICITY_WITH_CONTACT_H
2 #define LINEAR_ELASTICITY_WITH_CONTACT_H
3 
4 #include "libmesh/dof_map.h"
5 #include "libmesh/nonlinear_implicit_system.h"
6 
7 // define the subdomain IDs
8 #define TOP_SUBDOMAIN 2
9 #define BOTTOM_SUBDOMAIN 1
10 
11 // define the boundary IDs in the mesh
12 #define MIN_Z_BOUNDARY 1
13 #define MAX_Z_BOUNDARY 2
14 #define CONTACT_BOUNDARY_LOWER 3
15 #define CONTACT_BOUNDARY_UPPER 4
16 
17 using libMesh::DofMap;
18 using libMesh::NonlinearImplicitSystem;
19 using libMesh::dof_id_type;
20 using libMesh::Point;
21 using libMesh::Real;
22 using libMesh::Number;
23 using libMesh::MeshBase;
24 using libMesh::NumericVector;
25 using libMesh::SparseMatrix;
26 
27 /**
28  * This class encapsulate all functionality required for assembling
29  * and solving a linear elastic model with contact.
30  */
31 class LinearElasticityWithContact :
32   public NonlinearImplicitSystem::ComputeResidualandJacobian
33 {
34 private:
35 
36   /**
37    * Keep a reference to the NonlinearImplicitSystem.
38    */
39   NonlinearImplicitSystem & _sys;
40 
41   /**
42    * Penalize overlapping elements.
43    */
44   Real _contact_penalty;
45 
46   /**
47    * Store the intermediate values of lambda plus penalty. The dof IDs refer
48    * to the nodes on the upper contact surface.
49    */
50   std::map<dof_id_type, Real> _lambda_plus_penalty_values;
51 
52   /**
53    * Augmented Lagrangian values at each contact node. The dof IDs refer
54    * to the nodes on the upper contact surface.
55    */
56   std::map<dof_id_type, Real> _lambdas;
57 
58   /**
59    * This provides a map between contact nodes.
60    */
61   std::map<dof_id_type, dof_id_type> _contact_node_map;
62 
63 public:
64 
65   /**
66    * Constructor.
67    */
68   LinearElasticityWithContact(NonlinearImplicitSystem & sys_in,
69                               Real contact_penalty_in);
70 
71   /**
72    * Update the penalty parameter.
73    */
74   void set_contact_penalty(Real contact_penalty_in);
75 
76   /**
77    * Get the penalty parameter.
78    */
79   Real get_contact_penalty() const;
80 
81   /**
82    * Kronecker delta function.
83    */
84   Real kronecker_delta(unsigned int i,
85                        unsigned int j);
86 
87   /**
88    * Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
89    */
90   Real elasticity_tensor(Real young_modulus,
91                          Real poisson_ratio,
92                          unsigned int i,
93                          unsigned int j,
94                          unsigned int k,
95                          unsigned int l);
96 
97   /**
98    * Move the mesh nodes of input_mesh based on the displacement field
99    * in input_solution.
100    */
101   void move_mesh(MeshBase & input_mesh,
102                  const NumericVector<Number> & input_solution);
103 
104   /**
105    * Set up the load paths on the contact surfaces.
106    */
107   void initialize_contact_load_paths();
108 
109   /**
110    * Add edge elements into the mesh based on the contact load paths.
111    * This ensure proper parallel communication of data, e.g. if a node
112    * on one side of the contact surface has a constraint on it, then
113    * adding contact elements into the mesh ensures that the constraint
114    * will be enforced properly.
115    */
116   void add_contact_edge_elements();
117 
118   /**
119    * Evaluate the Jacobian of the nonlinear system.
120    */
121   virtual void residual_and_jacobian (const NumericVector<Number> & soln,
122                                       NumericVector<Number> * residual,
123                                       SparseMatrix<Number> * jacobian,
124                                       NonlinearImplicitSystem & /*sys*/);
125 
126   /**
127    * Compute the Cauchy stress for the current solution.
128    */
129   void compute_stresses();
130 
131   /**
132    * Update the lambda parameters in the augmented Lagrangian
133    * method.
134    * @return the largest change in the lambdas, and the largest
135    * lambda value.
136    */
137   std::pair<Real, Real> update_lambdas();
138 
139   /**
140    * @return the least and max gap function values for the current solution.
141    */
142   std::pair<Real, Real> get_least_and_max_gap_function();
143 };
144 
145 #endif
146