load "PETSc" macro dimension()2// EOM // 2D or 3D //macro vectorialfe()P1// EOM include "macro_ddm.idp" // additional DDM functions macro def(i)[i, i#B]// EOM // vector field definition macro init(i)[i, i]// EOM // vector field initialization func Pk = [P1, P1]; // finite element space //modeling //background //mesh Th = square(200, 200); border Top(t=0,1){x = t; y = 1;label=1;} border Right(t=0,1){x = 1; y = 1-t;label=2;} border Bottom(t=0,1){x = 1-t; y = 0;label=3;} border Left(t=0,1){x = 0; y = t;label=4;} //source_kong border Top1(t=0,1){x = 0.45+0.1*t; y = 0.99;label=5;} border Right1(t=0,1){x = 0.55; y = 0.99-t/20;label=6;} border Bottom1(t=0,1){x = 0.55-0.1*t; y = 0.94;label=7;} border Left1(t=0,1){x = 0.45; y = 0.94+t/20;label=8;} //Fracture border Top2(t=0,1){x = 0.2+0.6*t; y = 0.55;label=9;} border Right2(t=0,1){x = 0.8; y = 0.55-0.1*t;label=10;} border Bottom2(t=0,1){x = 0.8-0.6*t; y = 0.45;label=11;} border Left2(t=0,1){x = 0.2; y = 0.45+0.1*t;label=12;} int n = 100, n1=10; //build mesh mesh solid = buildmesh(Top(-n)+Left(-n)+Right(-n)+Bottom(-n)+Top1(n1)+Left1(n1/2)+Right1(n1/2)+Bottom1(n1) +Top2(6*n1)+Left2(n1)+Right2(n1)+Bottom2(6*n1)); mesh source1=buildmesh(Top1(-n1)+Left1(-n1)+Right1(-n1)+Bottom1(-n1)); mesh fracture=buildmesh(Top2(-6*n1)+Left2(-n1)+Right2(-n1)+Bottom2(-6*n1)); mesh Th=source1+solid+fracture;//+fracture //int[int] reg(2); //fespace Ph(Th,P0); //Ph reg=region; //int Background=reg(0.05,0.05); //int Fracture2=reg(0.3,0.51); //int Background = Th(0.5, 0.1).region, Fracture2 = Th(0.5, 0.5).region; //reg(0)=Th(0.5, 0.1).region; //reg(1)=Th(0.5, 0.5).region; //cout << " Background= " << reg(0) << endl; //cout << " Fracture2= " << reg(1) << endl; //plot(reg,fill=true,value=true); plot(Th,wait=1); //parameter of background real density=2650.0,//kg/m^3 lambda=6*10^9,//pa mu= 45.0*10^9;//pa //parameter of water real density2=1035.0,//kg/m^3 lambda2= 2.3*10^9,//Pa mu2=0, nu=1.11e-3;//Pa·S,dynamic viscosity,动力粘度-The engineering toolbx //Ph density=density1*(region==Background)+density2*(region==Fracture2); //Ph lambda=lambda1*(region==Background) + lambda2*(region==Fracture2); //Ph mu=mu1*(region==Background)+mu2*(region==Fracture2); //plot(density,fill=true,value=true); fespace Vh(Th,Pk); Vh [u,uB], [v,vB],[u1f,u2f]=[0,0],[u1b,u2b]=[0,0]; fespace Wh(Th,P1); Wh sigma11=0, sigma22=0, sigma12=0;//force real sqrt2=sqrt(2.); macro epsilon(u,uB) [dx(u),dy(uB),(dy(u)+dx(uB))/sqrt2] // EOM macro div(u,uB) ( dx(u)+dy(uB))// EOM real dt=0.000001,t=0,idt2=1./(dt*dt),amplitude=1;//0.001;0.000001 real f=80000., //Hz tt=30.*dt, //m idt=1./dt; int i=0; Mat A; createMat(Th, A, Pk) //solve varf Lame([u,uB],[v,vB])= int2d(Th)(lambda*div(u,uB)*div(v,vB)+2.*mu*( epsilon(u,uB)'*epsilon(v,vB))) +int2d(Th)(density*idt2*u*v)+int2d(Th)(density*2*idt2*u1f*v)-int2d(Th)(density*idt2*u1b*v)+ int2d(Th)(density*idt2*uB*vB)+int2d(Th)(density*2*idt2*u2f*vB)-int2d(Th)(density*idt2*u2b*vB) -int1d(Th,9,10,11,12)(-0.2*(sigma11*N.x*v + sigma22*N.y*vB + sigma12*(N.y*v+N.x*vB))) +on(5,6,7,8,u=0.,uB=amplitude*(1.-2.*(pi*f*(t-tt))^2)*exp(-(pi*f*(t-tt))^2)) +on(1,2,3,4,u=0.,uB=0.); //func Pm=[P1,P1]; macro def1(i)[i, i#B,i#C]// EOM // vector field definition macro init1(i)[i, i,i]// EOM // vector field initialization func Pm = [P2, P2,P1]; // finite element space fespace Mh(fracture, P1); //definition of the pressure space Mh up1=0.,up2=0.; fespace Xh(fracture, Pm); //definition of the velocity component space Xh [uu,uuB,p]=[0.,0.,0.],[vv,vvB,q]=[0.,0.,0.];// //def1(uu),def1(vv); Mat B; //buildMat(fracture, getARGV("-split", 1), B, Pk, mpiCommWorld) createMat(fracture,B,Pm) macro grad(uu)[dx(uu),dy(uu)]//EOM macro div1(uu)( dx(uu)+dy(uu#B) )//EOM varf NS ([uu,uuB,p],[vv,vvB,q])= int2d(fracture)(idt*(uu*vv+uuB*vvB)+(nu*grad(uu)'*grad(vv)+grad(uuB)'*grad(vvB)) -div1(uu)*q-div1(vv)*p + p*q*1e-10) //1e-6这个值太大可能会卡bug +int2d(fracture)(-idt*convect([up1,up2],-dt,up1)*vv-idt*convect([up1,up2],-dt,up2)*vvB) +on(9,10,11,12, uu=(u-u1f)*idt, uuB=(uB-u2f)*idt)//boundary condition ; [u1b,u2b]=[0,0]; [u1f,u2f]=[u1b,u2b]; [u,uB]=[u1f,u2f]; ofstream point3("uB(0.5,0.5).txt"); int iter=0; int numcolo=2;//5 real [int] viso(numcolo*10+1); for (int r=0;r<(numcolo*10+1);r++) {viso[r]=(-numcolo*0.5+0.1*r)*amplitude; } for(int i=0;i<150;++i) { t=t+dt; matrix Loc = Lame(Vh, Vh); //solve the problem real[int] rhs = Lame(0, Vh); Vh def(u); // local solution set(A, sparams = "-ksp_view -ksp_max_it 100 -pc_type ksp", bs=1);//, nearnullspace=u A = Loc; A=Lame(Vh,Vh); u[] = A^-1 * rhs; macro def2(u)u // EOM matrix Locns = NS(Xh, Xh); //solve the problem real[int] rhsns = NS(0, Xh); Xh def1(uu); // local solution [uu,uuB,p] set(B, sparams = "-ksp_view -ksp_max_it 100 -pc_type ksp", bs=1);//, nearnullspace=u B = Locns; uu[] = B^-1 * rhsns; macro def3(uu)uu // EOM sigma11([x, y]) = (2*dx(uu) - p); sigma22([x, y]) = (2*dy(uuB) - p); sigma12([x, y]) = (dx(uu) + dy(uuB)); up1=uu; up2=uuB; [u1b,u2b]= [u1f,u2f]; [u1f,u2f]=[u,uB]; point3 << uB(0.5,0.5)<< endl; macro params()cmm = "Global solution time="+t, fill = 1// EOM plotMPI(Th, uB, P1, def2, real, params) //def1(u)[P1,P1] // plot(u2,cmm="time="+t,fill=true,value=false,viso=viso(0:viso.n-1),coef=0.003); // int[int] Order=[1,1]; //savevtk("1.vtu", Th, [u1,u2], dataname="u1,u2", order=Order,append = true); }