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