Norm Calculation for Tensor / Matrix (VectorFunctionSpace).

Question

Hi,

I am working on implementation of Turbulence Models. Most of my code is written in a library (private) and I am using fenics to verify most of my code. (For eg. Element Matrix etc.)

The turbulence production term as per the, $$k - \epsilon $$ Model is computed as

$$ \nu_t \frac{ | \nabla u + \nabla u^T | }{2} $$

where,
$u = \text{Velocity} $ and $\nabla u^T$ represents transpose of the tensor

I am computing the term ( i.e. $ | \nabla u + \nabla u^T | $) via Forbenius Norm

I am stuck with the implementation for fenics. I am not sure as how to compute the norm for
the tensor term. So far,

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from   dolfin import *
import numpy  as     np

co = np.asarray([[0,0],[1,0],[0,1]], dtype=np.float64)

mesh = Mesh()
editor = MeshEditor()

# Create Reference Element
editor.open(mesh, 2, 2)
editor.init_vertices(3)
editor.init_cells(1)
for i in range(len(co)):
    editor.add_vertex(i, co[i, :])
editor.add_cell(0, np.asarray([0,1,2], dtype=np.uintp))
editor.close()

# Plot to check
# plot(mesh, interactive=True)

V = VectorFunctionSpace(mesh, "CG", 2)

# Velocity
U = Function(V)

# Production
P = grad(U) + grad(U).T

# HOW TO COMPUTE THE NORM FOR THE MATRIX

Can anyone please help me with this.

Answer 1

P will have a different value in each vertex. You can project it into an approproate function space like this:

Q = FunctionSpace(mesh, "DG", 1) # derivative of CG2 yields a DG1 function
Pnorm = project(sqrt(tr(P.T*P)), Q)

for x in co:
    print x, Pnorm(x)

If you instead use CG1, you will get one value per cell for the frobenius norm. You could also of course project directly do a DG0-function if you wish.