How can I use grad(h) when h is a Python "float"?

Question

I know this sounds peculiar, but I'm wondering if there's a way to do this. I'm solving a problem in which the source term depends on the gradient of a function defined over the entire mesh. This function, h, is computed as a Python expression. Most often this expression uses a FEniCS object defined over the mesh as an argument, and so grad knows the domain that it's working with. However, I'd like to have the flexibility to allow for this function, h, to be a constant. I can include a conditional statement to achieve the desired behavior, but I was curious if there is a more elegant approach.

To make my question more concrete, here is a simple example:

from dolfin import *
from mshr import *

N = 32
mesh = generate_mesh( Rectangle( dolfin.Point( 0., 0. ), 
                                 dolfin.Point( 1., 1. ) ), N )

g = Expression("x[0] + x[1]*x[1]", domain=mesh)
h = 10.
f = dot( grad(g), grad(h) )

This results in an error with grad, and returns

ufl.log.UFLException: Cannot get geometric dimension from an expression with no domains!

On the other hand, if my function is h = 10.*g, everything works as expected.

Answer 1

You could simply multiply with a Constant with the appropriate domain before applying grad. See example below.

from dolfin import *

mesh = UnitSquareMesh(16, 16)
V = FunctionSpace(mesh, "CG", 1)

def F(x, y, u): return x * y**2
def G(x, y, u): return u + sin(x)*y + 2
def H(x, y, u): return 1.0

x, y = MeshCoordinates(mesh)
u = interpolate(Expression("x[0]"), V)

for py_function in [F, G, H]:
    f = Constant(1.0, cell = mesh.ufl_cell()) * py_function(x, y, u)
    print assemble(grad(f)**2 * dx(domain = mesh))

Comments

Thank you. That's just what I was looking for.

Answer 2

You could simply write

h = Expression("0.", domain=mesh)

Comments

Thanks. However, what I'd like is to be able to use a user-specified Python function for h. The reason for this is that the form of h changes for different problems. So I'm looking for an elegant way to handle both situations such as h = 10.*g as well as h = 10..

Answer 3

I think that's impossible to do what you are asking without defining a domain for h (i.e. without to evaluate h on the entire domain), and if h is a type float or int isn't possible to do that (because a float or int object has no the attribute dx).

Comments

Thanks. I think you're probably right that there's no way to do exactly what I'd like. I was hoping that someone could come up with a clever way to convert any h to a UFL form.

---

You can use Sympy, i.e.

from sympy import Symbol, diff
sx = Symbol('sx'); sy = Symbol('sy')
h = sx+sy**2
hx = str(diff(h,sx)).replace("sx","x[0]").replace("sy","x[1]")
hy = str(diff(h,sy)).replace("sx","x[0]").replace("sy","x[1]")
gradh = Expression((hx,hy)) #use in place of grad(h)

---

Sympy has functionality for generating C code:

from sympy import var
from sympy.utilities.codegen import ccode

x, y = var("x[0] x[1]")
f_sympy = x**3 * y**2
f_dolfin = Expression(ccode(f_sympy))

I think UFL should be used instead of sympy if possible though.