Freefem++: Solving poisson equation with numerical function

Question

I am using Freefem++ to solve the poisson equation

Grad^2 u(x,y,z) = -f(x,y,z)

It works well when I have an analytical expression for f, but now I have an f numerically defined (i.e. a set of data defined on a mesh) and I am wondering if I can still use Freefem++.

I.e. typical code (for a 2D problem in this case), looks like the following

mesh Sh= square(10,10); // mesh generation of a square
fespace Vh(Sh,P1); // space of P1 Finite Elements
Vh u,v; // u and v belongs to Vh

func f=cos(x)*y; // analytical function

problem Poisson(u,v)= // Definition of the problem
    int2d(Sh)(dx(u)*dx(v)+dy(u)*dy(v)) // bilinear form
    -int2d(Sh)(f*v) // linear form
    +on(1,2,3,4,u=0); // Dirichlet Conditions
Poisson; // Solve Poisson Equation
plot(u); // Plot the result

I am wondering if I can define f numerically, rather than analytically.

Answer 1

Mesh & space Definition

We define a square unit with Nx=10 mesh and Ny=10 this provides 11 nodes on x axis and the same for y axis.

int Nx=10,Ny=10;
int Lx=1,Ly=1;
mesh Sh= square(Nx,Ny,[Lx*x,Ly*y]); //this is the same as square(10,10)
fespace Vh(Sh,P1); // a space of P1 Finite Elements to use for u definition

Conditions and problem statement

We are not going to use solve but we ll handle matrix (a more sophisticated way to solve with FreeFem).

First we define CL for our problem (Dirichlet ones).

varf CL(u,psi)=on(1,2,3,4,u=0); //you can eliminate border according to your problem state
    Vh u=0;u[]=CL(0,Vh);
    matrix GD=CL(Vh,Vh);

Then we define the problem. Instead of writing dx(u)*dx(v)+dy(u)*dy(v) I suggest to use macro, so we define div as following but pay attention macro finishes by // NOT ;.

 macro div(u) (dx(u[0])+dy(u[1])) //

So Poisson bilinear form becomes:

varf Poisson(u,v)= int2d(Sh)(div(u)*div(v));

After we extract Stifness Matrix

matrix K=Poisson(Vh,Vh);
matrix KD=K+GD; //we add CL defined above

We proceed for solving, UMFPACK is a solver in FreeFem no much attention to this.

set(KD,solver=UMFPACK);

And here what you need. You want to define a value of function f on some specific nodes. I'm going to give you the secret, the poisson linear form.

real[int] b=Poisson(0,Vh);

You define value of the function f at any node you want to do.

b[100]+=20; //for example at node 100 we want that f equals to 20
b[50]+=50; //and at node 50 , f equals to 50 

We solve our system.

u[]=KD^-1*b; 

Finally we get the plot.

plot(u,wait=1);

I hope this will help you, thanks to my internship supervisor Olivier, he always gives to me tricks specially on FreeFem. I tested it, it works very well. Good luck.

Answer 2

The method by afaf works in the case when the function f is a free-standing one. For the terms like int2d(Sh)(f*u*v), another solution is required. I propose (actually I have red it somewhere in Hecht's manual) an approach that covers both cases. However, it works only for P1 finite elements, for which the degrees of freedom are coincided with the mesh nodes.

fespace Vh(Th,P1);
Vh f;

real[int] pot(Vh.ndof);

for(int i=0;i<Vh.ndof;i++){
    pot[i]=something;   //assign values or read them from a file 
}

f[]=pot;