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 »