Element-wise strain energy

1

I’m working with structures under plane stress state and for that I have to solve a linear elasticity problem. I’m currently able to solve the PDE for the displacement and then compute the structure’s compliance (strain energy), which is a real number (scalar). The problem is that for the post processing I’m intersted, I need a element or nodal value of energy and I’m having trouble doing that.

When solving the problem without using FEniCS, It works fine, because I actually find a element wise value and then sum them all for the structures compliance. However, since FEniCS returns me the “sum” directly, I don’t know how to proceed. I tried projecting the value on a mesh, but I don’t know if that is correct. Could anyone please help me?

This is how I am currently getting the compliance value (assuming that U is the structure’s displacement, soluting of the PDE):

eU = sym(grad(U))
S1 = 2.0*mu*inner(eU,eU) + lmbda*tr(eU)**2 
compliance = assemble( eps_er*S1* dx(0) + S1*dx(1)) ##returns me a scalar

I have two components of dx because I’m working with a rectangular domain and the sctructure is a subset of this domain. Then, dx(1) refers to the inside of the structure and eps_er is a small value associate with the complementary part.

2

You can get the cell-wise average of a quantity f by project-ing it onto FunctionSpace(mesh,"DG",0), i.e., the space of constant functions on each element. (The project function is just a short-hand for L^2 projection.) To get the integral of the quantity over each element (instead of the average), you can project CellVolume(mesh)*f, to cancel the 1.0/CellVolume(mesh) in the average.

Some other notes:

  • The formula for compliance in the given code snippet is twice the strain energy (but of course a constant factor won’t affect the minimizer, if compliance is being used as an objective function in optimization).

  • If mu and lmbda are the 3D Lamé parameters, the formula for S1 is (twice) the energy density in plane strain when applied in a 2D domain, not plane stress. (Consider, e.g., that tr(eU) is missing the zz component of the strain tensor, which is nonzero in plane stress.) To convert a plane strain solution to a plane stress solution, one must replace the material parameters. In terms of Young’s modulus and Poisson ratio, the plane strain \to plane stress substitutions are

    E~\to~\frac{(1+2\nu)E}{(1+\nu)^2}~~~\text{and}~~~\nu~\to~\frac{\nu}{1+\nu}\text{ .}

    (An equivalent transformation of \mu and \lambda can be obtained after some additional algebra, but the above are more readily available in references.)

2 Likes
3

Thank you for the reply! Also, thank you for the notes, they are very helpful since I’m new to the area. I tried running both functions you suggested, but I still have a question.

For ilustration, I computed the compliance for an exemple structure and got 3768463. However, when I project S1 on FunctionSpace(mesh,"DG",0) and use .vector() to see the values, they are all 3768463 (the same as the total compliance!). I expected to get smaller values that summed up 3768463. Also, I expected each element to have a different compliance value, since they refer to different regions of the structure, that have different displacement values.

This is what I’m doing for the projection:

V1 =FunctionSpace(mesh,"DG",0)
eleComp = project(S1, V1)

I also tried using CellVolume(mesh)*S1, as you sugested, but then I get an object of this type:

Product(FloatValue(3768463.419263469), CellVolume(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 109)))

for which I cannot acess .vector() to actually see the values.

If its relevant, for solving the equilibrium problem, I’m working with this space:

V = FunctionSpace(mesh, 'CG', 1)

This seems to be correct, because at least in terms of dimensions I’m getting something right: when I project I get a vector with number of components equal to the number of elements and the solution U in vector form has a number of componentes equal to the degrees of freedom of the problem (twice the number of nodes).

Am I doing something wrong? Thank you so much for the help!

5

For ilustration, I computed the compliance for an exemple structure and got 3768463. However, when I project S1 on FunctionSpace(mesh,"DG",0) and use .vector() to see the values, they are all 3768463 (the same as the total compliance!). I expected to get smaller values that summed up 3768463.

This would be consistent with the projection giving the average energy on each cell if the domain has unit measure (e.g., from UnitSquareMesh) and the example problem has a uniform strain field.

I also tried using CellVolume(mesh)*S1 , as you sugested, but then I get an object of this type:

What I was suggesting was to include the scaling within the projection, i.e.,

eleComp = project(CellVolume(mesh)*S1,V1)

The cell-wise average of a quantity scaled by cell volume will be its integral over each cell.

1 Like
6

Thank you for the reply and for your patience. It works now!
One last question: what would be the equivalent expression/correct space for projecting in each node instead of element?

7

