Element-wise norm of a Function defined over a VectorFunctionSpace

Question

Is there a way to compute the norm of a vector Function without doing something like

V     = VectorFunctionSpace(mesh)
U     = Function(V)
u,v,w = U.split()
norm  = sqrt(u**2 + v**2 + w**2)

??

Answer 1

Hi, I don't think you need to split.

from dolfin import *                                                             

mesh = Mesh(Rectangle(-10, -10, 10, 10) - Circle(0, 0, 0.1), 100)                
V = FunctionSpace(mesh, "CG", 1)                                                 
W = VectorFunctionSpace(mesh, "CG", 1)                                           

u = interpolate(Expression(("x[0]/sqrt(pow(x[0], 2) + pow(x[1], 2))",            
                            "x[1]/sqrt(pow(x[0], 2) + pow(x[1], 2))")), W)       

u_norm = project(sqrt(inner(u, u)), V)      

Comments

This method does not appear to work as well when taking the gradient of a function. For example,

from dolfin import *
import numpy as np

n    = 3
mesh = BoxMesh(-1, -1, 0, 1, 1, 2, n, n, n)

# define function space :
Q = FunctionSpace(mesh, "CG", 2)

# function to be evaluated :
f = Expression('sqrt(pow(x[0],2) + pow(x[1], 2) + pow(x[2], 2))')
f = interpolate(f, Q)

# formulate variables :
gradf   = grad(f)
u,v,w   = gradf

u       = project(u)
v       = project(v)
w       = project(w)

u_v     = u.vector().array()
v_v     = v.vector().array()
w_v     = w.vector().array()

norm_U  = np.sqrt(u_v**2 + v_v**2 + w_v**2)
norm_Uf = project(sqrt(inner(gradf, gradf)), Q)

n1 = np.sqrt(sum(norm_U**2))
n2 = np.sqrt(sum(norm_Uf.vector().array()**2))
print n1, '?=', n2

Here, n1 should be equal to n2, but instead it prints:

8.49310912965 ?= 19.0029436322

which makes it seem that the 'project' method greatly overestimates the norm. Can you help?

---

Hi, sorry for late response. With

norm_Uf = project(sqrt(inner(gradf, gradf)))

the results are much closer

8.49310912965 ?= 8.01289207971

I consider this a more fair comparison as in this way, all functions are in CG1. Before
you had components in CG1 and the norm in CG2. Still, I don't see why the numbers
should be the same.

---

Yes, this makes sense with CG1 vs. CG2, thanks.

I would like the two numbers to be the same so I can create a unit vector - I can create unit vectors with norms equal to one with the 'array' method, but the 'project-inner' method does not. This makes me wonder why I would want to use the non-array method if it introduces errors.

In essence, I know that the numbers do not have to be the same, but for accuracy, I'd like them to be.