 ## Solid and Structural Mechanics Demo 3: A unilateral contact problem

This code solves an elastostatic problem involving unilateral contact.

• The code essentially looks like lelas2d. The only difference is that since the problem is nonlinear, due to the contact condition, an iterative procedure must be supplied. We describe here the part of the code that performs iterations:
• We make use in this code of some global variables of to monitor iterations. These variables are MaxNbIterations, theIteration, MaxNbIterations, theTolerance and theDiscrepancy with obvious meaning. Note that all these variables have default values.

 ``` MaxNbIterations = data.getNbIter(); theTolerance = data.getTolerance(); ```

• The iteration loop can be drived by the macro IterationLoop which acts like a loop and which updates iteration index and stops iterations when the value of theDiscrepancy is under theTolerance or when theIteration reaches the value of MaxNbIterations.

 ``` IterationLoop { v = 0.; A = 0.; MeshElements(ms) { Elas2DT3 eq(theElement); eq.Deviator(); eq.Dilatation(); eq.BodyRHS(body_f); eq.ElementAssembly(A); eq.ElementAssembly(v); } MeshSides(ms) { Elas2DT3 eq(theSide,u); eq.BoundaryRHS(bound_f); eq.SignoriniContact(bound_f,penal); eq.SideAssembly(A); eq.SideAssembly(v); } A.Prescribe(v,bc); A.solve(v); theDiscrepancy = Discrepancy(u,v,2,1); if (verbose > 0) cout << "Iteration: " << it << ". Discrepancy: " << theDiscrepancy << endl; } ```

Let us give some comments on the source code here above:

• The contact condition is implemented via a penalty formulation as a boundary traction. For this, we make a loop over sides. The vector bound_f contains the gap at nodes (the distance between boundary nodes and the obstacle. Note that the only concerned sides are those who have a code equal to CONTACT_A or CONTACT_B
• We have used global variables of to monitor iterations. These variables are theIteration

## An example

Let us run this program with the data presented in the following project file:

 ``` 1 1 1 1 1 beam.m beam.d beam.pr ```

The only thing that may differ from other seen project files is that prescription data is given in a separate file: beam.pr. Let us give this one:

 ``` 0.0 -1000.0 ```

Clearly, we have prescribed by this file a (Dirichlet) boundary condition that enforces a value of 0.0 for degrees of freedom with code 1. In addition, we have prescribed a traction of value -1000 for the second degrees of freedom (y-traction) having code 1.