ufl team mailing list archive

[Martin Sandve Alnæs <martinal@xxxxxxxxx>]

Thanks Kristian.

In the abs case, would this change be sufficient?

@@ -388,7 +388,8 @@
     def abs(self, o, a):
         f, fprime = a
         o = self.reuse_if_possible(o, f)
-        oprime = conditional(eq(f, 0), 0, Product(sign(f), fprime))
+        oprime = sign(f)*fprime
         return (o, oprime)

sign returns 0 if f is 0:

def sign(x):
    "UFL operator: Take the sign (+1 or -1) of x."
    return conditional(eq(x, 0), 0, conditional(lt(x, 0), -1, +1))

In the generic conditional case, I'm not sure if this is
easily solvable in the current UFL implementation.

Martin

On 6 August 2011 17:00, Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx> wrote:
> 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.
>>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~ufl
> Post to     : ufl@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~ufl
> More help   : https://help.launchpad.net/ListHelp
>