Search code examples
numerical-methodsfinite-element-analysisfreefem++

Freefem++: Solving poisson equation with numerical function


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.


Solution

  • 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.