### Trace of a tensor field in Hdiv x Hdiv

164

views

0

Hello,

I'm interested in approximating the solution of the so-called Pseudostress-velocity formulation of the Stokes problem. Denoting by "sigma" and "u" the Pseudostress (tensor field) and the velocity (vector field), respectively, my code is given by:

Thanks!

I'm interested in approximating the solution of the so-called Pseudostress-velocity formulation of the Stokes problem. Denoting by "sigma" and "u" the Pseudostress (tensor field) and the velocity (vector field), respectively, my code is given by:

```
from dolfin import *
from mshr import *
from math import log
# Exact solution
sigma_exact = Expression((("1/3 - mu*pow(pi,2)*cos(pi*x[0])*cos(pi*x[1]) - pow(x[0],2)", \
"mu*pow(pi,2)*sin(pi*x[0])*sin(pi*x[1])"), \
("-mu*pow(pi,2)*sin(pi*x[0])*sin(pi*x[1])", \
"mu*pow(pi,2)*cos(pi*x[0])*cos(pi*x[1]) - pow(x[0],2) + 1/3")), degree = 8,mu = 1)
u_exact = Expression(("-pi*cos(pi*x[1])*sin(pi*x[0])","pi*cos(pi*x[0])*sin(pi*x[1])"), degree = 8)
p_exact = Expression("pow(x[0],2)-1/3", degree = 2)
z_exact = Constant(0)
# Define source term
f = Expression(("2*x[0] - 2*mu*pow(pi,3)*cos(pi*x[1])*sin(pi*x[0])", \
"2*mu*pow(pi,3)*cos(pi*x[0])*sin(pi*x[1])"), degree = 8, mu = 1)
# lists to hold errors, meshsize, and number of elements
err_sigma = []
err_u = []
err_p = []
err_z = []
hmax = []
N = [4, 8, 16, 32]
for nn in N:
# Create mesh and define function space
mesh = RectangleMesh(Point(0.,0.), Point(1.,1.), nn, nn, "right/left")
# Function spaces
order = 1
H = VectorElement("RT",mesh.ufl_cell(), order) # Hdiv
Q = VectorElement("DG",mesh.ufl_cell(), order-1) # L2
Z = FiniteElement("R",mesh.ufl_cell(), 0) # Lagrange multiplier
W = FunctionSpace(mesh,MixedElement([H, Q, Z]))
# Define trial and test functions
(sigma, u, z) = TrialFunctions(W)
(tau, v, r) = TestFunctions(W)
# Define deviatoric tensor
def deviatoric(sigma):
return sigma-(1/2)*tr(sigma)*Identity(2)
# Outward unit normal vector to the boundary
n = FacetNormal(mesh)
# Define variational form
a = inner(deviatoric(sigma), deviatoric(tau))*dx \
+ inner(u,div(tau))*dx + inner(v,div(sigma))*dx \
+ z*tr(tau)*dx +r*tr(sigma)*dx
F = -inner(f,v)*dx + inner(tau*n,u_exact)*ds
# Function to solve for
X = Function(W)
solve(a == F, X)
(sigma, u, z) = X.split()
# Postprocessed pressure
pressure = -(1/2)*tr(sigma)
# Mesh-size
h = mesh.hmax()
hmax.append(h)
# Compute the errors
u_error = u - u_exact
z_error = z - z_exact
p_error = p_exact - pressure
err_sigma.append(errornorm(sigma_exact, sigma, norm_type = 'Hdiv'))
err_u.append(errornorm(u_exact, u, norm_type = 'L2'))
err_p.append(sqrt(abs(assemble(p_error*p_error*dx))))
err_z.append(sqrt(abs(assemble(z_error*z_error*dx))))
# Rate of convergence
i = 1
comm = mpi_comm_world()
print("h, ,err_S, rat_S, err_u, rat_u, err_c, rat_c")
while i < len(err_u) :
rate_p = log(err_p[i-1]/err_p[i])/log(hmax[i-1]/hmax[i])
rate_u = log(err_u[i-1]/err_u[i])/log(hmax[i-1]/hmax[i])
rate_sigma = log(err_sigma[i-1]/err_sigma[i])/log(hmax[i-1]/hmax[i])
rate_z = log(err_z[i-1]/err_z[i])/log(hmax[i-1]/hmax[i])
if MPI.rank(comm) == 0:
print("%2.4f %2.2e %2.2f %2.2e %2.2f %2.2e %2.2f" % \
(hmax[i], err_sigma[i], rate_sigma,err_u[i], rate_u, err_z[i], rate_z))
i = i+1
```

The pressure can be computed by a postprocessing formula in terms of the trace of the Pseudoestress tensor "sigma"; however, I think this does not work quite well in FeniCS because the L2-error of the pressure gives the following message:```
Solving linear variational problem.
This integral is missing an integration domain.
Traceback (most recent call last):
File "Pseudostress_based_Formulation.py", line 96, in <module>
err_p.append(sqrt(abs(assemble(p_error*p_error*dx))))
File "/usr/local/lib/python2.7/dist-packages/ufl/measure.py", line 446, in __rmul__
error("This integral is missing an integration domain.")
File "/usr/local/lib/python2.7/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: This integral is missing an integration domain.
Aborted (core dumped)
```

Recall that the degrees of freedom for the Raviart-Thomas space of lowest order is given by the normal components of sigma on the edges of the mesh, so what does tr(sigma) mean? Thanks!

Community: FEniCS Project

Please login to add an answer/comment or follow this question.