Integral of Function in C++

Question

Dear all,

I am porting a software from Python to C++, but I am missing one little piece: how to integrate a function in C++.

In Python, I did the following:

theta = Function(V) 
number = assemble(q(u, theta)*dx)

Both theta and u are defined over a VectorFunctionSpace, and q returns a scalar (it's basically an inner product). The result is, obviously, a single number.

How can I do this in C++?

I'm not particularly skilled in UFL to know this, but I see that assmeble could do exactly what I'd need. However, I cannot imagine the UFL part...

Thanks for any hint!

david.bernstein · Answer 1 · Jul 16, 2015

Best answer

In your UFL file add the lines

f = Coefficient(element)
functionTotal = f*dx
forms = [your other foms, functionTotal]

Here "element" is from the (already existing) line element = FiniteElement(...)

The "forms" line will export all your forms along with the new one "functionTotal", i.e., make them available in the header file that is generated when the UFL file is compiled.

To integrate a function, in your C++ do the following:
UFL_file_name::Form_functionTotal total(mesh);
total.f = function_you_want_to_integrate
double result = assemble(total);

I'm sure there is a less clumsy way of doing this but it will work!

answered Jul 16, 2015 by david.bernstein FEniCS User (2,000 points)
selected Jul 17, 2015 by senseiwa

So, as I understand, the integration will be on the C++ side, with a function as an expression like this:

// Just a zero constant expression
class ZeroExpression : public Expression
{
public:

    ZeroExpression() : Expression(2) { }

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

Am I right?

commented Jul 17, 2015 by senseiwa FEniCS User (2,620 points)

I've tried, but failed. Reading the documentation, I've created a class derived from Expression (Expression(): Create scalar expression.):

class ConstantExpression : public Expression
{
public:

    ConstantExpression(double value) : Expression(), value_(value)
    {
        // NOP
    }

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

public:
    double value_;
};

Then the UFL:

element_t = FiniteElement("Discontinuous Lagrange", cell, 0)
testf = Coefficient(element_t)
tot_f = testf * dx

Finally the C++:

auto integral_form = ElasticityPoisson2D::Form_tot_f(mesh);

integral_form.testf = ConstantExpression(1.0);

double area = assemble(integral_form);

I have a runtime error saying Reason: Missing eval() function (must be overloaded)..

Am I missing something here?

Thanks!

commented Jul 17, 2015 by senseiwa FEniCS User (2,620 points)

Hmm, I don't see the error. In the C++ code have you tried

auto integral_form = ElasticityPoisson2D::Form_tot_f(mesh);

auto f = ConstantExpression(1.0);
integral_form.testf = f;

double area = assemble(integral_form);

commented Jul 17, 2015 by david.bernstein FEniCS User (2,000 points)

Yes, it works!

Why does an automatic storage variable work, while and a temporary doesn't? In theory, they are equivalent.

Can you explain this weird behavior? I am missing something here...

Thanks a lot!

commented Jul 17, 2015 by senseiwa FEniCS User (2,620 points)

Great! I don't know to be honest, it seems like it should work the way you had it. A quirk (feature) of the form compiler maybe?

commented Jul 17, 2015 by david.bernstein FEniCS User (2,000 points)

Integral of Function in C++

1 Answer