hello. It’s been a while since I’ve posted.

I recently posted to address a question I had while exploring legacy DOLFIN a bit more.

I’m having an issue when I’m trying to get the strain energy in the linear elasticity example provided by this link, and when I add the equation defined by `UFL Form`

to the entire mesh domain with both the result obtained by `assemble()`

and the vector element values obtained by `project()`

, the values are different, and I’m wondering why the difference occurs and how I can reduce the difference.

In the FEA process, I only changed the parameter values in the example code as shown below.

(3D to 2D, number of mesh elements and `f`

adjustment)

```
from fenics import *
# Scaled variables
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma
# Create mesh and define function space
mesh = RectangleMesh(Point(0, 0), Point(L, W), 2, 2)
V = VectorFunctionSpace(mesh, 'P', 1)
# Define boundary condition
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < 1e-14
bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)
# Define strain and stress
def epsilon(u):
return 0.5*(grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u):
return lambda_*div(u)*Identity(d) + 2*mu*epsilon(u)
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, -100)) # f = Constant((0, -rho*g))
T = Constant((0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds
# Compute solution
u = Function(V)
solve(a == L, u, bc)
```

And the strain energy function expression was defined in `UFL format`

as follows.

Strain energy equation: \psi(\mathbf{u})=\frac{\lambda}{2}\{\textrm{tr}(\epsilon(\mathbf{u}))\}^{2}\ + \mu\textrm{tr}(\epsilon(\mathbf{u})^{2})

(\epsilon(\mathbf{u})=\frac{1}{2}(\nabla\mathbf{u}\ + (\nabla\mathbf{u})^{T})=\rm{sym}(\nabla\mathbf{u}))

```
def strain_energy(u):
return 0.5*lambda_*(tr(epsilon(u))**2) + mu*tr(epsilon(u)*epsilon(u))
```

When trying to find the total strain energy value for a given mesh domain, the result obtained with `assemble()`

and `"dx"`

, (type: `ufl.algebra.Sum`

) was different from the result obtained when `sum()`

of the vector results `project()`

in `DG0`

space as shown below.

**(I’d actually like to see the UFL Form data, but it’s hidden, so I can’t…)**

```
print(f"type of strain energy: {type(strain_energy(u))}")
# same as print(f"strain energy using assemble() quad: {assemble(strain_energy(u)*Measure("dx", mesh, metadata={"quadrature_degree": 0}))}\n")
print(f"strain energy using assemble(): {assemble(strain_energy(u)*dx)}\n")
print(f"strain energy using project() to DG0: {sum(project(strain_energy(u), FunctionSpace(mesh, "DG", 0)).vector()[:])}")
type of strain energy: <class 'ufl.algebra.Sum'>
strain energy using assemble(): -0.008689827046779201
strain energy using project() to DG0: -0.3475930818711679
```

(In the case of dx used in `assemble()`

, there is a difference even after adjusting the `"quadrature_degree"`

of `metadata`

, so my guess is that it’s an interpolation error when `project()`

with `DG0`

.)

**My questions are**

A) What mesh points(or functionspace?!) is that output obtained with `assemble()`

using `"dx"`

based on? I would like to know the underlying reason for the difference between these two values. Am I correct in my guess above as to why the error occurs?

B) If my guess is correct, I believe that the result obtained with `assemble()`

and `"dx"`

is more reliable value. If so, I would like to get the result obtained using `project()`

: (value: `-0.3475930818711679`

) to be the same value as the result using `assemble()`

: (value: `-0.008689827046779201`

), how can I do that?

**(In fact, I wanted to get the strain energy vector array term corresponding to each element).**

Any opinions are welcome. I can’t resolve my question.