load "PETSc" macro dimension()2// EOM include "macro_ddm.idp" load "iovtk" load "Curvature" load "isoline" bool importMesh = false; string meshFileName = "mesh/bend_pipe_2d.mesh"; string outputFileName = "output/"; int NN=1000, saveEach = 10; real muf=10, c1=1.e7, rhof=100, rhos=100, g=-0, rhod=rhos-rhof; real T=1, t, dt=T/NN, penal=0; int[int] orderOut = [1, 1, 1],orderOuts = [1, 1]; func PV1=[P2,P2,P1]; func PV2=[P2,P2]; mesh Th, Ths, Ths0; if (importMesh) { Th = readmesh(meshFileName); } else { Th = square(160,40, [4*x,y]); //plot(Th, wait=1); Ths = square(1,40, [0.0212*x+2-0.0106,0.8*y]); //plot(Ths, Th, wait=1); } cout << "Number of Elements on Rank before buildDmesh(Th)"<< mpirank<<": "<< Th.nt << endl; Mat A,B,P; buildDmesh(Th); { macro def(i)[i, i#B, i#C]// macro init(i)[i, i, i]// createMat(Th, A, PV1) } buildDmesh(Ths); { macro def(i)[i, i#B, i#C]// macro init(i)[i, i, i]// createMat(Ths, B, PV1) } Ths0=Ths; cout << "Number of Elements on Rank after buildDmesh(Th)"<< mpirank<<": "<< Th.nt << endl; /* The fespace has to be defined after real mesh is buit */ fespace Rh(Th, PV1); fespace Rhs(Ths, PV1); fespace Vh (Th,P2); Vh uu,ux,uy,uhx,uhy,uoldx=0,uoldy=0; fespace Ph (Th,P1); Ph p,ph; fespace Vhs(Ths, P2); Vhs wx,wy,whx,why,woldx=0,woldy=0; Vhs s11=0,s12=0,s22=0; fespace Phs (Ths,P1); Phs q,qh; fespace Vhs0(Ths0, P2); Vhs0 w0x,w0y,dispx=0,dispy=0; Vhs0 s011,s012,s022,f11,f12,f21,f22; macro div(u) (dx(u#x) + dy(u#y))// macro grad(u) [dx(u), dy(u)]// macro Grad(u) [grad(u#x), grad(u#y)]// macro DD(u) [[2*dx(u#x),(dx(u#y) + dy(u#x))],[(dx(u#y) + dy(u#x)),2*dy(u#y)]]// macro S[[s11,s12],[s12,s22]] // EOM varf fluid([ux,uy,p],[uhx,uhy,ph]) = int2d(Th)( rhof*[ux,uy]'*[uhx,uhy]/dt-div(uh)*p-div(u)*ph + penal*p*ph + muf/2*trace(DD(u)'*DD(uh))) + int2d(Th)(g*rhof*uhy + rhof*[convect([uoldx,uoldy],-dt,uoldx), convect([uoldx,uoldy],-dt,uoldy)]'*[uhx,uhy]/dt ) + on(1,ux=0,uy=0) + on(3,uy=0) + on(4,ux=15*y*(2-y)*sin(2*pi*t),uy=0) + on(2, p=0); varf solid([wx,wy,q],[whx,why,qh]) = int2d(Ths)( rhod*[wx,wy]'*[whx,why]/dt + 0*q*qh + dt*c1*trace((DD(w)+ dt*Grad(w)*Grad(wold)'+ dt*Grad(wold)*Grad(w)')*Grad(wh)') + dt*trace((Grad(w)*S + S*Grad(w))*Grad(wh)') + dt*dt*trace((Grad(w)*S*Grad(wold)' + Grad(wold)*S*Grad(w)')*Grad(wh)')) + int2d(Ths)(g*rhod*why + trace((c1*dt*dt*Grad(wold)*Grad(wold)' + dt*dt*Grad(wold)*S*Grad(wold)'-S)*Grad(wh)') + rhod*[woldx, woldy]'*[whx,why]/dt ); set(A, sparams = "-pc_type lu"); set(B, sparams = "-pc_type lu"); for(int n = 1; n <= NN; n++) { if(mpirank == 0) cout << "Time Step: " << n << endl; t=n*dt; { macro def(i)[i, i#B, i#C]// macro init(i)[i, i, i]// createMat(Th, A, PV1) createMat(Ths, B, PV1) transferMat(Th, PV1, A, Ths, PV1, B, P) A = fluid(Rh, Rh); B = solid(Rhs,Rhs); real[int] rhs1 = fluid(0,Rh); real[int] rhs2 = solid(0,Rhs); real[int] rhs(Rh.ndof); rhs = P' * rhs2; rhs += rhs1; Mat C; MatPtAP(B, P, C); A += C; Rh [utx,uty,pt]; real[int] sol(Rh.ndof); sol= utx[]; sol = A^-1 * rhs; utx[]=sol; ux=utx; uy=uty; p=pt; } wx=ux; wy=uy; w0x=0; w0y=0; w0x[] = wx[]; w0y[] = wy[]; dispx[] += w0x[]*dt; dispy[] += w0y[]*dt; f11=dx(dispx)+1; f12=dy(dispx); f21=dx(dispy); f22=dy(dispy)+1; //S0=c1*(F*F'-II); s011=c1*(f11*f11+f12*f12-1); s012=c1*(f11*f21+f12*f22); s022=c1*(f21*f21+f22*f22-1); s11=0;s12=0;s22=0; s11[]=s011[]; s12[]=s012[]; s22[]=s022[]; Ths = movemesh(Ths0, [x+dispx,y+dispy]); wx=wx; wy=wy; q=q; s11=s11;s12=s12;s22=s22; uoldx=ux; uoldy=uy; woldx=wx; woldy=wy; uu=sqrt(ux^2+uy^2); plotD(Th,uu,wait=0); //plotD(Ths,wx,wait=0); if (n % saveEach == 0){ savevtk(outputFileName+"fluid.vtu", Th, ux, uy, p, dataname="ux uy p", order=orderOut,append=true); savevtk(outputFileName+"solid.vtu", Ths, dispx, dispy, dataname="dx dy", order=orderOuts,append=true); } }