Average nodal stress and principal stress (convert tensor function from DG1 to CG2)

1

How can I compute the averaged nodal stress tensor components and principal stresses for an elasticity problem?

In other words, I have a DG(n-1) tensor function of the stress obtained by projection from CG(n) displacements
and I want to “convert” this to a CG(n) tensor function by averaging the stress tensor contribution from each element connected to each node (same as in Abaqus). Using these average nodal stresses I want to compute the principal stress values (eigenvalues) and directions (eigenvectors) at the nodes (one value for each direction component at each node).

Here is a simple example of what I have so far based on the elasticity and hyperelasticity demos:

from dolfin import *

# Degree of finite element for displacements (P1, P2 etc.)
degree = 2

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

# Create mesh and define function space
mesh = UnitCubeMesh(6, 6, 6)
V = VectorFunctionSpace(mesh, "CG", degree)

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

# Define Dirichlet boundary (x = 0 or x = 1)
c = Constant((0.0, 0.0, 0.0))
r = Expression(("scale*0.0",
                "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
                "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
                scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3, degree=2)

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]

# Elasticity parameters
E = 1.0e9
nu = 0.3
mu = E/(2.0*(1.0 + nu))
lmbda = E*nu/((1.0 + nu)*(1.0 - 2.0*nu))

# Stress as a function of displacement
def stress(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(len(v))

# Define variational problem and solve
u = TrialFunction(V)
v = TestFunction(V)
U = Function(V)
a = inner(stress(u), grad(v))*dx
f = Constant((0.0, 0.0, 0.0))
L = inner(f, v)*dx
solve(a == L, U, bcs)

# Discontinuous element nodal stress
VDGtensor = TensorFunctionSpace(mesh, "DG", degree - 1)
S = project(stress(U), V=VDGtensor, solver_type='cg', preconditioner_type='jacobi')

# TODO: compute average nodal stress

# TODO: compute principal stresses and directions using average nodal stress

# Rename output variables
U.rename("U", "U")
S.rename("S", "S")

# Save output in XDMF format
with XDMFFile("displacement.xdmf") as f:
    f.parameters["flush_output"] = True
    f.parameters["functions_share_mesh"] = True
    f.write(U, 0.0)
    f.write(S, 0.0)

EDIT: I guess projecting onto a continuous space of the same degree as the displacements would give me something similar to the averaged nodal stress:

VCGtensor = TensorFunctionSpace(mesh, "CG", degree)
S = project(stress(U), V=VCGtensor, solver_type='cg', preconditioner_type='jacobi')

But there are two problems with this:

  1. I need to compare the results with Abaqus which uses the averaged nodal stress.
  2. I cannot use global projection because my model is too big and I run out of memory. To obtain the discontinuous stress I use the local projection method described here which gives the same results as the simplified example above.