please i don’t know if you can also help me because i have almost the same problem just that in my case it is anisometric and when i use the dot function i get the errors. i am really stuck on this point and i would be very happy if you can unblock me.

8

Defining energy per element is more natural, especially for constant-strain elements. You could come up with some nodal proxy for it in the space FunctionSpace(mesh,"CG",1) (which has a degree of freedom at each mesh vertex), but the details would depend on how exactly you want to define the projection.

It occurred to me (from the language of the post), though, that perhaps you actually want nodal compliance (i.e., the work done by external forces), not nodal (or elemental) contributions to strain energy. Distributions of compliance and strain energy can actually be quite different locally, even if they both sum to the same total value. (E.g., the density of work done by a distributed transverse load on a cantilever beam will be greatest at the free end, while the strain energy density will be concentrated toward the support.) The nodal compliance values can be obtained from the algebraic system, as demonstrated in the following example using the scalar Poisson equation:

from dolfin import *

# Set up and solve a simple Poisson problem:
N = 32
#mesh = UnitSquareMesh(N,N)
mesh = UnitIntervalMesh(N)
V = FunctionSpace(mesh,"CG",1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u),grad(v))*dx
L = v*dx
uh = Function(V)
solve(a==L,uh,DirichletBC(V,Constant(0),"on_boundary"))

# Stiffness matrix without BCs:
A = assemble(a)
# Solution vector:
U = uh.vector()

# Function whose nodal values have interpretation of compliance contributions:
nodal_compliance = Function(V)
as_backend_type(nodal_compliance.vector())\
    .vec().pointwiseMult(as_backend_type(A*U).vec(),as_backend_type(U).vec())

# Check sum of nodal compliance against integral for total energy and sum
# of element contributions to energy:
total_compliance_sum = as_backend_type(nodal_compliance.vector())\
                       .vec().sum()
total_energy_integral = assemble(inner(grad(uh),grad(uh))*dx)
element_energy = project(CellVolume(mesh)*inner(grad(uh),grad(uh)),
                           FunctionSpace(mesh,"DG",0))
total_energy_sum = as_backend_type(element_energy.vector()).vec().sum()
# All give same value:
print(total_compliance_sum)
print(total_energy_integral)
print(total_energy_sum)

# Plot the element contributions to density and the 
from matplotlib import pyplot as plt

# Quite different distributions!
plot(element_energy)
plot(nodal_compliance)

plt.show()

(The use of PETSc linear algebra operations (pointwiseMult and sum) is a bit messy, due to DOLFIN’s layer of abstraction above PETSc; for a DOLFIN vector U, the corresponding PETSc object is accessed via as_backend_type(U).vec().)

9

It’s not clear to me what you mean; perhaps it would be better to create a separate discussion topic about your problem, with more details.

10

that is the topic

11

Hi, thank you for your reply and code. I hadn’t realised that they could be different locally.
What i need is indeed nodal compliance. But why do you interpret AU as compliance? For instance, fo the discrete version of elastic linear system, wouldn’t AU be the force vector? This might be a more theoretical question, but any tip or reference would really help me!

12

The total compliance is \mathbf{U}^T\mathbf{A}\mathbf{U} = \sum_i U_i(\mathbf{A}\mathbf{U})_i, and what I’m considering “nodal compliance” for node i is U_i (\mathbf{A}\mathbf{U})_i (with no summation on the repeated index); that is the interpretation of the pointwiseMult PETSc function.

13

Thank you for your answer! I tried that same analysis you proposed for the 2d linear elastic code i’m working and I didn’t get the same value, though. Also, I noticed that the nodal compliance’s distribution is really weird, concentrated in one region. Could you help me identify what I am doing wrong? I’m copying a snippet of my code here. Thank you!

from dolfin import *
from matplotlib import pyplot as plt
import numpy as np
%matplotlib inline

Nx,Ny,lx,ly = [62,31,2.0,1.0]    
#define mesh
mesh  = RectangleMesh(Point(0.0,0.0),Point(lx,ly),Nx,Ny,'crossed') 
Vvec = VectorFunctionSpace(mesh, 'CG', 1) 
nNodes = (Nx+1)*(Ny+1) + Nx*Ny

class DirBd(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],.0) #left  border    
dirBd = DirBd()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)    
dirBd.mark(boundaries, 1)
ds = Measure("ds")(subdomain_data=boundaries)  
bcd  = DirichletBC(Vvec, (0.0,0.0), boundaries, 1)   #states that a function on Vvec should be igual to (0,0) on the domain defined by boundaries (with value 1)
Load = Point(lx, 0.5*ly)

