Contact Mechanics - 2D problem


I am trying to solve a 2D contact problem based on the 3D code available here:

But I can’t figure out why my code is not running. My code is the following

Ny=30 # Refinamento da malha na direção y

Nx=30 # Refinamento da malha na direção x

mesh = fn.RectangleMesh(fn.Point(0, -1), fn.Point(1, 0),Nx,Ny) # gera a malha 



fn.plot(mesh, title="Malha")

class Top(fn.SubDomain):

    def inside(self, x, on_boundary):

        return fn.near(x[1], 0.) and on_boundary

def symmetry_x(x, on_boundary):

        return fn.near(x[0], 0.) and on_boundary

def bottom(x, on_boundary):

        return fn.near(x[1], -1.) and on_boundary

# exterior facets MeshFunction

facets = fn.MeshFunction("size_t", mesh, 1)


Top().mark(facets, 1)

ds = fn.Measure('ds', subdomain_data=facets)

R = 0.5

d = 0.02

obstacle = fn.Expression("-d+pow(x[0],2)/2/R", d=d, R=R, degree=2)

V = fn.VectorFunctionSpace(mesh, "CG", 1)

V2 = fn.FunctionSpace(mesh, "CG", 1)

V0 = fn.FunctionSpace(mesh, "DG", 0)

u = fn.Function(V, name="Displacement")

du = fn.TrialFunction(V)

u_ = fn.TestFunction(V)

bc =[fn.DirichletBC(V, fn.Constant((0., 0.)), bottom),

     fn.DirichletBC(V.sub(0), fn.Constant(0.), symmetry_x)]

E = fn.Constant(10.)

nu = fn.Constant(0.3)

mu = E/2/(1+nu)

lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):

    return fn.sym(fn.grad(v))

def sigma(v):

    return lmbda**fn.Identity(2) + 2.0*mu*eps(v)

def ppos(x):

    return (x+np.abs(x))/2.

pen = fn.Constant(1e4)

form = fn.inner(sigma(u), eps(u_))*fn.dx + pen*[1], ppos(u[1]-obstacle))*fn.ds(1)

J = fn.derivative(form, u, du)

problem = fn.NonlinearVariationalProblem(form, u, bc, J=J)

solver = fn.NonlinearVariationalSolver(problem)

solver.parameters["newton_solver"]["linear_solver"] = "cg"

solver.parameters["newton_solver"]["preconditioner"] = "ilu"


import matplotlib.pyplot as plt

p = fn.plot(sigma(u)[1,1], mode='color') 



I tried different nonlinear solvers but the result is the same. What is wrong?

Could you please supply the error message you are getting, as well as specifying what version of fenics you are using.

The difference comes from the fact that the tutorial is using a wildcard import (i.e., from dolfin import *). The problem is that you set subdomain_data in ds, but then use fn.ds in form. You can either make the change

fn.ds = fn.Measure('ds', subdomain_data=facets)

(where I’ve added an fn. at the beginning) or use ds(1) (without the fn.) to integrate the boundary term.


That’s true. Thank you so much. It solved my problem. Now it works.


Hey, how would one model the indenter explicitly instead of using the Expression without any mesh?

I am interested in this problem as well. You need to discretize both bodies and consider some penalty/augmented Lagrange formulation in order to impose displacement constraints and force equilibrium on the contact surfaces of both bodies. I don’t believe that there is any example available in the web and I don’t know if it is possible to solve this problem with Fenics yet. The best discussion that I found in the web was this:

It was from many years ago (2013). At that time, to solve this problem was very difficult or even impossible with Fenics. As far as I know, few things have changed since then.

I discretized both bodies and also marked them differently as done in the fenics tutorial. My idea was to define a distance function between the lower boundary of the upper geometry and the upper boundary of the bottom geometry and use it for the penalty approach. But I am not sure how to do that. My problem is that I do not understand how to use the boundary values of the two boundaries as a function input.

I suggest you to have a look at this link:

It is not Fenics but it is well documented and has a code with the implementation. It is the best finite element code applied to contact mechanics that I found. I am trying to follow a similar approach in order to implement it in Fenics. But I am still in the beginning.