Function as parameter for another form

Question

Hello,

I would like to calculate the volume under a function I have obtained by solving a linear problem:

Suppose I have a and L defined in an external .ufl file. I compute my solution u with

Poisson::BilinearForm a(V, V);
Poisson::LinearForm L(V);

Function u(V);

solve(a == L, u, bc);

Now, I have an integral.ufl file

element = FiniteElement("Lagrange", triangle, 2)

g = Coefficient(element)
integral = g * dx

forms = [integral]

which is supposed to evaluate the volume under a function g.
In the .cpp file I have

integral::Form_integral integrator(mesh);
integrator.g = u;
double area = dolfin::assemble(integrator);

But this confuses the compiler, since I get the error

error: no match for ‘operator=’ (operand types are ‘dolfin::CoefficientAssigner’ and ‘dolfin::Function’)
       integrator.g = u;
                    ^
note: candidate: void dolfin::CoefficientAssigner::operator=(std::shared_ptr<const dolfin::GenericFunction>)
         void operator= (std::shared_ptr<const GenericFunction> coefficient);
              ^
/usr/include/dolfin/function/CoefficientAssigner.h:50:10: note:   no known conversion for argument 1 from ‘dolfin::Function’ to ‘std::shared_ptr<const dolfin::GenericFunction>’

However, the whole thing works if I set

auto one = std::make_shared<One>();
integral::Form_integral integrator(mesh);
integrator.g = one;
double area = dolfin::assemble(integrator);

where

class One : public Expression
{
    public:
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = 1.0;
  }
};

Any idea how I can convert my Function u to an acceptable parameter for the integrator form?

Thanks a lot!

Answer 1

I think u needs to be a shared pointer. This based on the error message:

note: candidate: void dolfin::CoefficientAssigner::operator=(std::shared_ptr<const dolfin::GenericFunction>)
         void operator= (std::shared_ptr<const GenericFunction> coefficient);
              ^

try:

auto u = std::make_shared<Function>(V);

solve(a == L, *u, bc);

Comments

Yes, that was the problem.

Thank you very much!