load "msh3"; load "Element_P1dc1"; real k= 1e-6; verbosity = 0; int itermax=2; real T=1.; real pena=10.; for ( int i=1;i<=itermax;i++) { real cpu = clock(); int NN = 2^i; mesh3 Th = cube(NN,NN,NN,flags=2); int[int] labs=labels(Th); real dt=1./(NN^1); real TN=T/dt; fespace Vh(Th,P1dc); Vh ux, uy,uz, uu1, uu2,uu3, dux,duy,duz, vx, vy, vz; uu1=0; uu2=0; uu3=0; for (int j=1; j<=TN; j++) { real t=j*dt; func U1b= 0; func U2b= 0; func U3b= 0; func f1 = x; func f2 = y ; func f3 = z; int n; real err=0; // update solution uu1[]=ux[]; uu2[]=uy[]; uu3[]=uz[]; solve SAMPLE([dux, duy, duz], [vx, vy, vz], solver=sparsesolver)= int3d(Th)(((1./dt)*k)*(dx(dux)*dx(vx) + dy(dux)*dy(vx) + dz(dux)*dz(vx) + dx(duy)*dx(vy) + dy(duy)*dy(vy) + dz(duy)*dz(vy) + dx(duz)*dx(vz) + dy(duz)*dy(vz) + dz(duz)*dz(vz)) ) + intallfaces(Th)( (((1./dt)*k)*(mean(dx(dux))*N.x*jump(vx)+mean(dy(dux))*N.y*jump(vx)+mean(dz(dux))*N.z*jump(vx)+mean(dx(duy))*N.x*jump(vy)+mean(dy(duy))*N.y*jump(vy)+mean(dz(duy))*N.z*jump(vy)+mean(dx(duz))*N.x*jump(vz)+mean(dy(duz))*N.y*jump(vz)+mean(dz(duz))*N.z*jump(vz) +mean(dx(vx))*N.x*jump(dux)+mean(dy(vx))*N.y*jump(dux)+mean(dz(vx))*N.z*jump(dux)+mean(dx(vy))*N.x*jump(duy)+mean(dy(vy))*N.y*jump(duy)+mean(dz(vy))*N.z*jump(duy)+mean(dx(vz))*N.x*jump(duz)+mean(dy(vz))*N.y*jump(duz)+mean(dz(vz))*N.z*jump(duz) + (pena*NN)*(jump(dux)*jump(vx)+jump(duy)*jump(vy)+jump(duz)*jump(vz)) ) )/ nElementonB ) -int3d(Th)( ((1./dt)*k) * (dx(uu1)*dx(vx) + dy(uu1)*dy(vx) + dz(uu1)*dz(vx) + dx(uu2)*dx(vy) + dy(uu2)*dy(vy) + dz(uu2)*dz(vy) + dx(uu3)*dx(vz) + dy(uu3)*dy(vz) + dz(uu3)*dz(vz)) ) -intallfaces(Th)( ((1./dt)*k *(mean(dx(uu1))*N.x*jump(vx)+mean(dy(uu1))*N.y*jump(vx)+mean(dz(uu1))*N.z*jump(vx)+mean(dx(uu2))*N.x*jump(vy)+mean(dy(uu2))*N.y*jump(vy)+mean(dz(uu2))*N.z*jump(vy)+mean(dx(uu3))*N.x*jump(vz)+mean(dy(uu3))*N.y*jump(vz)+mean(dz(uu3))*N.z*jump(vz) +mean(dx(vx))*N.x*jump(uu1)+mean(dy(vx))*N.y*jump(uu1)+mean(dz(vx))*N.z*jump(uu1)+mean(dx(vy))*N.x*jump(uu2)+mean(dy(vy))*N.y*jump(uu2)+mean(dz(vy))*N.z*jump(uu2)+mean(dx(vz))*N.x*jump(uu3)+mean(dy(vz))*N.y*jump(uu3)+mean(dz(vz))*N.z*jump(uu3) + (pena*NN)*(jump(uu1)*jump(vx)+jump(uu2)*jump(vy)+jump(uu3)*jump(vz)) ))/ nElementonB ) -int3d(Th)([f1,f2,f3]'*[vx,vy,vz]) + on(labs, dux=U1b, duy=U2b, duz=U3b) ; } // Plot plot( cmm="[ux] ", ux ); }