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. 
