Hey dokken, sorry for bothering you again. I would have thought I have solved the problem with your help. But I meet problem again when I try to solve this linear elasticity problem:

I guess that the solution is

That solution satisfies the displacement boundary condition. Besides, we can calculate that

So, it also satisfies the balance equation and force boundary condition:

If this solution is right, the displacement at point(L,L) should be

But the result I get from FEniCSx is

```
array([ 0.01343142, -0.1 ])
```

Could you please help find whether there is any error in my calculation or code?

```
import numpy as np
from mpi4py import MPI
from ufl import Identity, Measure, TestFunction, TrialFunction, VectorElement, dot, dx, inner, grad, nabla_div, sym
from dolfinx.generation import RectangleMesh, UnitSquareMesh
from dolfinx.mesh import CellType, locate_entities_boundary
from dolfinx.fem import (Constant, DirichletBC, Function, LinearProblem, FunctionSpace, VectorFunctionSpace,
locate_dofs_geometrical)
L = 1
pixel = 20
mesh = RectangleMesh(MPI.COMM_WORLD, np.array([[0, 0, 0], [L, L, 0]]), [pixel, pixel], cell_type=CellType.quadrilateral)
V = VectorFunctionSpace(mesh, ("CG", 1))
nu = 0.2
mu_m = 1
mu_f = 5*mu_m
import ufl
x = ufl.SpatialCoordinate(mesh)
mu = ufl.exp(mu_m + 4*x[0]*(1-x[0])*(mu_f - mu_m))
lambda_ = ufl.exp((mu_m + 4*x[0]*(1-x[0])*(mu_f - mu_m))*2*nu/(1-2*nu))
#define boundary conditions
Vy = V.sub(1).collapse()
uDy_high = Function(Vy)
with uDy_high.vector.localForm() as uDy_high_loc:
uDy_high_loc.set(-0.1*L)
def highside(x):
return np.isclose(x[1], L)
boundary_high = locate_dofs_geometrical((V.sub(1), Vy), highside)
bch = DirichletBC(uDy_high, boundary_high, V.sub(1))
uDy_low = Function(Vy)
with uDy_low.vector.localForm() as uDy_low_loc:
uDy_low_loc.set(0)
def lowside(x):
return np.isclose(x[1], 0)
boundary_low = locate_dofs_geometrical((V.sub(1), Vy), lowside)
bcl = DirichletBC(uDy_low, boundary_low, V.sub(1))
Vx = V.sub(0).collapse()
uDx = Function(Vx)
with uDx.vector.localForm() as uDx_loc:
uDx_loc.set(0)
def origin(x):
return np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0))
boundary_origin = locate_dofs_geometrical((V.sub(0), Vx), origin)
bco = DirichletBC(uDx, boundary_origin, V.sub(0))
u_D = Function(V)
with u_D.vector.localForm() as loc:
loc.set(0)
bco = DirichletBC(u_D, locate_dofs_geometrical(V, origin))
bcs = [bco, bcl, bch]
def epsilon(u):
return sym(grad(u))
def sigma(u):
return lambda_ * nabla_div(u) * Identity(u.geometric_dimension()) + 2*mu*epsilon(u)
u = TrialFunction(V)
v = TestFunction(V)
from petsc4py import PETSc
f = Constant(mesh, PETSc.ScalarType((0, 0)))
T = Constant(mesh, PETSc.ScalarType((0, 0)))
ds = Measure("ds", domain=mesh)
a = inner(sigma(u), epsilon(v)) * dx
L = inner(f, v) * dx + inner(T, v) * ds
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
import pyvista
import dolfinx.plot
# Start virtual framebuffer
pyvista.start_xvfb(wait=0.05)
# Create plotter and pyvista grid
p = pyvista.Plotter(title="Deflection", window_size=[800, 800])
topology, cell_types = dolfinx.plot.create_vtk_topology(mesh, mesh.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, mesh.geometry.x)
# Attach vector values to grid and warp grid by vector
vals_2D = uh.compute_point_values().real
vals = np.zeros((vals_2D.shape[0], 3))
vals[:,:2] = vals_2D
grid["u"] = vals
actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1)
actor_1 = p.add_mesh(warped, show_edges=True)
p.view_xy()
if not pyvista.OFF_SCREEN:
p.show()
else:
fig_array = p.screenshot(f"component.png")
end = pixel + 1
vals_2D[end*end-1, :]
```

Thanks very much!