add new file

Author: frederichecht

2022/pdf/S-CD.pdf

@@ -0,0 +1,171 @@
+// calcul du couplage entre deux domaines de section rectangulaire Ex*Ey
+// a une distance dist
+// en contact avec une résistance de contact 1/c
+// 
+// La face opposee de chaque rectangle est a temperature imposee T1 resp T2
+// Le couplage est intégré directement dans le système d'équation
+//
+// version avec les espaces composite 
+
+real x0=0,x2=4;
+real y0=0,y1=1;
+real L= x2-x0; 
+real k1=1;
+real k2=1;
+
+real c=1;
+
+real T0=100.0, T1=00.0;
+
+
+//  mesh construction.  
+
+//  the curve in y .
+ func  F = 2+0.5*sin(2*pi*(y)); 
+ func  dF=  -1*0.5*pi*2*cos(2*pi*(y));  
+ func  dux =  (T1-T0)/(L+ (k1/c)/sqrt(1.+0*square(dF))); // 
+ 
+ //  k du/dn = c[u]
+//   if we neglect the du/dy 
+//  the solution 
+//  du/dx is constant in x and
+//  du/dx*L + [u] = T1-T0 
+//  du/dx*L + N.x k/c  du/dx = T1-T0
+//  => 
+//  du/dx = (T1-T0) /(L+nx*k/c)
+
+ //  t ->( F,t )  t' = ( dF, 1)  n = (-1, dF) / sqrt( 1+ dF^2) 
+ // N.x = 1 /  sqrt( 1+ dF^2)
+
+ func  Te1= T0+x*dux;
+ func  Te2= T1+(x-L)*dux;
+ 
+real xb1=  F(0,0);
+real xu1=  F(0,1);
+
+border left(t=y1,y0) {x=x0; y=t;label=3;} 
+border right(t=y0,y1) {x=x2; y=t;label=4;}
+border bot1(t=x0,xb1) {y=y0; x=t;label=2;}
+border bot2(t=xb1,x2) {y=y0; x=t;label=2;}
+border up1(t=x2,xu1) {y=y1; x=t;label=2;}
+border up2(t=xu1,x0) {y=y1; x=t;label=2;}
+border curve(t=0,1) {y=t;x=F(0,t) ; label=10;}
+
+
+// Maillages coincidents
+int n=30;
+int ny=n*(y1-y0);
+func Bord = left(ny)+right(ny)+curve(ny*2)
+                  +bot1(n*(xb1-x0))+bot2(n*(x2-xb1))
+                  +up1(n*(x2-xu1))+up2(n*(x2-xu1));
+plot(Bord,wait=1);                  
+mesh Th=  buildmesh( Bord );
+// get region numbion to buid 2 meshes:
+int rg1=Th(x0+0.1,y0+0.1).region;
+int rg2=Th(x2-0.1,y1-0.1).region;
+//  2 meshes conform. 
+mesh Th1=trunc(Th,region==rg1);
+mesh Th2=trunc(Th,region==rg2);
+// now the 2 meshes are not conform. 
+//Th2 = adaptmesh(Th2,0.05,IsMetric=1,nbvx=1000000);
+//Th1 = adaptmesh(Th1,0.03,IsMetric=1,nbvx=1000000);
+
+
+plot(Th1,Th2,wait=1);
+
+fespace Vh1(Th1,P1); //Vh1 u1=T0, v1,Vh1x=x,Vh1y=y,Vh1lab=label;
+fespace Vh2(Th2,P1);
+
+//  two good approximation of the exact solution. 
+Vh1 ue1=Te1;
+Vh2 ue2=Te2;
+plot(ue1,ue2,wait=1,cmm="analytic approximation");
+
+varf V11(u,v) = int2d(Th1)(k1*dx(u)*dx(v) + k1*dy(u)*dy(v)) + int1d(Th1,10,mortar=1)(c*u*v)  + on(3,u=T0); 
+varf V12(u,v) = int1d(Th1,10,mortar=1)(-c*u*v); 
+varf V21(u,v) = int1d(Th2,10,mortar=1)(-c*u*v);
+varf V22(u,v) = int2d(Th2)(k2*dx(u)*dx(v) + k2*dy(u)*dy(v)) + int1d(Th2,10,mortar=1)(c*u*v)   + on(4,u=T1);
+
+matrix A11 = V11(Vh1,Vh1);
+matrix A12 = V12(Vh2,Vh1);
+matrix A21 = V21(Vh1,Vh2);
+matrix A22 = V22(Vh2,Vh2);
+//  test the integrale of on curve line.
+
+{ //  verification of matrix A12 and A21... 
+
+macro verif1(A,vv) {
+	real[int] b(A.m);
+	b=1.;
+	vv[]=A*b;
+	cout << " sum = " << vv[].sum << "== " << int1d(Th1,10)(-c) << " == " << int1d(Th2,10)(-c) <<endl;
+	assert ( abs(vv[].sum-int1d(Th1,10)(-c) ) < 1e-3);
+}//	
+
+Vh1 v1;
+Vh2 v2;
+verif1(A12,v1);
+verif1(A21,v2);
+} //  end of verification ..
+
+//matrix A = [[A11, A12],[A12', A22]];  // valable dans le cas conforme uniquement (autrement A12' different de A21)
+matrix A = [[A11, A12],[A21, A22]];     
+
+//  rhs compuation 
+real[int]  B1= V11(0,Vh1);
+real[int] B2 = V22(0,Vh2);
+real[int] B = [B1,B2];
+
+Vh1 u1;
+Vh2 u2;
+set(A, solver=UMFPACK);
+real[int] T = A^-1*B;
+[u1[],u2[]]=T;  // dispatch the solution. 
+
+plot(u1,u2,value=1,fill=1,wait=1,nbiso=100);
+plot(u1,ue1,value=1,fill=0,wait=1);
+plot(u2,ue2,value=1,fill=0,wait=1);
+cout << " T1=" << int1d(Th1,10)(u1)/int1d(Th1,10)(1.0) << " T2=" << int1d(Th2,10)(u2)/int1d(Th2,10)(1.0) << endl;
+
+// composite version 
+fespace cVh(<Vh1,Vh2>);
+
+varf compositeVall([u1,u2],[v1,v2]) = 
+  int2d(Th1)(k1*dx(u1)*dx(v1) + k1*dy(u1)*dy(v1)) 
++ int2d(Th2)(k2*dx(u2)*dx(v2) + k2*dy(u2)*dy(v2)) 
++ int1d(Th1,10,mortar=1)(c*(u1-u2)*(v1-v2) )
++ on(3,u1=T0)
++ on(4,u2=T1);
+//+ int1d(Th1,10,mortar=1)(-c*u2*v1) 
+//+ int1d(Th2,10,mortar=1)(-c*u1*v2)
+//+ int1d(Th2,10,mortar=1)(c*u2*v2) 
+
+real[int] cB = compositeVall(0,cVh);
+matrix    cA = compositeVall(cVh,cVh);
+set(cA, solver=UMFPACK);
+
+real[int] cT = cA^-1*cB;
+
+Vh1 uc1;
+Vh2 uc2;
+[uc1[],uc2[]]=cT;  // dispatch the solution. 
+
+Vh1 u1diff;
+u1diff[] = uc1[]-u1[];
+Vh2 u2diff;
+u2diff[] = uc2[]-u2[];
+
+plot(uc1,uc2,value=1,fill=1,wait=1,nbiso=100);
+plot(uc1,ue1,value=1,fill=0,wait=1);
+plot(uc2,ue2,value=1,fill=0,wait=1);
+plot(u1diff,u2diff,value=1,fill=1,wait=1,nbiso=100);
+cout << " T1=" << int1d(Th1,10)(u1)/int1d(Th1,10)(1.0) << " T2=" << int1d(Th2,10)(u2)/int1d(Th2,10)(1.0) << endl;
+
+matrix D=A -cA;
+cout << "D.linfty=" << D.linfty << endl;
+real[int] db(B.n);
+db=B-cB;
+cout << "db.linfty=" << db.linfty << endl;
+end;
+
+/// FIN DU PROGRAMME

