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

164
views
0
5 months ago by
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:
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.

Similar posts:
Search »