I'm trying to solve the steady-state heat equation with a space-dependent thermal diffusivity in FreeFem. I.e. $$\nabla\cdot(\kappa(x,y)\nabla T) = 0$$ I'm hoping to read in from a file locations $(x_i,y_i)$ around which diffusivity changes. So, for the sake of this example, you might think of something of the form $$\kappa(x,y) = \sum_i \frac 1{\|(x,y)-(x_i,y_i)\|}$$ So what I would like would be something along the lines of (this obviously won't compile, because I haven't declared all of the variables etc. I'm just outlining)

func real kappa() {
   for (i=0;i<n;++i)
      k += 1/sqrt((x-x(i)) + (y-y(i))^2);
   return k;

problem heat(uh,vh,solver=GMRES) =
   int2d(Th)( (dx(uh)*dx(vh) + dy(uh)*dy(vh)) * kappa ) // * kappa doesn't work!
   + on(1,uh=1) + on(2,uh=0);

Where I've taken Dirichlet b.c. of 1 on one side, and 0 elsewhere. The problem is that I can't write this function on one line (because it needs a loop), and FreeFem doesn't let me call a user-defined function within the problem definition like this. So either I'm missing something, or I need another way to do this. Pointers?


Okay, I think the reason is that the objects dx(uh) and co. are of type Vh and so I need an object of the same type. Which means I must build my function values onto that type.

  • $\begingroup$ What might be wrong is the name associated to the points $(x(i),y(i))$. This is also the name of the variables $x$ and $y$. Have you tried $(pointx(i),pointy(i))$ instead? I do not think your edit is useful here. By the way, what is the error message? $\endgroup$
    – pluton
    May 4, 2016 at 0:43
  • $\begingroup$ I think you should create the freefem tag. $\endgroup$
    – pluton
    May 4, 2016 at 0:50
  • 1
    $\begingroup$ Actually, software-specific questions are off-topic here; you should rather ask them on the FreeFEM++ mailing list. $\endgroup$ May 4, 2016 at 10:38

1 Answer 1


Yeah, I got the type wrong. This is possible, but the function needs to be evaluated over an object of the same type. For example,

mesh Th=buildmesh(...);
fespace Vh(Th,P1);
Vh uh, vh;

Vh kappa = 0;
for(int i = 0; i < n; ++i) {
   kappa += 1/sqrt((x - c(i,0))^2 + (y - c(i,1))^2);

problem heat(uh,vh) =
int2d(Th)((dx(uh) * dx(vh) + dy(uh) * dy(vh)) * kappa)
+ on(1,uh=1) + on(2,uh=0);

Not the answer you're looking for? Browse other questions tagged or ask your own question.