Compute element-wise inverse of a function

Question

Hi,

I would like to compute the element-wise inverse of a Function (in order to compute relative error between two fields) :

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "Lagrange", 1)
u1 = Function(V)
u2 = Function(V)
uerr = project(u1-u2,V)     #error
uerr_rel = project(uerr/u1,V)   #relative error : does not work...

Many thanks in advance !

Answer 1

Hi, what exactly is the error that you are getting? The obvious problem is that in your example your are dividing by zero. If you u1 is nonzero then the code works, e.g.

from dolfin import * 

parameters['linear_algebra_backend'] = 'PETSc'

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'CG', 1) 

u1 = interpolate(Expression('1+x[0]+x[1]'), V)
u2 = interpolate(Expression('1+x[0]'), V)
uerr = project(u1 - u2,V)
uerr_rel = project(uerr/u1, V)

print uerr_rel.vector().array() 

As an alternative to projection you could consider doing the operation pointwise.

Uerr = Function(V)
Uerr.assign(FunctionAXPY([(1., u1), (-1., u2)]))

Uerr_rel = Function(V)
# Make Uerr_rel have values of Uerr
Uerr_rel.vector().axpy(1, Uerr.vector())

# Now divide the vectors
as_backend_type(Uerr_rel.vector()).vec().pointwiseDivide(
        as_backend_type(Uerr_rel.vector()).vec(),
        as_backend_type(u1.vector()).vec())

print Uerr_rel.vector().array() 

You might also want to update the ghost values as done here.