2D Elasticity code in fenicsx does not match the code from fenics

dolfinx
1

It would be great to have some advice on this code which confirms numerically the theoretical maximum deflection.

I’m trying to write the code in fenicsx, but the result does not confirm the theoretical deflections as the the preivous fenics code.

Here is my version of the code in fenicsx

import numpy as np
import dolfinx
from dolfinx import mesh
import ufl
from petsc4py import PETSc
from mpi4py import MPI
import numpy


# --------------------
# Domain
# --------------------


d_x0 = 25.0  #  size along x0, m
d_x1 = 1.0  #  size along x1, m

n_x0 = 250
n_x1 = 10
domain = mesh.create_rectangle(
    MPI.COMM_WORLD,
    [[0, 0], [d_x0, d_x1]],
    [n_x0, n_x1])



from dolfinx import fem

# --------------------
# Parameters
# --------------------

# Density
rho = 1.0e-3 #fem.Constant(domain, 1.0)

# Young's modulus and Poisson's ratio
E = 1.0e5
nu = 0.3

# Lame's constants
mu = E/(2*(1.0+nu))
lambda_ = E*nu/((1.0+nu)*(1.0-2.0*nu))



model = "plane_stress"

if model == "plane_stress":
    lambda_ = 2.0*mu*lambda_/(lambda_+2.0*mu)




from dolfinx.fem import FunctionSpace
#V = fem.VectorFunctionSpace(domain, ("Lagrange", 2))
element = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2)
V = FunctionSpace(domain, element)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)


# --------------------
# Boundary Functions
# --------------------
tdim = domain.topology.dim
fdim = tdim - 1


def bottom(x):
    return np.isclose(x[1], 0)
def top(x):
    return np.isclose(x[1], 1)
def left(x):
    return np.isclose(x[0], 0)
def right(x):
    return np.isclose(x[0], 25)


bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
bottom_values = np.full_like(bottom_facets, 1)
top_facets = mesh.locate_entities_boundary(domain, fdim, top)
top_values = np.full_like(top_facets, 2)
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
left_values = np.full_like(left_facets, 3)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
right_values = np.full_like(right_facets, 4)

facets = np.hstack([bottom_facets,top_facets, left_facets, right_facets])
values = np.hstack([bottom_values,top_values, left_values, right_values])
sort_order = np.argsort(facets)

facet_tags = dolfinx.mesh.meshtags(domain, fdim, facets[sort_order], values[sort_order])

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
dx = ufl.Measure("dx", domain=domain)

u_bc = fem.Function(V)
u_bc.interpolate(lambda x: ([x[0]*0,x[1]*0]))
#u_bc=fem.Constant(domain,(0.0, 0.0))
bcl = fem.dirichletbc(u_bc, left_facets)
bcr = fem.dirichletbc(u_bc, right_facets)


# Load

g_z = 0.0
b_z = -1.0

g  = fem.Constant(domain,(0,g_z))
b  = fem.Constant(domain,(0,b_z))

# ------
# Forms
# ------

def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


def sigma(u):
    #return lambda_ * ufl.tr(epsilon(u)) * ufl.Identity(2) + 2 * mu * epsilon(u)
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(2) + 2.0 * mu * epsilon(u)


a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = rho*ufl.dot(b, v)*ufl.dx + ufl.dot(g, v)*ds((2))



uh = dolfinx.fem.Function(V)
from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a, L, bcs=[bcl], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

print("Max vertical deflection",np.max(-uh.sub(1).x.array))
print("Theoretical vertical deflection",3*(rho)*d_x0**4/2/E/d_x1**3)