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!

Answer 1

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!

Comments

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?

---

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!

---

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);

---

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!

---

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?