using abs in forms
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]'}
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(
g = Expression(
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 = VariationalProb
u = problem.solve()
# Save solution in VTK format
file = File("poisson.pvd")
file << u
# Plot solution
plot(u, interactive=True)
Question information
- Language:
- English Edit question
- Status:
- Solved
- For:
- DOLFIN Edit question
- Assignee:
- No assignee Edit question
- Solved by:
- Chaffra Affouda
- Solved:
- Last query:
- Last reply:
|
Revision history for this message
|
#1 |
Maybe there is a better solution, but the following seems to work
L = f*v*dx + g*v*ds + sqrt(u**2 + 1.0e-15)*v*dx
and you also need to use
problem.solve(u)
|
Revision history for this message
|
#2 |
Nice trick. This works, but I am still wondering why abs bugs.
Thanks Praveen.
|
Revision history for this message
|
#3 |
On 4 August 2011 14:56, Chaffra <email address hidden> wrote:
> New question #166940 on DOLFIN:
> https:/
>
> 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(
u = Coefficient(V)
du = TrialFunction(V)
v = TestFunction(V)
L1 = abs(u)*v*dx
J = derivative(L1,u,du)
Found Argument in Conditional(
Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
(), {}), Product(
Space(2)), 1, None), 1),
Conditional(
Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
Conditional(
Space(2)), 1, None), 0), Zero((), (), {})), IntValue(-1, (), (), {}),
IntValue(1, (), (), {}))))), this is an invalid expression.
Traceback (most recent call last):
File "/home/
sys.
File "/home/
compile_
File "/home/
line 150, in compile_form
analysis = analyze_
File "/home/
line 64, in analyze_forms
common_cell) for form in forms)
File "/home/
line 64, in <genexpr>
common_cell) for form in forms)
File "/home/
line 151, in _analyze_form
ffc_
File "/home/
line 332, in compute_
parts = compute_
File "/home/
line 317, in compute_
res = transform_
File "/home/
line 799, in transform_
integrand = transform(
File "/home/
line 313, in _transform
e, provides = pe.visit(e)
File "/home/
line 144, in visit
r = h(o, *map(self.visit, o.operands()))
File "/home/
line 148, in visit
r = h(o)
File "/home/
line 68, in expr
error("Found Argument in %s, this is an invalid expression." % repr(x))
File "/home/
line 148, in error
raise self._exception
ufl.log.
Conditional(
Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
Product(
1, None), 1), Conditional(
Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
(), {}), Conditional(
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(
u = Coefficient(V)
du = TrialFunction(V)
v = TestFunction(V)
L1 = abs(u)*v*dx
J1 = derivative(L1,u,du)
print expand_
{ 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(
conditional(
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(
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(
> g = Expression(
> 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 = VariationalProb
> 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.
>
|
Revision history for this message
|
#4 |
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_
- 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 <email address hidden> wrote:
> On 4 August 2011 14:56, Chaffra <email address hidden> wrote:
>> New question #166940 on DOLFIN:
>> https:/
>>
>> 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(
> u = Coefficient(V)
> du = TrialFunction(V)
> v = TestFunction(V)
> L1 = abs(u)*v*dx
> J = derivative(L1,u,du)
>
> Found Argument in Conditional(
> Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
> (), {}), Product(
> Space(2)), 1, None), 1),
> Conditional(
> Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
> Conditional(
> Space(2)), 1, None), 0), Zero((), (), {})), IntValue(-1, (), (), {}),
> IntValue(1, (), (), {}))))), this is an invalid expression.
> Traceback (most recent call last):
> File "/home/
> sys.exit(
> File "/home/
> compile_
> File "/home/
> line 150, in compile_form
> analysis = analyze_
> File "/home/
> line 64, in analyze_forms
> common_cell) for form in forms)
> File "/home/
> line 64, in <genexpr>
> common_cell) for form in forms)
> File "/home/
> line 151, in _analyze_form
> ffc_assert(
> File "/home/
> line 332, in compute_
> parts = compute_
> File "/home/
> line 317, in compute_
> res = transform_
> File "/home/
> line 799, in transform_
> integrand = transform(
> File "/home/
> line 313, in _transform
> e, provides = pe.visit(e)
> File "/home/
> line 144, in visit
> r = h(o, *map(self.visit, o.operands()))
> File "/home/
> line 148, in visit
> r = h(o)
> File "/home/
> line 68, in expr
> error("Found Argument in %s, this is an invalid expression." % repr(x))
> File "/home/
> line 148, in error
> raise self._exception
> ufl.log.
> Conditional(
> Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
> Product(
> 1, None), 1), Conditional(
> Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
> (), {}), Conditional(
> 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(
> u = Coefficient(V)
> du = TrialFunction(V)
> v = TestFunction(V)
> L1 = abs(u)*v*dx
> J1 = derivative(L1,u,du)
> print expand_
>
> { 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(
> conditional(
> 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(
> 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(
>> g = Expression(
>> 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 = VariationalProb
>> 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:/
> Post to : <email address hidden>
> Unsubscribe : https:/
> More help : https:/
>
|
Revision history for this message
|
#5 |
On 8 August 2011 15:19, Martin Sandve Alnæs <email address hidden> wrote:
> 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_
> - 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.
For the given problem the above change should do the trick.
If other combinations of conditional turns out to crash compilation we
might need to look at other ways of solving the issue as you point
out.
It's difficult to predict so let's wait and see if some one hits a
corner case we can investigate.
Kristian
> Martin
>
> On 6 August 2011 17:00, Kristian Ølgaard <email address hidden> wrote:
>> On 4 August 2011 14:56, Chaffra <email address hidden> wrote:
>>> New question #166940 on DOLFIN:
>>> https:/
>>>
>>> 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(
>> u = Coefficient(V)
>> du = TrialFunction(V)
>> v = TestFunction(V)
>> L1 = abs(u)*v*dx
>> J = derivative(L1,u,du)
>>
>> Found Argument in Conditional(
>> Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
>> (), {}), Product(
>> Space(2)), 1, None), 1),
>> Conditional(
>> Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
>> Conditional(
>> Space(2)), 1, None), 0), Zero((), (), {})), IntValue(-1, (), (), {}),
>> IntValue(1, (), (), {}))))), this is an invalid expression.
>> Traceback (most recent call last):
>> File "/home/
>> sys.exit(
>> File "/home/
>> compile_
>> File "/home/
>> line 150, in compile_form
>> analysis = analyze_
>> File "/home/
>> line 64, in analyze_forms
>> common_cell) for form in forms)
>> File "/home/
>> line 64, in <genexpr>
>> common_cell) for form in forms)
>> File "/home/
>> line 151, in _analyze_form
>> ffc_assert(
>> File "/home/
>> line 332, in compute_
>> parts = compute_
>> File "/home/
>> line 317, in compute_
>> res = transform_
>> File "/home/
>> line 799, in transform_
>> integrand = transform(
>> File "/home/
>> line 313, in _transform
>> e, provides = pe.visit(e)
>> File "/home/
>> line 144, in visit
>> r = h(o, *map(self.visit, o.operands()))
>> File "/home/
>> line 148, in visit
>> r = h(o)
>> File "/home/
>> line 68, in expr
>> error("Found Argument in %s, this is an invalid expression." % repr(x))
>> File "/home/
>> line 148, in error
>> raise self._exception
>> ufl.log.
>> Conditional(
>> Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
>> Product(
>> 1, None), 1), Conditional(
>> Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
>> (), {}), Conditional(
>> 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(
>> u = Coefficient(V)
>> du = TrialFunction(V)
>> v = TestFunction(V)
>> L1 = abs(u)*v*dx
>> J1 = derivative(L1,u,du)
>> print expand_
>>
>> { 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(
>> conditional(
>> 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(
>> 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(
>>> g = Expression(
>>> 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 = VariationalProb
>>> 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:/
>> Post to : <email address hidden>
>> Unsubscribe : https:/
>> More help : https:/
>>
>
|
Revision history for this message
|
#6 |
A workaround for such future problems can be
c = conditional(..., 1, 0)
f = c*t + (1-c)*f
which should work fine.
Martin
Den 8. aug. 2011 16.06 skrev "Kristian B. Ølgaard" <
<email address hidden>> følgende:
> Question #166940 on DOLFIN changed:
> https:/
>
> Kristian B. Ølgaard posted a new comment:
> On 8 August 2011 15:19, Martin Sandve Alnæs <email address hidden> wrote:
>> 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_
>> - 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.
>
> For the given problem the above change should do the trick.
> If other combinations of conditional turns out to crash compilation we
> might need to look at other ways of solving the issue as you point
> out.
> It's difficult to predict so let's wait and see if some one hits a
> corner case we can investigate.
>
> Kristian
>
>> Martin
>>
>> On 6 August 2011 17:00, Kristian Ølgaard <email address hidden> wrote:
>>> On 4 August 2011 14:56, Chaffra <email address hidden>
wrote:
>>>> New question #166940 on DOLFIN:
>>>> https:/
>>>>
>>>> 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(
>>> u = Coefficient(V)
>>> du = TrialFunction(V)
>>> v = TestFunction(V)
>>> L1 = abs(u)*v*dx
>>> J = derivative(L1,u,du)
>>>
>>> Found Argument in Conditional(
>>> Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
>>> (), {}), Product(
>>> Space(2)), 1, None), 1),
>>> Conditional(
>>> Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
>>> Conditional(
>>> Space(2)), 1, None), 0), Zero((), (), {})), IntValue(-1, (), (), {}),
>>> IntValue(1, (), (), {}))))), this is an invalid expression.
>>> Traceback (most recent call last):
>>> File "/home/
>>> sys.exit(
>>> File "/home/
>>> compile_
>>> File
"/home/
>>> line 150, in compile_form
>>> analysis = analyze_
common_cell)
>>> File
"/home/
>>> line 64, in analyze_forms
>>> common_cell) for form in forms)
>>> File
"/home/
>>> line 64, in <genexpr>
>>> common_cell) for form in forms)
>>> File
"/home/
>>> line 151, in _analyze_form
>>> ffc_assert(
>>> File
"/home/
>>> line 332, in compute_
>>> parts = compute_
>>> File
"/home/
>>> line 317, in compute_
>>> res = transform_
>>> File
"/home/
>>> line 799, in transform_
>>> integrand = transform(
>>> File
"/home/
>>> line 313, in _transform
>>> e, provides = pe.visit(e)
>>> File
"/home/
>>> line 144, in visit
>>> r = h(o, *map(self.visit, o.operands()))
>>> File
"/home/
>>> line 148, in visit
>>> r = h(o)
>>> File
"/home/
>>> line 68, in expr
>>> error("Found Argument in %s, this is an invalid expression." %
repr(x))
>>> File
"/home/
>>> line 148, in error
>>> raise self._exception
>>> ufl.log.
>>> Conditional(
>>> Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
>>> Product(
>>> 1, None), 1), Conditional(
>>> Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
>>> (), {}), Conditional(
>>> 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(
>>> u = Coefficient(V)
>>> du = TrialFunction(V)
>>> v = TestFunction(V)
>>> L1 = abs(u)*v*dx
>>> J1 = derivative(L1,u,du)
>>> print expand_
>>>
>>> { 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(
>>> conditional(
>>> 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(
>>> 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(
0.02)")
>>>> g = Expression(
>>>> 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 = VariationalProb
>>>> 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:/
>>> Post to : <email address hidden>
>>> Unsubscribe : https:/
>>> More help : https:/
>>>
>>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
