# Modify a freefem code

Hello, I want to upgrade my code to obtain the velocity and pressure for several points that are independent of the mapping, but I don’t know how to do that. Can someone help me figure it out? Here is my code:

``````// -*- C++ -*-

int[int] fforder2 = [1, 1, 1];

int N = 10000;
int Tolerance = 1.e-4;
real dt = 1e-1; // time step
int meshSize = 100;

mesh Th = square(meshSize, meshSize);

// plot(Th); exit(0);

fespace Vh(Th, P2); // Velocity space
fespace Mh(Th, P1); // Pressure space

Vh u1, u2, u1p, u2p, v1, v2; // Declaration of velocity components and test functions
Mh p, pp, q;               // Declaration of pressure and test function q

real Reynolds = 0.;     // Reynolds
// real  nu       = 1./Reynolds;       // Dynamic viscosity
real U0 = 1.;           // Value imposed on the top side of the square

// Maximum iteration number and declaration of errors

real Err, errP, errU;

// Time loop for the Euler resolution algorithm (first order in time)

real t;

u1 = 0.;
u2 = 0.;
p = 0.;

real nu;
real ReynlodsMax = 10000.;
real incremReynolds = 200;

int NbContinuationReynolds = ReynlodsMax / 25;

//  cout << "  NbContinuationReynolds "  <<  NbContinuationReynolds << endl; exit(0);

for (int j = 0; j <= NbContinuationReynolds; j++) {

Reynolds = Reynolds + incremReynolds;
nu = 1. / Reynolds;

ofstream gnu("Error_Cavity2D_Reynolds=" + Reynolds + ".txt");

for (int i = 0; i <= N; i++) {
t = i * dt;

u1p = u1;
u2p = u2;
pp = p;

solve NS([u1, u2, p], [v1, v2, q]) =

int2d(Th)((u1 * v1 + u2 * v2) / dt) -
int2d(Th)((u1p * v1 + u2p * v2) / dt) +

int2d(Th)(nu * (2 * dx(u1) * dx(v1) + 2 * dy(u2) * dy(v2) +
dy(u1) * dy(v1) + dx(u2) * dx(v2) +
dx(u2) * dy(v1) + dy(u1) * dx(v2))) +

int2d(Th)((u1p * dx(u1) + u2p * dy(u1)) * v1 +
(u1p * dx(u2) + u2p * dy(u2)) * v2) -

int2d(Th)(p * (dx(v1) + dy(v2)) + q * (dx(u1) + dy(u2)) - 1e-7 * p * q) +

on(3, u1 = U0) + on(1, 2, 4, u1 = 0) + on(1, 2, 3, 4, u2 = 0);

errU = sqrt(int2d(Th)(square(u1 - u1p) + square(u2 - u2p)));
errP = sqrt(int2d(Th)(square(p - pp)));
Err = sqrt(errP^2 + errU^2);

// Store error values in an array for convergence analysis:
gnu << t
<< "       " << errU
<< "       " << errP
<< "       " << Err
<< endl;

gnu.flush;

cout << " Display *************.  iter " << i << "    Reynolds === " << Reynolds << "   Err L2 = " << Err << "    " << errU << "  " << errP << endl;

// Save the solution every 10 time steps to avoid filling the disk with vtu files
if (i % 10 == 0) {
savevtk("Sol_Cavity2D_Reynolds=" + Reynolds + "_" + i + ".vtu", Th, p, [u1, u2, 0],
order = fforder2, dataname = "Pressure Velocity", bin = true);
savevtk("velocity_fields=" + Reynolds + "_" + i + ".vtu", Vh, p, [u1, u2, 0],
order = fforder2, dataname = "Pressure Velocity", bin = true);
}

// Termination test as soon as the error is < a given tolerance:
if (Err < Tolerance) {
cout << " ****************** Convergence is achieved, so we stop the calculation at Tfinal ==  **************** " << t << endl;

savevtk("FinalSolution_Cavity2D_Reynolds=" + Reynolds + ".vtu", Th, p, [u1, u2, 0],
order = fforder2, dataname = "P V", bin = true);
break;
}
} // End of the time loop (i)
} // End of the Reynolds continuation loop (j)
``````

IT is very simple,

put all the coordination in a array,

``````int ncoor = 10; // nb of coordinate
real[int] X(ncoor), Y(ncoor);
//..  set the coordinate HERE

macro savedata(filename,p, u1,u2)
{
ofstream ff(filename);
for(int i=0; i< ncoor;++i)
ff << p(X[i],Y(i]) << " " << u1((X[i],Y(i]) << " " << u2(X[i],Y(i]) << endl;
}

``````

to save the solution

``````savedata("path.tex",p,u1,u2);
``````

remark the code is not test so some erreurs are possible.

Thanks a lot. I will test it tomorrow, and I will come back if there are any bugs that I can’t solve myself.