2022/pdf/S-GS.pdf

@@ -0,0 +1,48 @@
+load "msh3"
+load "medit"
+int nn=10;
+// a potentiel flow on a lac .. 
+//  a first freefem++ 3d example 
+// ------  not to bad ......
+verbosity=3;
+mesh Th2("lac-leman-v4.msh");
+fespace Vh2(Th2,P1);
+Vh2 deep;
+{  Vh2 v; 
+	macro Grad(u) [dx(u),dy(u)] //
+	solve P(deep,v)= int2d(Th2)(Grad(deep)'*Grad(v))+int2d(Th2)(v)
+	+on(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,deep=-1);
+	deep = deep*5/abs(deep[].min);
+	plot(deep,wait=1,cmm="min: "+deep[].min+ " max: "+deep[].max );
+}
+Vh2 ux,uz,p2;
+int[int] rup=[0,200],  rdown=[0,100],rmid(17*2);
+for(int i=0;i<rmid.n;++i)
+  rmid[i]=1+i/2;
+cout << rmid << endl;
+real maxdeep = deep[].min;
+mesh3 Th=buildlayers(Th2,nn,
+  coef= deep/maxdeep,
+  zbound=[deep,0],
+  labelmid=rmid, 
+  labelup = rup,
+  labeldown = rdown);
+//medit("Leman",Th);
+fespace Vh(Th,P13d);
+Vh p,q;
+//  (-deep[].min)*c = 0.5 Km
+//  c = 
+real cc=(0.5/-deep[].min);
+cout << cc << " cc = " << endl;
+cc=1; //  otherwise bug in bounding condition ... 
+macro Grad(u) [dx(u),dy(u),cc*dz(u)] //
+
+real ain=int2d(Th,1)(1.);
+real aout=int2d(Th,2)(1.);
+cout << " area " << ain << " " << aout << endl;
+real din=1./ain;
+real dout=-1./aout;
+solve P(p,q)= int3d(Th)(Grad(p)'*Grad(q)+1e-5*p*q)-int2d(Th,1)(q*din)+int2d(Th,2)(q*dout);
+
+plot(p,wait=1,nbiso=30,value=1);
+medit("potentiel",Th,p);

2022/pdf/S-MB.pdf

@@ -0,0 +1,27 @@
+load "msh3" load "medit"
+mesh3 Th;
+try {Th=readmesh3("/tmp/3daxi.mesh");} // try to read 
+catch(...) { // if not build ...
+ func f=2*((0.1+(((x/3))*(x-1)*(x-1)/1+x/100))^(1/3.)-(0.1)^(1/3.));
+ real yf=f(1.2,0); 
+ border up(t=1.2,0.){ x=t;y=f;label=0;}
+ border axe2(t=0.2,1.15) { x=t;y=0;label=0;}
+ border hole(t=pi,0) { x= 0.15 + 0.05*cos(t);y= 0.05*sin(t); label=1;}
+ border axe1(t=0,0.1) { x=t;y=0;label=0;}
+ border queue(t=0,1) { x= 1.15 + 0.05*t; y = yf*t; label =0;}
+ int np= 100;
+ func bord= up(np)+axe1(np/10)+hole(np/10)+axe2(8*np/10)+ queue(np/10);
+ plot( bord); 
+ mesh Th2=buildmesh(bord);
+ plot(Th2,wait=1);
+ int[int] l23=[0,0,1,1]; 
+ Th=buildlayers(Th2,coef= max(.15,y/max(f,0.05)), 50 ,zbound=[0,2*pi]
+   ,transfo=[x,y*cos(z),y*sin(z)],facemerge=1,labelmid=l23);
+ savemesh(Th,"/tmp/3daxi.mesh");   
+}
+
+macro Grad(u) [dx(u),dy(u),dz(u)] //
+fespace Vh(Th,P1);  Vh u,v;
+solve Poisson(u,v) = int3d(Th)( Grad(u)'*Grad(v) ) - int3d(Th)( v) + on(1,u=1);
+plot(u,wait=1,nbiso=20,value=1);
+medit("u",Th,u);

2022/pdf/S-PV.pdf

@@ -0,0 +1,34 @@
+border C(t=0,2*pi){x=cos(t);y=sin(t);};
+
+mesh Th= buildmesh(C(20));
+func  u = ( (10*x^3+y^3) +  tanh(500*(sin(5*y)-2*x) )); 
+func  v = ( (10*y^3+x^3) +  tanh(5000*(sin(5*x)-2*y)  ));
+
+real xz = 0.5, yz = 0.3 , dz=0.05;
+func bb=[[xz-dz,yz-dz],[xz+dz,yz+dz]];
+int nn = Th.nv;
+ plot(Th, cmm = 0 + " nv =" + Th.nv, ps="../plots/aTh-"+0+".eps" );
+ 
+for(int iter=1;iter < 20; ++iter )
+{
+  Th = adaptmesh(Th,[u,v],err=0.06,nbvx=50000,thetamax=30);
+  if(abs(Th.nv - nn) < nn*0.005) break; 
+  nn = Th.nv;
+  plot(Th, cmm = iter + " nv =" + Th.nv, ps="../plots/aTh-"+iter+".eps" );
+  plot(Th, cmm = iter + " nv =" + Th.nv, ps="../plots/zTh-"+iter+".eps",   bb = bb,wait=1 ,WindowIndex=1);
+
+}
+//   plot ... 
+fespace Vh(Th,P1);
+Vh uh=u,vh=v;
+
+load "msh3" 
+meshS Thu= movemesh23(Th,transfo=[x,y,u/10]);Thu=change(Thu,fregion=1);
+meshS Thv= movemesh23(Th,transfo=[x,y,v/10]);Thv=change(Thv,fregion=2);
+
+
+plot(uh,dim=3, cmm="uh",wait=1,fill=1);
+plot(vh,dim=3, cmm="vh",wait=1,fill=1);
+plot(Thu,Thv,wait=1);
+//plot(Thv,wait=1);
+load "medit" medit("Thu",Thu,Thv); 

2022/pdf/edp-c/composite_contact_thermique.edp

@@ -0,0 +1,25 @@
+load "MetricPk" load "Element_P3"
+int[int] lab=[1,1,1,1];
+mesh Th = square(10,10,label=lab);
+Th=trunc(Th,x<0.2 | y<0.2, label=1);
+fespace Vh(Th,P3);
+fespace Metric(Th,[P1,P1,P1]);	
+Vh u,v;
+real error=0.01;
+solve  Poisson(u,v,solver=CG,eps=1.0e-6) =
+    int2d(Th,qforder=5)( u*v*1.0e-10+  dx(u)*dx(v) + dy(u)*dy(v)) 
+  + int2d(Th,qforder=5)( (x-y)*v);
+// Do sub optimal adaptation in Norm $W^{r,p}$ for $Pk$   
+real  Nt= 1000,r = 1, k =3, p=1; 
+ 
+for (int i=0;i< 4;i++)
+{   
+   Metric [m11,m12,m22];	
+   m11[]=MetricPk(Th,u ,   kDeg=k,rDeg=r,pExp=p, mass=Nt/2);
+   Th = adaptmesh(Th,m11,m12,m22,IsMetric=true, nbvx= Nt*2);
+   u=u;// reinterpolation of u on new mesh Th. 
+   Poisson; 
+   plot(Th,wait=1);
+   plot(u,wait=1);
+ 
+} ;