add missing files

Author: frederichecht

2023/pdf/FH-ffday23.zip

@@ -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);

2023/pdf/FH-ffday23/edp2-full/3d-leman.edp

@@ -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);

2023/pdf/FH-ffday23/edp2-full/3daximesh.edp

@@ -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); 

2023/pdf/FH-ffday23/edp2-full/Adapt-uv.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);
+ 
+} ;