FreeFEM webplot

border ba(t=0,1.0) { x=t;   y=0;   label=1; };
border bb(t=0,0.5) { x=1;   y=t;   label=1; };
border bc(t=0,0.5) { x=1-t; y=0.5; label=1; };
border bd(t=0.5,1) { x=0.5; y=t;   label=1; };
border be(t=0.5,1) { x=1-t; y=1;   label=1; };
border bf(t=0.0,1) { x=0;   y=1-t; label=1; };

mesh Th = buildmesh( ba(6)+bb(4)+bc(4)+bd(4)+be(4)+bf(6) );

fespace Vh(Th,P1);//P1 P1nc
Vh u,v;
func f = 1;

solve Function(u,v)
  = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
    - int2d(Th)(f*v) + on(1,u=0);

fespace Vh0(Th,P0);
Vh0 u0 = u;

load "webplot"
webplot(Th);
webplot(u,Th);
webplot(u0,Th);
server();

macro  e11(u) (dx(u#1)) //  $\varepsilon_{11}$
macro  e22(u) (dy(u#2)) //  $\varepsilon_{22}$
macro  e12(u) ((dy(u#1) + dx(u#2))/2.0) 
macro  div(u) (dx(u#1) + dy(u#2))

int i=0;
real t=0;
real nu = 1.0/400.0;
real dt=0.1;        // time step
int imax=100;         // max iteration steps
int iplt0=imax/40;    // interval for plot
int iplt1=imax/10;    // interval for ps file
real alpha=1.0/dt;


int wall=1, slide=2, inflow=4, outflow=5; 
border  ba(t=0,1.0){x=t*5.0-1.0;y=-1.0;label=slide;}; 
border  bb(t=0,1.0){x=4.0;y=2.0*t-1.0;label=outflow;}; 
border  bc(t=0,1.0){x=4.0-5.0*t;y=1.0;label=slide;};  
border  bd(t=0,1.0){x=-1.0;y=1.0-2.0*t;label=inflow;}; 
border cc(t=0,2*pi){x=cos(t)*0.25;y=sin(t)*0.25;label=wall;};
mesh Th=buildmesh(ba(30)+bb(30)+bc(30)+bd(30)+cc(-30));

fespace Vh(Th,[P2,P2]),Qh(Th,P1),Xh(Th,P2);

Vh [u1,u2],[v1,v2],[up1,up2];
Qh p,q;
Xh psi,phi,phii;

problem  NS([u1,u2,p], [v1,v2,q],solver=UMFPACK,init=i) =
    int2d(Th)(  // $\int_{\Omega}\{$
           alpha * (u1*v1 + u2*v2)  // $\frac{1}{\delta t}(\vec{u}\cdot \vec{v})$
       + nu*2.0*(e11(u)*e11(v)+2.0*e12(u)*e12(v)+e22(u)*e22(v))
	    - div(v)*p-div(u)*q  //$-p\textrm{div}\vec{v}-q\textrm{div}\vec{u}$
    	- p*q*1.e-10 )  // $-\epsilon pq$
     - int2d(Th)( alpha * convect([up1,up2],-dt,up1)*v1
                + alpha * convect([up1,up2],-dt,up2)*v2 )
        + on(slide,u2=0) 
        + on(inflow,u1=1.0-y*y,u2=0) 
        + on(wall,u1=0,u2=0); 

phii = y-y*y*y/3.0;
problem streamlines(phi,psi,solver=Crout,init=i) = // from cavity.edp
  int2d(Th)(dx(phi)*dx(psi)+dy(phi)*dy(psi))
   - int2d(Th)((dx(u2)-dy(u1))*psi)
   + on(slide,phi=phii)
   + on(inflow,phi=phii)
   + on(wall,phi=0);

real  area = int2d(Th)(1.0);
real  meanp; 

[u1,u2] = [0,0];
load "webplot"
for ( i = 1; i <= imax; i++) {
   t=dt*i;
   [up1,up2]=[u1,u2]; 
   NS; 
   meanp = int2d(Th)(p) / area; p = p - meanp;  // $\int_{\Omega}p^{n}\, d\Omega=0$
   streamlines;
   webplot(phi,Th,cmm="Streamlines t="+t);
}

server();

load "msh3"
// Parameters
int nn = 10; // Mesh quality

// Mesh
int[int] labs = [1, 2, 2, 1, 1, 2]; // Label numbering
mesh3 Th = cube(nn, nn, nn, label=labs);
// Remove the ]0.5,1[^3 domain of the cube
Th = trunc(Th, (x < 0.5) | (y < 0.5) | (z < 0.5), label=1);

// Fespace
fespace Vh(Th, P1);
Vh u, v;

// Macro
macro Grad(u) [dx(u), dy(u), dz(u)] //

// Define the weak form and solve
solve Poisson(u, v, solver=CG)
    = int3d(Th)(
          Grad(u)' * Grad(v)
    )
    -int3d(Th)(
          1 * v
    )
    + on(1, u=0)
    ;

// Plot
load "webplot"
webplot(Th);
webplot(u, Th);
server();

//  run with MPI:  ff-mpirun -np 4 script.edp
// NBPROC 4

load "PETSc"                        // PETSc plugin
macro dimension()2// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
func Pk = P2;                       // finite element space

real dt = getARGV("-dt", 0.1);                                 // time step
real mu = 0.01;
int iMax = getARGV("-iMax", 10);                                // number of iterations

int[int] l = [1, 1, 1, 1];
mesh Th = square(getARGV("-global", 40), getARGV("-global", 40), label = l); // global mesh
Mat A;
buildMat(Th, getARGV("-split", 1), A, Pk, mpiCommWorld)


fespace Wh(Th, Pk);                 // local finite element space
real[int] rhs(Wh.ndof);             // local right-hand side
matrix<real> M;                     // local mass matrix
{                                   // local weak form
    fespace Ph(Th, P0);
    Ph f = x < 0.25 ? 1.0e+3 : 1.0e-3;
    varf vPb(u, v) = int2d(Th)(u * v + dt * (grad(u)' * grad(v))) + int2d(Th)(dt * f *  v) + on(1, u = 0.0);
    A = vPb(Wh, Wh, tgv = -1);
    rhs = vPb(0, Wh, tgv = -1);
    varf vPbM(u, v) = int2d(Th)(u * v);
    M = vPbM(Wh, Wh);
}

set(A, sparams = "-ksp_view");
Wh<real> u = 0.0;                   // local solution
plotDmesh(Th, cmm = "Global solution")

for(int i = 0; i < iMax; ++i) {
    real[int] newRhs(rhs.n);
    newRhs = M * u[];
    newRhs += rhs;

    u[] = A^-1 * newRhs;

    macro def(u)u//
    // plotMPI(Th, u, Pk, def, real, cmm = "Global solution")
    plotD(Th, u, cmm = "Global solution")

}
serverMPI()

//  run with MPI:  ff-mpirun -np 4 script.edp
// NBPROC 4

load "PETSc"                        // PETSc plugin
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

macro grad(u)[dx(u), dy(u), dz(u)]// // three-dimensional gradient
func Pk = P1;                       // finite element space

int[int] LL = [1,2, 1,2, 1,1];
mesh3 Th = cube(getARGV("-global", 10), getARGV("-global", 10), getARGV("-global", 10), [x, y, z], label = LL);    // global mesh
Mat A;
buildMat(Th, getARGV("-split", 1), A, Pk, mpiCommWorld)

fespace Wh(Th, Pk);                 // local finite element space
varf vPb(u, v) = int3d(Th)(grad(u)' * grad(v)) + int3d(Th)(v) + on(1, u = 0.0);
real[int] rhs = vPb(0, Wh);

set(A, sparams = "-ksp_view");
Wh<real> u;                         // local solution

A = vPb(Wh, Wh);
u[] = A^-1 * rhs;

real[int] err = A * u[];            // global matrix-vector product
real[int] transpose = A' * u[];
exchange(A, rhs, scaled = true);
err -= rhs;
plotDmesh(Th, cmm = "Mesh")

macro def(u)u//
plotMPI(Th, u, Pk, def, real, cmm = "Global solution")
u[] = err;
plotMPI(Th, u, Pk, def, real, cmm = "Global residual")

Wh<real> Rb[1];
Rb[0] = 1;
set(A, sparams = "-pc_type gamg -ksp_type gmres -ksp_max_it 200", nearnullspace = Rb);
u[] = 0.0;
u[] = A^-1 * rhs;
plotMPI(Th, u, Pk, def, real, cmm = "Global solution")
serverMPI();

load "msh3" 
load "tetgen"
load "webplot"

real R = 3, r=1; 
real h = 0.2; // 
int nx = R*2*pi/h;
int ny = r*2*pi/h;
func torex= (R+r*cos(y*pi*2))*cos(x*pi*2);
func torey= (R+r*cos(y*pi*2))*sin(x*pi*2);
func torez= r*sin(y*pi*2);


meshS ThS=square3(nx,ny,[torex,torey,torez]) ;
mesh3 Th3=tetg(ThS,switch="paAAQYY"); //,nbofregions=1,regionlist=domain);

webplot(ThS);

fespace VhS(ThS,P1);
VhS u,v;
macro grad3(u) [dx(u),dy(u),dz(u)]  // EOM

solve Lap(u,v) = int2d(ThS)( u*v+grad3(u)'*grad3(v)) -int2d(ThS)((x+y)*v);
webplot(u,ThS);
server();