/*-- 
Aimant cylindrique en chute libre dans un tube en cuivre
1A ENSEM -- Apprentissage -- Base de l'électromagnétisme
G. Vinsard
*/

//-- Coefficients et dimensions physiques

real M=1.0e+6,Ha=0.03,              // densité d'aimantation (A/m) et hauteur de l'aimant (m)
     Re=0.02/2,Ri=0.01/2,           // rayons intérieur et extérieur du tube  (rayon aimant = Ri) (m)
     Rx=6*Re,Hx=Rx*2,               // rayon et auteur pour borner le domaine de calcul (m)
     sigma=58e+6,                   // conductivité tube (S/m)
     v=0.1,                         // vitesse de l'aimant (m/s)
     mu0=4.0e-7*pi,                 // perméabilité magnétique du vide (H/m)
     dens=1000.0;                   // densité de maillage (nb de triangles par m)


macro colind(ind,m){fespace mP0(m,P0);mP0 reg;real [int] tb(m.nt);for (int el=0;el<m.nt;el++){if (m[el].region==ind) {tb[el]=1.0;} else {tb[el]=0.0;};};reg[]=tb;plot(reg,fill=1);}//
macro minmax(phi,vphi) cout<<phi+" : min = "<<vphi[].min<<" ; "<<" max = "<<vphi[].max<<endl//
macro InfoM(m)cout << "nbno = " << m.nv << " ; " << " nbelv = " << m.nt;{int [int] inds(1000);int Ninds=-1;for (int el=0;el<m.nt;el++) {int ind=m[el].label;int ok=1;for (int i=0;i<=Ninds;i++) if (inds[i] == ind) ok=0;if (ok == 1) {Ninds=Ninds+1;inds[Ninds]=ind;};};cout << " ; nbindv= " << Ninds +1 << " , soit : ";for (int i=0;i<=Ninds;i++) cout << inds[i] << " " ;};{int [int] indl(1000);int Nindl=-1;for (int el=0;el<m.nbe;el++) {int ind=m.be(el).label;int ok=1;for (int i=0;i<=Nindl;i++) if (indl[i] == ind) ok=0;if (ok == 1) {Nindl=Nindl+1;indl[Nindl]=ind;};};cout << " ; nbindl= " << Nindl +1 << " , soit : ";for (int i=0;i<=Nindl;i++) cout << indl[i] << " " ;cout << endl;}//

macro maillage (Ha,Ri,Re,Rx,Hx)
 real Has2=Ha/2,Hxs2=Hx/2,Rxs2=Rx/2;
 border b1(t=0,1){x=Ri*t;y=-Hxs2;label=2;};int l1=Ri*dens;
 border b2(t=0,1){x=Ri+(Re-Ri)*t;y=-Hxs2;label=2;};int l2=(Re-Ri)*dens;
 border b3(t=0,1){x=Re+(Rx-Re)*t;y=-Hxs2;label=2;};int l3=(Rx-Re)*dens;
 border b4(t=0,1){x=Ri*t;y=Has2;label=1;};int l4=l1;
 border b5(t=0,1){x=Ri*t;y=Hxs2;label=2;};int l5=Ri*dens;
 border b6(t=0,1){x=Ri+(Re-Ri)*t;y=Hxs2;label=2;};int l6=(Re-Ri)*dens;
 border b7(t=0,1){x=Re+(Rx-Re)*t;y=Hxs2;label=2;};int l7=(Rx-Re)*dens;
 border b8(t=-1,1){x=0;y=Has2*t;label=4;};int l8=Ha*dens;
 border b9(t=0,1){x=0;y=Has2+(Hxs2-Has2)*t;label=4;};int l9=(Hxs2-Has2)*dens;
 border b10(t=-1,1){x=Ri;y=Has2*t;label=3;};int l10=2*l8;
 border b11(t=0,1){x=Ri;y=Has2+(Hxs2-Has2)*t;label=1;};int l11=l9;
 border b12(t=-1,1){x=Re;y=Hxs2*t;label=1;};int l12=Hxs2*dens;
 border b13(t=-1,1){x=Rx;y=Hxs2*t;label=2;};int l13=l12;
 border b14(t=0,1){x=Ri*t;y=-Has2;label=1;};int l14=l4;
 border b15(t=0,1){x=0;y=-Has2-(Hxs2-Has2)*t;label=4;};int l15=l9;
 border b16(t=0,1){x=Ri;y=-Has2-(Hxs2-Has2)*t;label=1;};int l16=l11;
 plot(b1(l1)+b2(l2)+b3(l3)+b4(l4)+b5(-l5)+b6(-l6)+b7(-l7)+b8(-l8)+b9(-l9)+b10(l10)+b11(l11)+b12(l12)+b13(l13)+b14(l14)+b15(-l15)+b16(l16));
 mesh m=buildmesh(b1(l1)+b2(l2)+b3(l3)+b4(l4)+b5(-l5)+b6(-l6)+b7(-l7)+b8(-l8)+b9(-l9)+b10(l10)+b11(l11)+b12(l12)+b13(l13)+b14(l14)+b15(l15)+b16(l16))//


maillage(Ha,Ri,Re,Rx,Hx);plot(m,cmm="maillage");

fespace P0(m,P0);P0 reg=region;
int Ia=reg(Ri/2,Ha/4),It=reg((Ri+Re)/2,Hx/4),Iai=reg(Ri/2,Ha/2+(Hx-Ha)/4),Iax=reg(Re+(Rx-Ri)/4,Ha/2+(Hx-Ha)/4);

fespace P1(m,P1); P1 a0,ra0,a1,ra1,ap,rda;

solve Pb0(a0,ap)=
 int2d(m)(x*dx(a0)*dx(ap)+dx(a0)*ap+a0*dx(ap)+a0*ap/x+x*dy(a0)*dy(ap))
+int2d(m,Ia)(x*dx(ap)+ap)
+on(2,4,a0=0);
ra0=x*a0;
plot(ra0,nbiso=11,bw=true,cmm="Lignes induction : aimant seul");


solve Pbf(a1,ap)=
 int2d(m)(x*dx(a1)*dx(ap)+dx(a1)*ap+a1*dx(ap)+a1*ap/x+x*dy(a1)*dy(ap))
+int2d(m,It)(sigma*mu0*v*x*dy(a1)*ap)
+int2d(m,Ia)(x*dx(ap)+ap)
+on(2,4,a1=0);
ra1=x*a1;
plot(ra1,nbiso=11,bw=true,cmm="Lignes induction : aimant + courants induits");
real Pj=int2d(m,It)(x*dy(a1)*dy(a1))*2*pi*v*v*sigma;

Pj=Pj*(M*mu0)^2;
cout <<  "Puissance Joule dans le tube = " << Pj << endl;

rda=ra1-ra0;
plot(rda,nbiso=11,bw=true,cmm="Lignes induction : courants induits seuls");