Including a function in C++ expression

Question

Hi community,

I have been trying to migrate a certain code I had to c++ because the call to python was very slow (very new to C++, answer might be trivial, so please explain in detail :) ). What I had, was a code like this:

        class FExpression(Expression):
        def __init__(self,element, u = Expression(("0.","0."), degree = 0)):
            self.element = element
            self.u = u
        def eval(self, value, x):
            p = x+self.u(x)
            val = alpha*(T(p) - R(x))
            value[0] = val*T(p,1,0)
            value[1] = val*T(p,0,1)

My C++ proposal goes like this:

class Load : public Expression
{
  public:

  Function u;

  Load() : Expression(2) {
      u = std::make_shared<ImageReg::FunctionSpace>(mesh)
      u->vector()->zero();
    }   

 void eval(Array<double>& values, const Array<double>& x) const
 {

    // Constants
    double beta1 = 0.7;
    double beta2 = 1.3;

    // Evaluation points
    double dx = x[0] - 0.5;
    double dy = x[1] - 0.5;

    // Get reference intensity
    double R = exp(-(dx*dx + dy*dy));

    // Warped evaluation points
    double temp [2];
    u->eval(values, x);
    dx += values[0];
    dy += values[1];

    // Target intensity
    double T = exp(-beta1*dx*dx - beta2*dy*dy);
    double k = -2*(R-T)*T;
    values[0] = k*beta1*(-beta1*dx*dx - beta2*dy*dy);
    values[1] = k*beta2*(-beta1*dx*dx - beta2*dy*dy);
    }  
};

The way I see it, I define a Function 'u' in the class, and calling the constructor creates the Function u and sets its coordinates to 0. The problem is that compiling this fails miserably. Below I copy the boundary and main, plus the compiling errors. This is a minimal working version, because in the future R and T will be interpolations of PNG images, but if I manage to load a Function to the expression, then it could be easier to load an array. Thanks in advance.

class DirichletBoundary : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return on_boundary;
  }
};


int main(){
auto mesh = std::make_shared<UnitSquareMesh>(50, 50);
auto V = std::make_shared<ImageReg::FunctionSpace>(mesh);
ImageReg::BilinearForm a(V, V);
ImageReg::LinearForm L(V);

// Parameters
auto boundary_value = std::make_shared<Constant>(0.,0.);
auto mu = std::make_shared<Constant>(1.E5);
auto lamb = std::make_shared<Constant>(0.);
auto alpha = std::make_shared<Constant>(1.);

// Boundary conditions
auto boundary = std::make_shared<DirichletBoundary>();
DirichletBC bc(V, boundary_value, boundary);

// LHS construction
a.mu = mu;
a.lamb = lamb;

// RHS construction
auto f  = std::make_shared<Load>();
auto u0 = std::make_shared<Function>(V);
// u0->vector()->zero();
// f.u = *u0;
L.f = f;

// Solve problem
Function u(V);
solve(a == L, u, bc);
}

[ 50%] Building CXX object CMakeFiles/ImageReg.dir/main.cpp.o
/home/nico/Dropbox/Magister/MFEM Implementation/Implementation2D/cpp/main.cpp: In constructor ‘Load::Load()’:
/home/nico/Dropbox/Magister/MFEM Implementation/Implementation2D/cpp/main.cpp:14:24: error: no matching function for call to ‘dolfin::Function::Function()’
Load() : Expression(2) {
^
In file included from /usr/include/dolfin/ale/MeshDisplacement.h:29:0,
from /usr/include/dolfin/mesh/Mesh.h:35,
from /usr/include/dolfin/mesh/dolfin_mesh.h:11,
from /usr/include/dolfin.h:15,
from /home/nico/Dropbox/Magister/MFEM Implementation/Implementation2D/cpp/main.cpp:1:
/usr/include/dolfin/function/Function.h:113:5: note: candidate: dolfin::Function::Function(const dolfin::Function&, std::size_t)
Function(const Function& v, std::size_t i);
^
/usr/include/dolfin/function/Function.h:113:5: note: candidate expects 2 arguments, 0 provided
/usr/include/dolfin/function/Function.h:102:5: note: candidate: dolfin::Function::Function(const dolfin::Function&)
Function(const Function& v);
^
/usr/include/dolfin/function/Function.h:102:5: note: candidate expects 1 argument, 0 provided
/usr/include/dolfin/function/Function.h:94:5: note: candidate: dolfin::Function::Function(std::shared_ptr, std::__cxx11::string)
Function(std::shared_ptr V,
^
/usr/include/dolfin/function/Function.h:94:5: note: candidate expects 2 arguments, 0 provided
/usr/include/dolfin/function/Function.h:84:5: note: candidate: dolfin::Function::Function(std::shared_ptr, std::shared_ptr)
Function(std::shared_ptr V,
^
/usr/include/dolfin/function/Function.h:84:5: note: candidate expects 2 arguments, 0 provided
/usr/include/dolfin/function/Function.h:72:14: note: candidate: dolfin::Function::Function(std::shared_ptr)
explicit Function(std::shared_ptr V);
^
/usr/include/dolfin/function/Function.h:72:14: note: candidate expects 1 argument, 0 provided
/home/nico/Dropbox/Magister/MFEM Implementation/Implementation2D/cpp/main.cpp:15:53: error: ‘mesh’ was not declared in this scope
u = std::make_shared(mesh)
^

Comments

Did you forget to include the compiling errors?

---

Yes, sorry, it is included now. Thanks for the correction.

---

The Function constructor takes at least one argument (see the docs). You might also want to look at the C++ demo programs if you haven't already.

---

Thanks for the reference. I have looked at some examples, but none include a function in an expression, and it is not quite clear to me why turn everything into a shared pointer, or even less likely to include some external array data. It would be great to have some references, thanks in advance.

---

Maybe the MixedPoisson demo can help you.

---

Thanks johann, I had checked it already but didn't find it very enlightening. I will read more about c++ to understand the language better and try again. I got a small expression running, compiled it with instant in the python code, so it is much faster than it was, but I still have issues figuring what object the expression should have in order to used an external image as data.

Best regards

---

Where do you create a Load() object?
Notes:
1. Your class Load() doesn't see the mesh.
2. u is a Function, but you use std::make_shared for it

---

Hi str, in the end I created a cpp expression and imported it in python, which solved the problem. For the sake of completeness, what I did was save Functions with the fileds R, T, dT and then load them into a cpp expression with

load = '''
    class test : public Expression
        {
        public:

        test() : Expression(2) { }

        void eval(Array<double>& values, const Array<double>& x) const {                
            R->eval(values, x);
            const double R_val = values[0];

            u->eval(values, x);

            Array<double> dx (2);
            dx[0] = values[0] + x[0];                
            dx[1] = values[1] + x[1];

            TT->eval(values, dx);
            const double T_val = values[0];

            dTx->eval(values, dx);
            const double dTx_val = values[0];

            dTy->eval(values, dx);
            const double dTy_val = values[0];

            values[0] = (T_val - R_val)*dTx_val;
            values[1] = (T_val - R_val)*dTy_val;
        }

        std::shared_ptr<const Function> u; 
        std::shared_ptr<const Function> R;
        std::shared_ptr<const Function> TT;
        std::shared_ptr<const Function> dTx;
        std::shared_ptr<const Function> dTy;
    };
    '''

and after creating and instance of it, say f = Expression(cpp_code, element), I set them with

f.R = Rfunction
f.TT = Tfuntion

where it is important to note that I used TT instead of T, because T gave errors (it is apparently a reserved variable). I am closing this to new answers.