This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

[dolfin 1.2.0] How to sum functions and multiply function by scalar?

+2 votes

How can i get function f3 like that:

Mesh mesh("./mesh.xml");
Function f1 (mesh,"./f1.xml");
Function f2 (mesh,"./f2.xml");
Function f3 (mesh);
double a = 2;
//f3=a*f1+f2;

using dolfin 1.2.0 via C++ interface?

related to an answer for: How to sum functions?
asked Jul 18, 2013 by mgrviper FEniCS Novice (370 points)

1 Answer

+3 votes

Generally

f3 = project(a*f1 + f2, V)

If all the functions are from same space then

f3.vector()[:] += f2.vector()
f3.vector().axpy(a, f1.vector())

should work. If you can afford mutate f2 then you can do this in one operation obviously.

answered Jul 18, 2013 by Jan Blechta FEniCS Expert (51,420 points)

Oh, I didn't notice C++. For first approach you define linear and bilinear form and solve. Second approach should work in C++ without [:].

Thanks for answer!
When i try to sum functions like this:

  f3.vector() += f2.vector();

I get an error: "no match for ‘operator+=’ in ‘dolfin::Function::vector()() += dolfin::Function::vector()()’"

and when i try this:

f3.vector()= f1.vector() + f2.vector();

I get similar error: "no match for ‘operator+’ in ‘dolfin::Function::vector()() + dolfin::Function::vector()()’"

Also f3.vector() doesn't have axpy member so

f3.vector().axpy(a, f1.vector());

gives me an error: "error: ‘class boost::shared_ptr’ has no member named ‘axpy’"

Then i tried first approach but failed as well. How exactly first approach would look in C++? I suggested:

f3 = project(a*f1 + f2, V);

But get an error: "‘project’ was not declared in this scope".

P.S.: I'm using current stable version of dolfin - 1.2.0 . If upgrade to newer version is only sane solution to my problems i will upgrade, but currently i want to stick with 1.2.0 for compatibility purposes.

In C++ Function::vector() returns pointer so you need to do rather f->vector() to get the Vector object which defines Vector::operator+= and so on. You should refer to documentation of Function and GenericVector or post short and running example. You could also consider using python iterface which is simpler from a programming perspective.

For first approach you need to solve a linear variational problem
$$ \int f_3 v \,\mathrm{d}x = \int (af_1+f_2) v \,\mathrm{d}x \qquad \forall v \in V$$
with suitable space $V$.

Thanks for pointing it out!
Finally i implement

f3 = a*f1+f2;

like this:

f3.vector()->operator =(0);
f3.vector()->axpy(a,f1.vector()->operator *=(1));
f3.vector()->operator +=(f2.vector()->operator *=(1));

it would be nice to produce cleaner and simpler code, but i don't know how.
Any suggestions?

What about

*f3.vector() = *f2.vector();
f3.vector()->axpy(a, *f1.vector());
...