How to store diferential operator in a Function object

1

Hello everyone,

I have a 1D function f which is in fact the stream function of a vector field \vec{v}.

f is the solution of a PDE and is stored in a Function() object. I would like to get the vector field \vec{v} and store it in a Function() object as well in order to use it in another advection PDE like \vec{v} \cdot \nabla q = 0 and be able to evaluate it in different points of domain, anyway all my attempts to do so have failed.

The main idea would be:

import ufl
from dolfin import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(20,20)

V1 = FunctionSpace(mesh, 'DG', 2)
Vdim = VectorFunctionSpace(mesh, 'DG', 2)


#Example streamfunction
f = Expression('1/(x[0]*x[0] + x[1]*x[1] + 0.01)', element = V1.ufl_element())
f_interp = interpolate(f, V1)

#Get the vector field from the stream function:
v = as_vector((-f_interp.dx(1), f_interp.dx(0)))

The only other option I came up with would be sampling f in the domain to compute the derivatives with numpy and then interpolating the values back to fenics but this seems like a long detour.

Thanks in advance.

2

You can project the differential operator into a suitable function space if you use dolfin.

If you want to use interpolation for such an operation, consider DOLFINx, and the example from:

It has also been explained for the curl operator in:
https://jorgensd.github.io/dolfinx-tutorial/chapter3/em.html

1 Like
3

Hi @dokken, thanks for your reply.

Projection should suffice but, how can one create a valid expression to project onto the function space?

I tried:

class Expr(UserExpression):

    def eval(self, value,x):
        value[0] = -f.dx(1)
        value[1] = f.dx(0)

    def value_shape(self):
                return (2,)

project(Expr(), Vdim)

However I get:

ValueError: setting an array element with a sequence.

4 
from dolfin import *

mesh = UnitSquareMesh(1, 1)
f = Expression("x[0]*x[1]", degree=1, domain=mesh)
grad_f = project(grad(f), VectorFunctionSpace(mesh, "DG", 0))
print(grad_f.vector()[:])

gives:

[ 0.  1.  1.  0.]

Make sure you use appropriate spaces.

1 Like
5

I guess in his case it would be:

grad_f = project(as_vector((-u.dx(1), u.dx(0))), VectorFunctionSpace(mesh, "DG", 0))
2 Likes
6

Indeed, I had not thought of that.

Thank you both.