V = FunctionSpace(mesh, 'CG', 1)
n  = FacetNormal(mesh)
dX = Measure('dx')

E = 1
nu = 0.3
mu,lmbda = Constant(E/(2*(1 + nu))),Constant(E*nu/((1+nu)*(1-2*nu)))

u,v = [TrialFunction(Vvec), TestFunction(Vvec)]
a = 2.0*mu*inner(sym(grad(u)),sym(grad(v)))*dx + lmbda*div(u)*div(v)*dx
K = assemble(a)   
F = assemble( inner(Expression(('0.0', '0.0'),degree=2) ,v) * ds(2))    
U = Function(Vvec)
delta = PointSource(Vvec.sub(1), Load, -1.0)
delta.apply(F) 
bcd.apply(K,F)    
solver = LUSolver(K)
solver.solve(U.vector(), F)  

# Stiffness matrix without BCs:
A = assemble(a)
# Solution vector:
uh = U.vector()

# Function whose nodal values have interpretation of compliance contributions:
nodal_compliance = Function(Vvec)
nodeComp = as_backend_type(nodal_compliance.vector()).vec()
nodeComp.pointwiseMult(as_backend_type(A*uh).vec(),as_backend_type(uh).vec())

# Check sum of nodal compliance against integral for total energy and sum
# of element contributions to energy:
total_compliance_sum = nodeComp.sum()
total_energy_integral = assemble(inner(grad(U),grad(U))*dx)
element_energy = project(CellVolume(mesh)*inner(grad(U),grad(U)),
                           FunctionSpace(mesh,"DG",0))
total_energy_sum = as_backend_type(element_energy.vector()).vec().sum()

print(total_compliance_sum)
print(total_energy_integral)
print(total_energy_sum)

plt.plot(np.array(nodeComp))```
14

I tried that same analysis you proposed for the 2d linear elastic code i’m working and I didn’t get the same value, though.

You need to update the internal energy density from the Poisson equation,

\psi = \frac{1}{2}\vert\nabla u\vert^2

to the elastic energy density,

\psi = \frac{1}{2}\boldsymbol{\sigma}:\boldsymbol{\epsilon}\text{ .}

In your MWE, this would look like the following:

sigma = 2.0*mu*sym(grad(U)) + lmbda*div(U)*Identity(len(U))
psi = 0.5*inner(sigma,grad(U)) # Elastic energy density
total_energy_integral = assemble(2.0*psi*dx)
element_energy = project(CellVolume(mesh)*2.0*psi,FunctionSpace(mesh,"DG",0))

(To understand the factors of 1/2 and 2, refer to my note on this in an earlier answer.)

Also, I noticed that the nodal compliance’s distribution is really weird, concentrated in one region.

If we interpret “nodal compliance” as the work done by an external load on each node, and the external load is only applied to certain nodes, then the nodal compliance will be likewise concentrated. (If that is not what you’re expecting, then perhaps you really are looking for local elastic energy density, not compliance.)

15

Hello @kamensky! After studying some more I’m back to dealing with this problem and I wonder if I could ask you one more thing. What I’m interested is indeed element compliance and I am successfully computing it as you advised me:

#U is the solution in Vvec = VectorFunctionSpace(mesh, 'CG', 1)
S1 = 2.0*mu*inner(sym(grad(U)),sym(grad(U))) + lmbda*div(U)**2 
#with constants for plane stress
compliance = (1/2)*assemble(eps_er*S1*dx(0) + S1*dx(1))

V1 = FunctionSpace(mesh,"DG",0) space for 
eleComp = project(CellVolume(mesh)*S1,V1)

Then I get one constant value for each element of my mesh, which is retangular and crossed

mesh  = RectangleMesh(Point(0.0,0.0),Point(lx,ly),Nx,Ny,'crossed')

But for the rest of my computations (which doesn’t rely on fenics) I need a value associated to each node. So I’m trying to take the mean of the elements around each node, using node coordinates, etc. I wonder however if there is a more straighfoward way of doing this changing the way of projecting it since I remembered your reply:

" You could come up with some nodal proxy for it in the space FunctionSpace(mesh,"CG",1) (which has a degree of freedom at each mesh vertex), but the details would depend on how exactly you want to define the projection."

Thank you!

16

Perhaps a lumped-mass projection (as demonstrated in the function lumped_project here) of the cell-wise energy to FunctionSpace(mesh,"CG",1) would be a good candidate. That would produce a sort of weighted average of the values in cells surrounding each mesh vertex.

1 Like
17

Thank you so much! That works :slight_smile: