Map ufl expression?

variational formulation
1

I can appy functions to ufl expression like I apply mynorm here:

from fenics import *
mesh = UnitIntervalMesh(10)
pt = SpatialCoordinate(mesh)
E = as_vector((pt[0],))
def mynorm(E):
    return sqrt(E**2)

assemble(mynorm(E)*dx)

This works, because mynorm is build out of primitives supported by fenics. How to do it, if I replace mynorm by an arbitrary function double -> double? Like

def mynorm(E):
    import math
    return math.sqrt(E**2)

In my real use case I want to apply a function written in Julia that knows nothing about fenics.

2

I guess you can write a user-expression. See for instance:
https://bitbucket.org/fenics-project/dolfin/src/946dbd3e268dc20c64778eb5b734941ca5c343e5/python/demo/documented/mixed-poisson/demo_mixed-poisson.py.rst?at=master#demo_mixed-poisson.py.rst-205

1 Like
4

I got it somewhat working. I can evaluate, but get a warning and it throws an error if I try to assemble:

from fenics import *

class MyNorm(UserExpression):
    def __init__(self, E, **kwargs):
        super(MyNorm, self).__init__(**kwargs)
        self.E = E
        
    def eval(self, out, pt):
        import math
        out[0] = math.sqrt((self.E**2)(pt))
    
    def value_shape(self):
        return ()
    
mesh = UnitIntervalMesh(10)
pt = SpatialCoordinate(mesh)
E = as_vector((pt[0],))

MyNorm(E)([0.2]) # FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
# if arg in ("+", "-"):

assemble(MyNorm(E)*dx) # Error: This integral is missing an integration domain.
5

Try assemble(MyNorm(E)*dx(domain=mesh))

1 Like
6

Thanks, that works for me. Do you know how to fix the warning? Also can MyNorm remember the domain, just like the fenics.sqrt(E**2) does?

7

I do not know how to remove the warning.
To get the expression to take a domain, you need to supply what function space the expression should be interpolated to:

from fenics import *

class MyNorm(UserExpression):
    def __init__(self, E, **kwargs):
        super(MyNorm, self).__init__(**kwargs)
        self.E = E

    def eval(self, out, pt):
        import math
        from IPython import embed;embed()
        out[0] = math.sqrt((self.E**2)(pt))

    def value_shape(self):
        return ()
mesh = UnitIntervalMesh(10)
pt = SpatialCoordinate(mesh)
E = as_vector((pt[0],))

element = FiniteElement("CG", mesh.ufl_cell(),  1)
MyNorm(E, domain=mesh, element=element)([0.2]) # FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
# if arg in ("+", "-"):

assemble(MyNorm(E, domain=mesh, element=element)*dx)
1 Like
8

Thanks! If I to assemble(sqrt(E**2)*dx) it just works without explicit domain or element I guess the domain is passed forward from E, which got it passed from

pt = SpatialCoordinate(mesh)
E = as_vector((pt[0],))

How about the element? I never specified any element in the orignal code. Is FiniteElement("CG", domain.ufl_cell(), 1) a default element that gets choosen by fenics for convenience?

9

The default is a CG 2 space. However, this should be chosen as matching the expressions nature. You can also select DG spaces or other spaces.

1 Like
10

It seems that derivatives of UserExpressions silently vanish:

from fenics import *

class MapScalar(UserExpression):
    def __init__(self, f, base_expr, **kwargs):
        self._base_expr = base_expr
        self._f = f
        super(MapScalar, self).__init__(**kwargs)
        
    def eval(self, out, pt):
        out[0] = self._f(self._base_expr(pt))
        
    def value_shape(self):
        return ()
    
mesh = UnitIntervalMesh(10)
pt = SpatialCoordinate(mesh)

import math
assert near(
    MapScalar(math.sqrt, pt[0])([4]),
    2)

assert near(
    assemble(MapScalar(math.sqrt, pt[0])**2*dx(domain=mesh)),
    0.5, 
    eps=0.01)


V = FunctionSpace(mesh, "P", 1)
u = Function(V)
v = TestFunction(V)

F = inner(MapScalar(math.sqrt, u**2), v) * dx
u.assign(Constant(0.1))
assemble(derivative(F, u)).array() # all zero

How to fix this? Is there a way to tell fenics to use finite differences for the UserExpression?

11

Implement the derivative yourself. There is no way for the user expression to know how to differentiate for u.

12

Where to implement it? Can I once and for all implement e.g. derivative(MapScalar(math.sqrt, x), x) and fenics will take care of the rest, or do I need to manually build the derivative of each ufl that involves MapScalar(math.sqrt, x) ?

13

This is slightly non-trival, as it requires usage of ufl functions not exposed directly in dolfin.
To make it simpler for you, I made a simpler example to illustrate how to do this (and verify that the method is correct). At the end of the example I am using your original form.

from fenics import *
import ufl


class MapScalar(UserExpression):
    def __init__(self, f, base_expr, **kwargs):
        self._base_expr = base_expr
        self._f = f
        self.derivatives = {}
        for operand in self._base_expr.ufl_operands:
            if isinstance(operand, ufl.Coefficient):
                self.derivatives[operand] = ufl.diff(self._f(self._base_expr), operand)
        super(MapScalar, self).__init__(**kwargs)

    def eval(self, out, pt):
        out[0] = self._f(self._base_expr(pt))

    def value_shape(self):
        return ()


mesh = UnitIntervalMesh(2)
pt = SpatialCoordinate(mesh)

V = FunctionSpace(mesh, "CG", 1)
u = Function(V, name="u")
v = TestFunction(V)
element = FiniteElement("CG", mesh.ufl_cell(), 2)
expression = MapScalar(sqrt, u**2, domain=mesh, element=element)

# Simple example with "F" being a functional
F = expression*dx
u.assign(project(pt[0]**2, V))
direct_application = assemble(inner(expression.derivatives[u],v)*dx).get_local()
print(direct_application)
analytical =  assemble(0.5*inner(2*u*(u**2)**(-0.5),v)*dx).get_local()
print(analytical)


derivative_dict = {expression:expression.derivatives[u]}
dFdu = ufl.algorithms.apply_derivatives.apply_derivatives(derivative(F, u, coefficient_derivatives=derivative_dict))
using_derivative = assemble(dFdu).get_local()
print("derivative:", using_derivative )


# Matrix example
F = inner(expression, v) * dx
dFdu = ufl.algorithms.apply_derivatives.apply_derivatives(derivative(F, u, coefficient_derivatives=derivative_dict))
print(assemble(dFdu).array())
2 Likes
14

Thanks a lot! So where would I pass the derivative of my function?
Would I do something like

self.derivatives[operand] = derivative_of_mysqrt(ufl.diff(self._base_expr, operand))
15

I am not sure what you are asking for. Please have a careful look at the example as everything to make you MWE work is supplied there.