Plane Stress Problem

Hi.

I am learning to use Fenics and after completing the basic tutorials, I am trying to solve a plane stress problem. After searching for some clues in this forum, I can now make it work, but with a wrong solution when comparing to an analytical solution.

The test case is a cantilever beam with a constant distributed load on the top. The analytical solution is known.

Thanks for your help !

My actual code:

import numpy as np
from  dolfin import *

def main():

    # Length, heigth and width
    L = 1.0
    H = 0.1
    W = 0.01

    # Number of divisions 
    NX = 100
    NY = 10

    # Young modulus. Null Poisson in order to compare to the analytical solution (beam)
    Ex = 210.0E9
    nu = 0.0

    # Distributed load (N/m)
    q = -100.0
    T = Constant((0.0,q))


    # Analytical solution
    I = (W*H**3)/12.0
    delta_max = (q*L**4)/(8.0*Ex*I)


    # Mesh 
    mesh = RectangleMesh(Point(0.0, 0.0), Point(L, H), NX, NY, "right/left")

    Degree = 1
    V = VectorFunctionSpace(mesh, "Lagrange", Degree, dim=2)
    u = TrialFunction(V)
    v = TestFunction(V)


    # Constitutive relation for the Plane Stress case
    c1 = Ex/(1.0-nu**2.0)
    c2 = (1.0-nu)/2.0
    C_numpy = c1* np.array([[1.0 ,  nu   , 0.0],
                           [ nu  ,  1.0  , 0.0],
                           [0.0  ,  0.0  , c2 ]])
    C = as_matrix(C_numpy)

    # Dirichlet boundary conditions
    def Engaste(x, on_boundary):
        tol = 1E-14
        return on_boundary and abs(x[0]-0.0) < tol

    bc = DirichletBC(V, Constant((0, 0)), Engaste)


    # Natural boundary condition
    boundary_markers = FacetFunction('size_t', mesh)

    ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

    class Face_Superior(SubDomain):
        def inside(self, x, on_boundary):
            tol = 1E-14 
            return on_boundary and near(x[1], H, tol) 

    by = Face_Superior()
    by.mark(boundary_markers, 0)

    # Strain - Voigt notation 
    def eps(u):
        """ Returns a vector of strains of size (3,1) in the Voigt notation
        layout {eps_xx, eps_yy, gamma_xy} where gamma_xy = 2*eps_xy"""
        return as_vector([u[i].dx(i) for i in range(2)] +
                         [u[i].dx(j) + u[j].dx(i) for (i,j) in [(0,1)]])


    # Equilibrium problem
    LHS = W*inner(eps(u), C*eps(v))*dx
    RHS = dot(T, v)*ds(0)

    # Solution
    U = Function(V)
    problema = LinearVariationalProblem(LHS, RHS, U, bcs=bc)
    solver = LinearVariationalSolver(problema)
    solver.parameters["symmetric"] = True
    solver.solve()

    print("Numerical solution: ", U( (L, 0.0) ) , " analytical solution: ", delta_max)



if __name__ == "__main__":
    main()