Define a piecewise constant source

Question

Hello everybody!

I am looking at the demo_poisson.py demo and instead of using a source term defined by:

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")

I would like to have a source that is 1.0 in the first half of the problem and 0 everywhere else. And I couldn't really see how this could be done using 'Expression' so I used:

def f(x):
    if x[0] < 0.5: 
        return 1.0
    else:
        return 0.0

But it returns an error:

TypeError: unsupported operand type(s) for *: 'function' and
'Argument'

It may be a very basic question so I'm sorry if it is obvious but I haven't been able to find any examples in the demos or answer in the Q&A explaining how to define functions for the variational form, except using 'Expression' or 'Constant'.

Could someone please give me a hint?

Thanks a lot!
Vincent

timm · Answer 1 · Apr 28, 2014

Best answer

Hi - I didn't plug this into the demo_poisson.py problem, but this snippet at least runs without error. I think that you need to subclass expression if you wanted to use it like that. There are probably much cooler ways of doing it, but that was what I came up with …

from dolfin import *

class MyExpression(Expression):
  def eval(self, value, x):
    if x[0] <= 0:
      value[0] = 1.0
    else:
      value[0] = 0.0

  def value_shape(self):
    return (1,)

e = MyExpression()

answered Apr 28, 2014 by timm FEniCS User (2,100 points)
selected Apr 28, 2014 by V_L

Yes, I used:

class MyExpression(Expression):
  def eval(self, value, x):
    if x[0] <= 0.5:
      value[0] = 1.0
    else:
      value[0] = 0.0

  def value_shape(self):
    return (1,)

f = MyExpression()
...
L = f[0]*v*dx + g*v*ds

and it seems to work!

Thanks a lot for this very quick and efficient answer! :)

commented Apr 28, 2014 by V_L FEniCS User (4,440 points)

Hi, in FEniCS 2016.2.0, your snippet produces the error

raise TypeError("Expression needs an integer argument 'degree' if no 'element' is provided.")

How should your snippet be modified to adapt to this change?

thanks

commented May 24, 2017 by Antoine FEniCS Novice (810 points)

In fact one just need to write

f = MyExpression(degree=2)

in FEniCS 2016.2.0

commented May 25, 2017 by Antoine FEniCS Novice (810 points)

I'm still stuck using 1.6.0, but it looks like you figured out that it wanted the degree to be passed in. Sorry it took so long for me to get back to you.

commented May 30, 2017 by timm FEniCS User (2,100 points)

Define a piecewise constant source

1 Answer