This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Define a piecewise constant source

+1 vote

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

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

1 Answer

+2 votes
 
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! :)

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

In fact one just need to write

f = MyExpression(degree=2)

in FEniCS 2016.2.0

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.

...