ufl team mailing list archive

[Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx>]

On 4 August 2011 14:56, Chaffra <question166940@xxxxxxxxxxxxxxxxxxxxx> wrote:
> New question #166940 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/166940
>
> Dear all,
>
> I want to be able to solve poisson's equation with some abs(absolute value) functions in it, but ffc breaks with
>
> Exception: False value of Condtional should only be one function: {((1, 'k', 3, 3),): 'FE0[ip][k]*C[1]'}

I tried the following simplified version of your code as a pure UFL
file and compiled it with FFC which results in a different, but
related crash (probably because I'm on the dev version of FEniCS):

V  = FiniteElement("Lagrange", triangle, 1)
u  = Coefficient(V)
du = TrialFunction(V)
v  = TestFunction(V)
L1 = abs(u)*v*dx
J = derivative(L1,u,du)

Found Argument in Conditional(EQ(Coefficient(FiniteElement('Lagrange',
Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
(), {}), Product(Argument(FiniteElement('Lagrange', Cell('triangle',
Space(2)), 1, None), 1),
Conditional(EQ(Coefficient(FiniteElement('Lagrange', Cell('triangle',
Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
Conditional(LT(Coefficient(FiniteElement('Lagrange', Cell('triangle',
Space(2)), 1, None), 0), Zero((), (), {})), IntValue(-1, (), (), {}),
IntValue(1, (), (), {}))))), this is an invalid expression.
Traceback (most recent call last):
  File "/home/oelgaard/software/fenics/bin/ffc", line 195, in <module>
    sys.exit(main(sys.argv[1:]))
  File "/home/oelgaard/software/fenics/bin/ffc", line 176, in main
    compile_form(ufd.forms, ufd.object_names, prefix, parameters)
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ffc/compiler.py",
line 150, in compile_form
    analysis = analyze_forms(forms, object_names, parameters, common_cell)
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ffc/analysis.py",
line 64, in analyze_forms
    common_cell) for form in forms)
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ffc/analysis.py",
line 64, in <genexpr>
    common_cell) for form in forms)
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ffc/analysis.py",
line 151, in _analyze_form
    ffc_assert(len(compute_form_arities(preprocessed_form)) == 1,
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/formtransformations.py",
line 332, in compute_form_arities
    parts = compute_form_with_arity(form, arity)
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/formtransformations.py",
line 317, in compute_form_with_arity
    res = transform_integrands(form, _transform)
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/transformations.py",
line 799, in transform_integrands
    integrand = transform(itg.integrand())
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/formtransformations.py",
line 313, in _transform
    e, provides = pe.visit(e)
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/transformations.py",
line 144, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/transformations.py",
line 148, in visit
    r = h(o)
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/formtransformations.py",
line 68, in expr
    error("Found Argument in %s, this is an invalid expression." % repr(x))
  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/log.py",
line 148, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Found Argument in
Conditional(EQ(Coefficient(FiniteElement('Lagrange', Cell('triangle',
Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
Product(Argument(FiniteElement('Lagrange', Cell('triangle', Space(2)),
1, None), 1), Conditional(EQ(Coefficient(FiniteElement('Lagrange',
Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
(), {}), Conditional(LT(Coefficient(FiniteElement('Lagrange',
Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})),
IntValue(-1, (), (), {}), IntValue(1, (), (), {}))))), this is an
invalid expression.

The problem is the UFL derivative of abs(u) which can be inspected by:

from ufl.algorithms import expand_derivatives
V  = FiniteElement("Lagrange", triangle, 1)
u  = Coefficient(V)
du = TrialFunction(V)
v  = TestFunction(V)
L1 = abs(u)*v*dx
J1 = derivative(L1,u,du)
print expand_derivatives(J1)

{ v_{-2} * (((w_0) == (0)) ? (0) : (v_{-1} * (((w_0) == (0)) ? (0) :
(((w_0) < (0)) ? (-1) : (1))))) } * dx0

Here, 'v_{-2}' is 'v', 'w_0' is 'u' and 'v_{-1}' is 'du'.

This is identical to:

c = conditional( eq(u, 0), 0, du*conditional(eq(u, 0), 0,
conditional(lt(u,0),-1,1) ) )
J = v*c*dx
print J

Notice that 'du' is present inside the first conditional (the second
return value) which triggers the error.

I'm not sure that it makes sense to allow 'du' in the condition which
is being tested e.g.:

c = conditional( eq(du, 0), 0, 1 )

but there should be no problems with allowing it in the return values like:

c = conditional( eq(u, 0), du, 2*du )

However, with the particular form in question, the return value of 'c'
will be either 0, du or -du, so J in the case of u == 0 will be a
linear form!!!

Things would work fine if UFL could move 'du' outside the conditional
such that it is equal to:

c = du*conditional( eq(u, 0), 0, conditional(eq(u, 0), 0,
conditional(lt(u,0),-1,1) ) )
J = v*c*dx

which compiles fine and looks correct.

Kristian

> Any idea?
>
> Here's an example based on demo_poisson.py:
>
> from dolfin import *
>
> # Create mesh and define function space
> mesh = UnitSquare(32, 32)
> V = FunctionSpace(mesh, "Lagrange", 1)
>
> # Define Dirichlet boundary (x = 0 or x = 1)
> def boundary(x):
>    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
>
> # Define boundary condition
> u0 = Constant(0.0)
> bc = DirichletBC(V, u0, boundary)
>
> # Define variational problem
> u = Function(V)
> du = TrialFunction(V)
>
> v = TestFunction(V)
> f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
> g = Expression("sin(5*x[0])")
> a = inner(grad(u), grad(v))*dx
> L = f*v*dx + g*v*ds + abs(u)*v*dx
> F = a-L
> J = derivative(F,u,du)
>
> # Compute solution
> problem = VariationalProblem(F, J, bc)
> u = problem.solve()
>
> # Save solution in VTK format
> file = File("poisson.pvd")
> file << u
>
> # Plot solution
> plot(u, interactive=True)
>
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>