DOLFIN

piecewise defined function

Asked by Achim Schroll

Hi!
My piecewise defined function alpha (see below) just returns zero.
What is wrong? How to define it instead?

Best, Achim

# set parameters
alpha1 = 0.2; alpha2 = 0.8

# create mesh and finite element
mesh = UnitSquare(20, 20)
V = FunctionSpace(mesh, "CG", 1)

# define trial & testfunctions
u_trial = TrialFunction(V); phi = TestFunction(V)

# define diffusion parameter
class alpha(Function):
    def __init__(self, V):
        Function.__init__(self, V)
    def eval(self, v, x):
        v[0] = alpha1
        if x[0] > 0.5: v[0] = alpha2
d = alpha(V)
print " |d| = %g" %(d.vector().norm("l2"))

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Anders Logg
Solved:
Last query:
Last reply:
Revision history for this message
Best Anders Logg (logg) said : 		#1

When using eval() you need to subclass Expression (not Function).

That gives you an Expression that can be used in a variational
problem. But an Expression does not have a vector so if you need the
norm of the vector (I guess just for debugging), you need to project
or interpolate the Expression to a Function.

class Alpha(Expression):
    def eval(self, v, x):
        v[0] = alpha1
        if x[0] > 0.5: v[0] = alpha2

alpha = Alpha()
d = interpolate(alpha, V)

print " |d| = %g" %(d.vector().norm("l2"))

--
Anders

On Fri, Mar 19, 2010 at 08:49:35AM -0000, Achim Schroll wrote:
> New question #104855 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/104855
>
> Hi!
> My piecewise defined function alpha (see below) just returns zero.
> What is wrong? How to define it instead?
>
> Best, Achim
>
> # set parameters
> alpha1 = 0.2; alpha2 = 0.8
>
> # create mesh and finite element
> mesh = UnitSquare(20, 20)
> V = FunctionSpace(mesh, "CG", 1)
>
> # define trial & testfunctions
> u_trial = TrialFunction(V); phi = TestFunction(V)
>
> # define diffusion parameter
> class alpha(Function):
> def __init__(self, V):
> Function.__init__(self, V)
> def eval(self, v, x):
> v[0] = alpha1
> if x[0] > 0.5: v[0] = alpha2
> d = alpha(V)
> print " |d| = %g" %(d.vector().norm("l2"))
>
>

Revision history for this message
Achim Schroll (achim-simula-deactivatedaccount) said : 		#2

Thanks Anders Logg, that solved my question.

Revision history for this message
Johan Hake (johan-hake) said : 		#3

On Friday 19 March 2010 02:30:40 Anders Logg wrote:
> Question #104855 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/104855
>
> Status: Open => Answered
>
> Anders Logg proposed the following answer:
> When using eval() you need to subclass Expression (not Function).
>
> That gives you an Expression that can be used in a variational
> problem. But an Expression does not have a vector so if you need the
> norm of the vector (I guess just for debugging), you need to project
> or interpolate the Expression to a Function.
>
> class Alpha(Expression):
> def eval(self, v, x):
> v[0] = alpha1
> if x[0] > 0.5: v[0] = alpha2
>
> alpha = Alpha()
> d = interpolate(alpha, V)
>
> print " |d| = %g" %(d.vector().norm("l2"))
>
> --
> Anders
>
> On Fri, Mar 19, 2010 at 08:49:35AM -0000, Achim Schroll wrote:
> > New question #104855 on DOLFIN:
> > https://answers.launchpad.net/dolfin/+question/104855
> >
> > Hi!
> > My piecewise defined function alpha (see below) just returns zero.
> > What is wrong? How to define it instead?
> >
> > Best, Achim
> >
> > # set parameters
> > alpha1 = 0.2; alpha2 = 0.8
> >
> > # create mesh and finite element
> > mesh = UnitSquare(20, 20)
> > V = FunctionSpace(mesh, "CG", 1)
> >
> > # define trial & testfunctions
> > u_trial = TrialFunction(V); phi = TestFunction(V)
> >
> > # define diffusion parameter
> > class alpha(Function):
> > def __init__(self, V):
> > Function.__init__(self, V)
> > def eval(self, v, x):
> > v[0] = alpha1
> > if x[0] > 0.5: v[0] = alpha2
> > d = alpha(V)
> > print " |d| = %g" %(d.vector().norm("l2"))

What version of DOLFIN are you using. It has not been possible to overload
Function.eval in PyDOLFIN since 0.9.5.

You might want to upgrade ;)

Johan

Revision history for this message
Achim Schroll (achim-simula-deactivatedaccount) said : 		#4

Hej again,
how can I see my version of DOLFIN and how to upgrade?
Best, Achim

Revision history for this message
Johannes Ring (johannr) said : 		#5

This command will show you the version of DOLFIN:

  pkg-config dolfin --modversion

If you use the Ubuntu packages, you should get upgrades automatically. If you use Dorsal, you should update Dorsal and rebuild.

Revision history for this message
Achim Schroll (achim-simula-deactivatedaccount) said : 		#6

it says dolfin version 0.9.4.
I thought I would get automatic updates, obviously not true :(

Revision history for this message
Johannes Ring (johannr) said : 		#7

So you are using the Ubuntu packages? If so, please follow the instructions at

  http://fenics.org/wiki/Download#Ubuntu_packages

and you should get the latest DOLFIN (and future upgrades).