Evaluating Solution at Arbitrary Points

General Discussion
1

Dear FreeFEM users and developers,

I am trying to evaluate the value of a solution at a specific point in order to compare it with a theoretical value. If the point lies exactly on a mesh node, the method described in this post works well.

However, if the point does not coincide with a mesh node, interpolation is necessary. I would like to ask whether FreeFEM provides any built-in functionality to perform this kind of interpolated evaluation at arbitrary spatial points.

Thank you in advance for your help.

2

Yes, it does, it’s explained in the manual.

1 Like
3

Thanks. I also got it from AI. I was just surprised since AI couldn’t give a correct answer about FreeFem++ a few months ago.

real value;
value = uh(x_coord, y_coord);

4

That code is wrong, yet another example of wrongful use of AI.

5

Thanks for the comment. I think you mean that the underscore causes the problem. However, u(x, y) should be the correct syntax to obtain the solution value at the point (x, y), right?

The following example solves the equation \\Delta u = 1 on a unit circle with zero boundary conditions. The exact solution is u(x, y) = 0.25, (1 - (x^2 + y^2)). I tested several points, and the numerical results are close to the exact solution.

By the way, I still have trouble finding an explanation of this syntax in the FreeFem++ documentation. I would appreciate it if you could point me to the relevant section or page.

// Define mesh boundary
border C(t=0, 2*pi){x=cos(t); y=sin(t);}
// The triangulated domain Th is on the left side of its boundary
mesh Th = buildmesh(C(50));

// The finite element space defined over Th is called here Vh
fespace Vh(Th, P2);
Vh u, v;// Define u and v as piecewise-P1 continuous functions

// Define a function f
func f= 1;

// Get the clock in second
real cpu=clock();

// Define the PDE
solve Poisson(u, v, solver=LU)
    = int2d(Th)(    // The bilinear part
          dx(u)*dx(v)
        + dy(u)*dy(v)
    )
    - int2d(Th)(    // The right hand side
          f*v
    )
    + on(C, u=0);   // The Dirichlet boundary condition

// Plot the result
// plot(u);
// Display the total computational time
cout << "CPU time = " << (clock()-cpu) << endl;


real xcoord, ycoord, value;
xcoord=0.5;
ycoord=0.5;
value = u(xcoord, ycoord);
func fun=0.25*(1-(x^2+y^2));
cout<<"solution value="<<u(xcoord, ycoord)<<endl;
cout<<"exact value="<<fun(xcoord, ycoord)<<endl;