Variational problem with user-defined function depending on u(x)

Question

Dear FEniCS community,

as far as I understand, one may include user-defined functions in x in the variational form via Expression objects, but regarding functions in u(x), in UFL only predefined functions (as sin, sqrt, conditional ...) are supported.
Is it possible in FEniCS to solve variational problems that include user-defined functions in trial functions? In particular, I have a function representing an imported distribution defined as

# fct simulating the distribution of an image
def g(x):
    #calculate y,z depending on x
    ... 
    return image.getpixel((math.floor(y),math.floor(z)))

Based on this I want to solve a variational problem like this

u = TrialFunction(V)
v = TestFunction(V)
F = ... + g(nabla_grad(u))*v*dx

F = action(F, u_)
J  = derivative(F, u_, u)   # Gateaux derivative in dir. of u

problem = NonlinearVariationalProblem(F, u_, None, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

Using this python-definition of g, the code crashes with

  Couldn't map 'v_1' to a float, returning ufl object without evaluation.
Traceback (most recent call last):
  File "test/test user defined.py", line 58, in <module>
    F = g(nabla_grad(u))*v*dx
  File "test/test user defined.py", line 38, in g
    pixel = image.getpixel((math.floor(y),math.floor(z)))
  File "/usr/lib/python2.6/site-packages/ufl/expr.py", line 133, in __float__
    return self(()) # No known x
  File "/usr/lib/python2.6/site-packages/ufl/exproperators.py", line 286, in _call
    return _eval(self, arg, mapping, component)
  File "/usr/lib/python2.6/site-packages/ufl/exproperators.py", line 278, in _eval
    return f.evaluate(coord, mapping, component, index_values)
  File "/usr/lib/python2.6/site-packages/ufl/algebra.py", line 215, in evaluate
    tmp *= o.evaluate(x, mapping, (), index_values)
  File "/usr/lib/python2.6/site-packages/ufl/algebra.py", line 303, in evaluate
    return float(a) / float(b)
TypeError: __float__ returned non-float (type Sum)

Answer 1

Your code can't possibly work, because you have defined g as a Python function and are calling it with the argument nabla_grad(u), which is a Python object of type ufl.differentiation.Derivative which I don't think your function g is supposed to handle...

You should try to define g by subclassing Expression, i.e. change your code to

class MyExpression(Expression):
    def eval(self, values, x):
        #calculate y,z depending on x
        ... 
        values[0] = image.getpixel((math.floor(y),math.floor(z)))

g = MyExpression() 

Comments

If I substitute the python-definition of g by your suggested expression object, how do I tell him to use nabla_grad(u(x))? eval only knows x (and is called implicitly).

---

Of course, you're right. Hm.. well I'm afraid I don't know either -.-'

---

Ah yes: I guess, the way to go is to first define u as a Function(V) (not TrialFunction), so u = u_ in your code. Then, interpolate nabla_grad(u) onto some new Finite Element Function, i.e.

gradu = Function(V)
gradu.interpolate(nabla_grad(u))

With this, you can use the following in your Expression:

class MyExpression(Expression):
    def eval(self, values, x):
        x = gradu(x)
        #calculate y,z depending on x
        ... 
        values[0] = image.getpixel((math.floor(y),math.floor(z)))

---

Mmmh ... Thanks, but I guess this does not solve the problem. I tried to make a little example code solving the Poisson equation with constant rhs and dirichlet zero boundary (using a nonlinear solver). to test your suggestion

from dolfin import *

mesh = UnitSquareMesh(10,10,'crossed')

#function spaces
V = FunctionSpace(mesh, 'CG', 2)  
W = VectorFunctionSpace(mesh, 'CG', 1)

u  = TrialFunction(V)
v = TestFunction(V)
u_ = project(Expression('x[0]*x[0]',degree=2), V)

gradu = project(nabla_grad(u_),W)

class MyExpression(Expression):
    def eval(self, v, x):
      v = gradu(x)
    def value_shape(self):
        return (2,)

g = MyExpression()

#define geometry
h = CellSize(mesh)

#variational formulation for the poisson equation
F = (inner(nabla_grad(u_),nabla_grad(v))-Constant(1.0)*v)*dx 
#F = (inner(g,nabla_grad(v))-Constant(1.0)*v)*dx # using the user-defined expression g

#boundary conditions
def u0_boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, Constant(0.0), u0_boundary)

J  = derivative(F, u_, u)   # Gateaux derivative in dir. of u

problem = NonlinearVariationalProblem(F, u_, bc, J)
solver  = NonlinearVariationalSolver(problem)

solver.solve()

plot(u_)
interactive()

Using
F = (inner(g,nabla_grad(v))-Constant(1.0)*v)*dx
instead of
F = (inner(nabla_grad(u_),nabla_grad(v))-Constant(1.0)*v)*dx
the code crashes('ufl.log.UFLException: Cannot replace unknown element domain without unique common domain in form.`). If you interpolate the gradient of u_ to a function gradu, the nonlinear solver regards gradu as constant (and independent of u_). Hence, the problem becomes unsolvable.

---

Yes, of course now you have to implement the nonlinear solving yourself, because you have to interpolate gradu in every step of the Newton iteration!

---

I had hoped to avoid to calculate the Jacobian myself. But I guess there is no easy workaround. Thank you!

Answer 2

I'm sorry I couldn't provide you with a quick and easy solution, but there probably isn't one, for very general reasons:

In the Jacobian of your nonlinear form, the derivative of your (original) function g has to enter somewhere. But how is fenics supposed to know it? The automatic derivative(F,u_,u) tries to apply symbolic differentiation to a symbolically defined ufl form. Hence it can only handle the limited number of nonlinear expressions that are predefined in ufl, as you already noticed in the beginning.

Hence, to be able to use a Newton method, you have to supply the derivative of g yourself (if this is even possible), build the Jacobian bilinear form yourself and implement the Newton iteration yourself, since the derivative of g probably is again a complicated Expression that you have to interpolate manually to use it inside a ufl form.

Maybe Newton isn't the way to go? If you're lucky, you maybe can get away with a simple fixed-point iteration (Picard iteration). This should be implementable rather straightforward in fenics, but still, you will have to write the iteration loop yourself...

Some things just can't be done automatically ;) I hope you figure it out! Cheers

Answer 3

Maybe this discussion is helpful? http://fenicsproject.org/qa/8169/expressions-automated-differentiation-custom-functions

Answer 4

You can use the "R" space + defining functions as expressions or via UFL.
See the following code:

from fenics import *

mesh = UnitIntervalMesh(100)
dim = 5
R = VectorFunctionSpace(mesh, "R", 0, dim=dim)
x = SpatialCoordinate(mesh)

Define reduced basis

rb = [sin(pix[0](i+1)) for i in range(dim)]

udofs = TrialFunction(R)
vdofs = TestFunction(R)
u = sum([udofs[i]rb[i] for i in range(dim)])
v = sum([vdofs[i]rb[i] for i in range(dim)])

L = v*dx(degree=3)

a = uvdx + inner(grad(u), grad(v))*dx(degree=3)

xdofs = Function(R)
solve(a == L, xdofs)

Get full solution

x = sum([xdofs[i]*rb[i] for i in range(dim)])
plot(x, interactive=True)