Linear elasticity example with unit disc mesh

Question

Hi,

I was attempting to redo the Linear Elasticity example found here with a unit disc mesh.

the code with the disc mesh:

---

# variables
mu = 1
rho = 1
beta = 1.25
lambda_ = beta
g = 9.81

# create mesh
comm = mpi_comm_world()
mesh = UnitDiscMesh(comm, 20, 2, 3) 
# create function space
V = VectorFunctionSpace(mesh, 'Lagrange', 3) # defining a vector valued function space 
over the mesh with a lagrangian finite elements of degree 1

 # boundary condition
 tol = 1E-14

def clamped_boundary(x, on_boundary):
     return on_boundary and x[0] < 0 # this should "clamp" the part of the disc left of the x axis

bc = DirichletBC(V, (0, 0, 0), clamped_boundary) # create dirichlet boundary condition 
where u = (0,0,0) on the clamped boundary

# define stress and strain
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T) # define the symmetric part of the gradient of a vector function u(a tensor)

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u) # defining the function for 
sigma dependent on the displacement u

# define trial and test functions
u = TrialFunction(V)
d = u.geometric_dimension() # not sure what this is doing
v = TestFunction(V) #defined both the test and trial function over the function space V

# variational problem
f = Constant((0, 0, -rho*g)) # force per unit body mass
T = Constant((0, 0, 0)) # sigma dot n
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds # variational equation

# compute solution
u = Function(V)
solve(a == L, u, bc, solver_parameters={'linear_solver':'mumps'})

# save solution for paraview
vtkfile = File('LE-solution.pvd')
vtkfile << u

---

I was hoping to get something similar to the example above but instead of having a box clamped at x = 0, have the unit disc clamped in the region where the x coordinate is negative.
instead I get this

I'm not exactly sure what I'm doing wrong here.
Any help would be appreciated.

Comments

Please indent your code snippets with four spaces, leaving newlines before and after it, so that it displays correctly, like this:

from dolfin import *
V = FunctionSpace(UnitSquareMesh(1,1), "CG", 2)
# ... 

Answer 1

I've copied and pasted your example and I see a disk clamped on the half of its boundary with x[0] < 0 as expected:

Maybe you are running some super-old and buggy version of FEniCS?

Comments

Hi,

Thanks for the reply.

I'm running FEniCS on windows through docker.
I tried to pull the latest stable release a second ago, and its telling me I already have it.

Its not paraview either because when I view the original linear elasticity example it shows what you'd expect.

---

I guess you are applying the "warp by vector" filter in ParaView to see the solution. Maybe the scaling factor is too large?