stokessolns.edp

verbosity=0;
int n=15;    // mesh size 
int imax=150; // number of snapshots

real L = 1;

border a(t=0,1){x=L*t;y=0;label=1;}; // bas
border b(t=0,1){x=L;y=t;label=2;}; // droite
border c(t=1,0){x=L*t ;y=1;label=3;}; // haut
border d(t=1,0){x=0;y=t;label=4;}; // gauche
mesh Th = buildmesh( a(n) + b(n) + c(n) + d(n) );

// radius and centers
real R= 0.05;

// 
real eps=1e-11;

// Fe space
fespace Uh(Th,P2);
Uh u1,u2,v1,v2;

fespace Ph(Th,P1);
Ph p,q;

// export dofs points
{
  func fx = x;
  func fy = y;
  Uh ux = fx;
  Uh uy = fy;
  cout << "Exporting dofs points in dofs.dat...";
  ofstream mfile("tmp/dofs.dat");
  for (int i = 0; i < Uh.ndof; i++) {
    mfile << ux[][i] << "\t" << uy[][i] << endl;
  }
  cout << "done!" << endl;
}


// energy error
varf ps(u,v) = int2d(Th)(u*v + dx(u)*dx(v) + dy(u)*dy(v));
matrix mPS;
matrix A = ps(Uh,Uh);
{
  ofstream fnrj("tmp/nrj");
  mPS = [[A,0],[0,A]];
  
  cout << "Writing energy matrix...";
  fnrj << mPS;
  cout << "done!" << endl;
}

int seed=454774;
randinit(seed);
{
  ofstream wi("./trial/input.dat");
  ofstream wo("./trial/output.dat");

  // obstacle
  real x1 = 0.5;
  real y1 = 0.5;
  func obstacles = ( (x-x1)*(x-x1) + (y-y1)*(y-y1) <= R*R );

  // snapshots generation
  for (int i = 1 ; i <= imax; i++) {
    real r=0.;
    r=randreal1();
    
    real xl=0.*r + (1-r)*2.; // [0,2]
    r=randreal1();
    real yl=-1.*r + (1-r)*2.; // [-1,2]
    r=randreal1();
    
    real xr=-1.*r + (1-r)*2.; // [-1,1]
    r=randreal1();
    real yr=-1.*r + (1-r)*2.; // [-1,2]
    

    solve stokes([u1,u2,p],[v1,v2,q]) = int2d(Th)( dx(u1)*dx(v1)+dy(u1)*dy(v1) + dx(u2)*dx(v2)+dy(u2)*dy(v2)
						   //+(1/eps)*obstacles*(u1*v1+u2*v2)
						   - p*q*(0.000001)
						   - p*(dx(v1)+dy(v2))
						   - q*(dx(u1)+dy(u2))) 
      + on(4,u1=xl,u2=yl)
      + on(2,u1=xr,u2=yr)
      + on(1,3, u1=0., u2=0.);

    //
    cout <<"INPUT\t" << xl << " " << yl << " " << xr << " " << yr << endl;
    wi << xl << " " << yl << " " << xr << " " << yr << endl;

    // energy
    /*
    {
      real[int] Aux = A*u1[];
      real[int] Auy = A*u2[];
      real uxAux = u1[]'*Aux;
      real uyAuy = u2[]'*Auy;
      
      real uAu = sqrt(uxAux + uyAuy);
      cout << " error = " << uAu;
    }
    */
    
    cout << " OUTPUT - Writing u1 and u2..." << endl;
    for (int i = 0; i < Uh.ndof; i++) {
      u1[][i] = tanh(u1[][i]);
      wo << u1[][i] << " ";
    }
    for (int i = 0; i < Uh.ndof; i++) {
      u2[][i] = tanh(u2[][i]);
      wo << u2[][i] << " ";
    }
    wo << endl;
    
    //Uh normuh = sqrt(u1*u1+u2*u2);
    plot([u1,u2],fill=1,value=1,cmm="Obstacle at x="+x1+" y="+y1+" i="+i+"/"+imax);
  }


  /*
  {
    real xl=1;
    real yl=1.;
    real xr=1.;
    real yr=1.;
    

    solve stokes([u1,u2,p],[v1,v2,q]) = int2d(Th)( dx(u1)*dx(v1)+dy(u1)*dy(v1) + dx(u2)*dx(v2)+dy(u2)*dy(v2)
						   +(1/eps)*obstacles*(u1*v1+u2*v2)
						   - p*q*(0.000001)
						   - p*(dx(v1)+dy(v2))
						   - q*(dx(u1)+dy(u2))) 
      + on(4,u1=xl,u2=yl)
      + on(2,u1=xr,u2=yr)
      + on(1,3, u1=0., u2=0.);

    //
    cout <<"TEST :\t" << xl << " " << yl << " " << xr << " " << yr << endl;

    for (int i = 0; i < Uh.ndof; i++) {
      u1[][i] = tanh(u1[][i]);
      u2[][i] = tanh(u2[][i]);
    }
    
    //Uh normuh = sqrt(u1*u1+u2*u2);
    plot([u1,u2],fill=1,value=1);
  }  
  */
  
}
cout << "dofs: vel=" << 2*Uh.ndof << " pre=" << Ph.ndof << endl;