load "PETSc" load "tetgen" include "cube.idp" macro dimension()3// EOM include "macro_ddm.idp" load "iovtk" load "Curvature" load "isoline" load "msh3" int nn = 10; real mu=1.e0; real u0x=1, u0y=0, u0z=0, u1x=0,u1y=0,u1z=0; mesh3 Th = cube(nn, nn, nn); mesh3 ThG=Th; func Pk=[P2,P2,P2,P1]; func Pm=[P2,P2,P2]; int[int] fn2o; macro ThN2O() fn2o // buildDmesh(Th); fespace Vh(Th, Pk); Vh [ux,uy,uz,p]; fespace VhG(ThG, Pk); VhG [uxG,uyG,uzG,pG],[uxR,uyR,uzR,pR]; int[int] findex = restrict(Vh, VhG, fn2o); mesh ThS=square(nn,nn); /* if(mpirank == 0) { // mesh3 Thg = buildBdMesh(ThG); //meshS ThS = Thg.Gamma; //ThS=trunc(ThSg,(x < 0.001),split=2); int[int] labs = [1]; ThS = extract(ThG, label=labs); //plot(ThS, wait=1); } broadcast(processor(0), ThS); */ fespace Sh(ThS, Pm); Sh [ax,ay,az]; macro Grad(u) [[dx(u#x), dy(u#x), dz(u#x)], [dx(u#y), dy(u#y), dz(u#y)], [dx(u#z), dy(u#z), dz(u#z)]] // macro div(u) (dx(u#x)+dy(u#y)+dz(u#z)) // macro vect(u) [u#x, u#y, u#z] // varf uu([ux,uy,uz,p],[uhx,uhy,uhz,ph]) = int3d(Th)(mu*(Grad(u):Grad(uh)) -div(uh)*p -div(u)*ph) +on(2,3,4,5,6, ux=0, uy=0, uz=0); varf ua([ux,uy,uz,p],[ahx,ahy,ahz]) = int2d(ThS)(vect(u)'*vect(ah)); varf aa([ax,ay,az],[ahx,ahy,ahz]) = int2d(ThS)(0.*vect(ah)'*vect(a)); varf af([ax,ay,az],[ahx,ahy,ahz]) = int2d(ThS)(vect(u0)'*vect(ah)); string[int] names(2); names[0] = "velocity"; names[1] = "pressure"; Mat UU; { macro def(i)[i, i#B, i#C, i#D]// macro init(i)[i, i, i, i]// createMat(Th, UU, Pk) } UU = uu(Vh, Vh, tgv=-1); matrix AA= aa(Sh,Sh); matrix UAtemp = ua(Vh,Sh); //Mat UA(UU,AA,UAtemp); /* matrix UU=uu(Vh, Vh, tgv=-1); matrix UA=ua(Vh, Sh); matrix AA=aa(Sh, Sh); */ /* real[int] Af = af(0,Sh); Mat M=[[UU, UA'], [UA, AA]]; int undof = Vh.ndof; int andof = Sh.ndof; int tndof = undof + andof; real[int] f(tndof); real[int] sol(tndof); f(0:undof - 1) = 0; f(undof:tndof - 1) = Af; set(M, solver=UMFPACK); sol = M^-1*f; ux[] = sol(0:undof - 1); ax[] = sol(undof:tndof - 1); plot(Th, ax, fill=1,value=1, nbiso=20); int[int] orderOut = [1, 1, 1, 1, 1, 1, 1]; savevtk("fluid.vtu", Th, ux,uy,uz,p, ax,ay,az, dataname="ux uy uz p ax ay az", order=orderOut); */