How to choose projection space for stress tensor (post-processing)?

variational formulation
1

Hi everyone,
I am a bit confused. I have seen in various examples that the correct post-processing calculation of the normal stress is given by

Vsig = TensorFunctionSpace(mesh, "DG", degree=0) # copied from texts but not understood
sig_num = Function(Vsig, name="Stress Numeric")
sig_num.assign(project(sigma(u, p), Vsig))

In this example I am solving a Stokes flow problem using Taylor Hood elements (P2 - velocity, P1 - pressure), constant viscosity. Why is Vsig chosen as Discontinuous Lagrange of degree zero? Since the stress tensor uses both u and p, I would have thought that it would be CG 1 (like the pressure). Is there a reason it should be piecewise scalar, or is it just a decision I have to make for visualisation purposes?

I think that is the correct space to use for Taylor Hood mixed elements, but I wanted to experiment with different combinations of spaces (something similar to what is done in Chapter 20 of the Fenics book), but I don’t know how to choose the stress tensor space in each case.

This is a general question about the spaces, but I attach a MWE for reference, where I am using the manufactured solution from Chapter 20 and impose velocity Dirichlet on all of the boundary.

from dolfin import *


# Define the boundary domains
class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary


def epsilon(v):
    return sym(as_tensor([[v[0].dx(0), v[0].dx(1)],
                          [v[1].dx(0), v[1].dx(1)]]))


# stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-p*Identity(2)


# Mesh
mesh = UnitSquareMesh(50, 50, "crossed")

x = SpatialCoordinate(mesh)

# Manufactured solution
u_solns_ex = ['sin(4*pi*x[0])*cos(4*pi*x[1])', '-cos(4*pi*x[0])*sin(4*pi*x[1])']
p_solns_ex = ['pi*cos(4*pi*x[0])*cos(4*pi*x[1])', '-9*pi*cos(4*pi*x[0])']
f_solns_ex = ['28*pi*pi*sin(4*pi*x[0])*cos(4*pi*x[1])', '-36*pi*pi*cos(4*pi*x[0])*sin(4*pi*x[1])']

mu = Constant(1.0) # constant viscosity
n = FacetNormal(mesh)

#Taylor-Hood
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = V * Q
W = FunctionSpace(mesh, TH)

# MINI
# P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
# B = FiniteElement("Bubble", mesh.ufl_cell(), mesh.topology().dim() + 1)
# V = VectorElement(NodalEnrichedElement(P1, B))
# Q = P1
# W = FunctionSpace(mesh, V * Q)

# CD
# V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
# Q = FiniteElement("DG", mesh.ufl_cell(), 0)
# CD = V * Q
# W = FunctionSpace(mesh, CD)

# Define trial and test functions
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))

# Manufactured solution of the easy case
u_ex = Expression((u_solns_ex[0], u_solns_ex[1]), element=V, domain=mesh)
f_ex = Expression((f_solns_ex[0], f_solns_ex[1]), element=V, domain=mesh)


bcs = DirichletBC(W.sub(0), u_ex, Boundary())

# Define the variational problem
a1 = (inner(sigma(u, p), epsilon(v))) * dx
a2 = (- div(u) * q) * dx
a3 = (- dot(f_ex, v)) * dx
F = a1 + a2 + a3

solve(F == 0, w, bcs)
(u, p) = w.split()

# Post-processing calculations
# Numerical values for the stress tensor
Vsig = TensorFunctionSpace(mesh, "DG", degree=0) # copied from texts but not understood
sig_num = Function(Vsig, name="Stress Numeric")
sig_num.assign(project(sigma(u, p), Vsig))

I commented out the other options of mixed spaces I am experimenting with. If anyone could shed light on this subject, I would really appreciate it.

2

I’m not sure what you’re asking here since it seems like you’ve answered your own question in your question.

Unless I’ve made a stupid mistake: If \sigma_h = \nabla \vec{u}_h - p_h \mathrm{I}, with your Taylor-Hood element, \vec{u}_h \in [V^{p}_{CG}]^{d}, \nabla \vec{u}_h \in [V^{p-1}_{DG}]^{d\times d}, and p_h \in V^{p-1}_{CG}. So, their linear combination is housed in \sigma_h \in [V^{p-1}_{DG}]^{d\times d}.

Things will get trickier if you incorporate a viscosity function.

3

But, with Taylor Hood, isn’t p = 2? so the linear combination should be housed in

Vsig = TensorFunctionSpace(mesh, "DG", degree=1) #

not degree=0? [I feel like I’m missing something really obvious…]

Also, could you elaborate a bit about the viscosity function? In my real problem, I couple a temperature problem via an exponential viscosity function, which makes everything very tricky. In that case, if the temperature function lives in the same space as the pressure, could I still use the same Vsig as I defined above?

Thanks again!

4

Indeed, if p=2 then one would aptly note p-1=1, and project into the appropriate space.

Should you incorporate a nontrivial viscosity function, your projection operation would yield an approximation computed from your approximate solution. So one should carefully consider whether this is an issue or not in terms of the precision and accuracy desired from the simulation.

1 Like
5

Thank you! It seems like I need to read a lot more about projections to understand the precision of my stress calculations, as I have been getting some weird results and that’s probably the issue.

Thanks for your patience. I asked again about the order because in the 2D elasticity example, they use P2 elements, but the stress is projected onto a piecewise constant (degree zero) space for post-processing.

6

As far as I am aware, that is a workaround from before XDMFFile.write_checkpoint was introduced. Prior to this function you could only visualize either:

  1. CG -1 functions (all other functions would be interpolated into CG 1
  2. DG-0 functions

However, with the introduction of write_checkpoint, you can in theory write arbitrary order CG and DG to file for visualization in Paraview

3 Likes
7

That clarifies things significantly, I really appreciate it.