Modify a freefem code

General Discussion
1

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++ -*-

load "iovtk"
int[int] fforder2 = [1, 1, 1];

// load "UMFPACK64"
// load "msh3"
load "medit"

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)
2

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.

3

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