C++ code for user-defined Expression

Asked by Øystein Sørensen

I have a problem in which I need to define an Expression using "complicated" C++ code. Before making the real code, I am trying to test it, with the following trivial code:

-----------------------------------------------------------

from dolfin import *

code = '''
class SomeExpr : public Expression
{
public:
    SomeExpr() : Expression(1) {}

    void eval(Array<double>& values, const Data& data) const
{
    values[0] = data.x[0];
}

};
'''

mesh = UnitInterval(10)
V = FunctionSpace(mesh, 'CG', 2)

u = TrialFunction(V)
v = TestFunction(V)
f = Expression(code)
a = f*inner(grad(u),grad(v))*dx
L = Constant(0)*v*dx

-----------------------------------------------------------

When running it, I get the error message:
Trying to integrate expression of rank 1 with free indices ().

On the other hand, if I just write
f = Expression('x[0]'),
which is what I try to do in the above code, no error messages are issued.

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Øystein Sørensen
Solved:
Last query:
Last reply:
Revision history for this message
Øystein Sørensen (oystein-sorensen) said :
#1

Seems like the problem was easily solved by replacing Expression(1) by Expression() in the C++ code. But if there exist examples or documentation on user-defined expressions defined in this way, I would be very glad to know.

Revision history for this message
Anders Logg (logg) said :
#2

I think your problem is the way you call the Expression constructor:

  Expression(1)

This creates a vector-valued expression (with only one component).

Try with Expression() and it should work.

--
Anders

On Sun, Sep 12, 2010 at 10:46:36AM -0000, Øystein Sørensen wrote:
> New question #125106 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/125106
>
> I have a problem in which I need to define an Expression using "complicated" C++ code. Before making the real code, I am trying to test it, with the following trivial code:
>
> -----------------------------------------------------------
>
> from dolfin import *
>
> code = '''
> class SomeExpr : public Expression
> {
> public:
> SomeExpr() : Expression(1) {}
>
> void eval(Array<double>& values, const Data& data) const
> {
> values[0] = data.x[0];
> }
>
> };
> '''
>
> mesh = UnitInterval(10)
> V = FunctionSpace(mesh, 'CG', 2)
>
> u = TrialFunction(V)
> v = TestFunction(V)
> f = Expression(code)
> a = f*inner(grad(u),grad(v))*dx
> L = Constant(0)*v*dx
>
> -----------------------------------------------------------
>
> When running it, I get the error message:
> Trying to integrate expression of rank 1 with free indices ().
>
> On the other hand, if I just write
> f = Expression('x[0]'),
> which is what I try to do in the above code, no error messages are issued.
>

--
Anders

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

On Sunday September 12 2010 05:50:15 Øystein Sørensen wrote:
> Question #125106 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/125106
>
> Status: Open => Solved
>
> Øystein Sørensen confirmed that the question is solved:
> Seems like the problem was easily solved by replacing Expression(1) by
> Expression() in the C++ code. But if there exist examples or
> documentation on user-defined expressions defined in this way, I would
> be very glad to know.

Hello!

One of few examples of a more complex compiled C++ expression resides in

  demo/undocumented/tensor-weighted-poisson

and in the docstring of Expression.

Johan

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

On Sunday September 12 2010 09:27:26 Anders Logg wrote:
> Question #125106 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/125106
>
> Anders Logg posted a new comment:
> I think your problem is the way you call the Expression constructor:
>
> Expression(1)
>
> This creates a vector-valued expression (with only one component).
>
> Try with Expression() and it should work.

It looked like this helped. But why are a Vector expression with only one
value different from a scalar expression?

Johan

> --
> Anders
>
> On Sun, Sep 12, 2010 at 10:46:36AM -0000, Øystein Sørensen wrote:
> > New question #125106 on DOLFIN:
> > https://answers.launchpad.net/dolfin/+question/125106
> >
> > I have a problem in which I need to define an Expression using
> > "complicated" C++ code. Before making the real code, I am trying to test
> > it, with the following trivial code:
> >
> > -----------------------------------------------------------
> >
> > from dolfin import *
> >
> > code = '''
> > class SomeExpr : public Expression
> > {
> >
> > public:
> > SomeExpr() : Expression(1) {}
> >
> > void eval(Array<double>& values, const Data& data) const
> >
> > {
> >
> > values[0] = data.x[0];
> >
> > }
> >
> > };
> > '''
> >
> > mesh = UnitInterval(10)
> > V = FunctionSpace(mesh, 'CG', 2)
> >
> > u = TrialFunction(V)
> > v = TestFunction(V)
> > f = Expression(code)
> > a = f*inner(grad(u),grad(v))*dx
> > L = Constant(0)*v*dx
> >
> > -----------------------------------------------------------
> >
> > When running it, I get the error message:
> > Trying to integrate expression of rank 1 with free indices ().
> >
> > On the other hand, if I just write
> > f = Expression('x[0]'),
> > which is what I try to do in the above code, no error messages are
> > issued.
>
> --
> Anders

Revision history for this message
Anders Logg (logg) said :
#5

On Sun, Sep 12, 2010 at 10:32:55PM -0700, Johan Hake wrote:
> On Sunday September 12 2010 09:27:26 Anders Logg wrote:
> > Question #125106 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/125106
> >
> > Anders Logg posted a new comment:
> > I think your problem is the way you call the Expression constructor:
> >
> > Expression(1)
> >
> > This creates a vector-valued expression (with only one component).
> >
> > Try with Expression() and it should work.
>
> It looked like this helped. But why are a Vector expression with only one
> value different from a scalar expression?

Because we have value_rank() and value_dimension(i), and checking for
value_dimension[0] == 1 would be a special case.

--
Anders