/*-- 
Chauffage par induction d'une billette
3A ENERGIE du module Transfert d'énergie par induction.
G. Vinsard
*/

//-- Coefficients et dimensions physiques
real ConductiviteBillette=1/(3*1.0e-6);            // S/m
real PermeabiliteBillette=10.0;                      // relative
real RayonBillette=0.1;                              // m
real HauteurBillette=0.4;                            // m
real RayonSpires=0.12;                               // m
real RayonSectionSpires=0.01;                        // m
real EspaceEntreSpires=0.01;                         // m 
real CourantInducteur=1000;                          // A
int  NombreSpires=10;                                //
real frequence=1000;                                 // Hz


//-- Parametres numeriques
real Lb= HauteurBillette/2.0;
real e=EspaceEntreSpires;
real a=RayonSectionSpires;
real Li=4*(2*a+e);
real LL=max(Li,Lb);
real Ri=RayonSpires;
real RayonExterieur=2*sqrt(Ri*Ri+LL*LL);             // m    
real DensiteMaillage=500;                             // m^-1

//-- Constantes
real mu0=4*3.141592653589793*1.0e-7;                 // H/m permeabilite du vide
real ConductiviteCuivre=1/(2.0*1.0e-8);              // S/m

//-- Abreviations
real sigma=ConductiviteBillette;
real mur=PermeabiliteBillette;
real R=RayonBillette;
real Rs=RayonSpires;
real f=frequence;
real I=CourantInducteur;
real Ns=NombreSpires;
real Rx=RayonExterieur;
real dens=DensiteMaillage;

//-- Coefficients intermédiaires
real omega=2*pi*f;
complex cf=1i*sigma*omega*mu0;
complex czero=0;
complex cjs=I*mu0/(pi*a^2);
real delta=sqrt(2.0/(omega*sigma*mu0*mur));
cout << endl << "delta = " << delta << " ; R = " << R << endl;


//-- Maillage
border b1(t=0,1){x=R*t;y=0;label=1;};         int l1=max(dens*R,2.0);
border b2(t=0,1){x=R+(Ri-R)*t;y=0;label=1;};  int l2=max(dens*(Ri-R),2.0);
border b3(t=0,1){x=Ri+(Rx-Ri)*t;y=0;label=1;};int l3=max(dens*(Rx-Ri),2.0);
border b5(t=0,1){x=0;y=t*Lb;label=2;};        int l5=max(dens*Lb,2.0);
border b6(t=0,1){x=0;y=Lb+(Rx-Lb)*t;label=2;};int l6=max(dens*(Rx-Lb),2.0);
border b8(t=0,1){x=R;y=Lb*t;label=3;};        int l8=l5;
border b9(t=0,1){x=R*t;y=Lb;label=3;};        int l9=l1;
border b21(t=0,2*pi){x=Rs+a*cos(t);y=e/2+a+0*(e+2*a)+a*sin(t);label=6;};int l21=max(dens*2*pi*a,3.0);
border b22(t=0,2*pi){x=Rs+a*cos(t);y=e/2+a+1*(e+2*a)+a*sin(t);label=7;};int l22=max(dens*2*pi*a,3.0);
border b23(t=0,2*pi){x=Rs+a*cos(t);y=e/2+a+2*(e+2*a)+a*sin(t);label=8;};int l23=max(dens*2*pi*a,3.0);
border b24(t=0,2*pi){x=Rs+a*cos(t);y=e/2+a+3*(e+2*a)+a*sin(t);label=9;};int l24=max(dens*2*pi*a,3.0);
border b11(t=0,pi/2){x=Rx*cos(t);y=Rx*sin(t);label=5;};int l11=max(dens*pi/2*Rx,2.0);
//plot(b1(l1)+b2(l2)+b3(l3)+b5(l5)+b6(l6)+b8(l8)+b9(l9)+b11(l11));
mesh m=buildmesh(b1(l1)+b2(l2)+b3(l3)+b5(-l5)+b6(-l6)+b8(l8)+b9(-l9)+b11(l11)
+b21(l21)+b22(l22)+b23(l23)+b24(l24));
plot(m);

//-- Indicateurs
fespace mP0(m,P0);mP0 reg=region; // Rituel ad hoc pour définir la fonction « reg » qui renvoie le label du sous-domaine dans lequel est le point dont les coordonnées sont ses arguments
int IndicateurBillette=reg(R/2.0,Lb/2.0);
int IndicateurInducteur=4;
int IndicateurAir=reg(R/2.0,(Lb+Rx)/2.0);
int Ib=IndicateurBillette;
int Ii1=reg(Rs,e/2+a+0*(e+2*a));
int Ii2=reg(Rs,e/2+a+1*(e+2*a));
int Ii3=reg(Rs,e/2+a+2*(e+2*a));
int Ii4=reg(Rs,e/2+a+3*(e+2*a));
int Ia=IndicateurAir;
cout << "Indicateurs : " << endl;
cout << "Billette = " << Ib << " ; inducteur = " << Ii1 << ", " << Ii2 << ", " << Ii3 << ", " << Ii4   << endl;
cout << "air =      " << Ia << endl;


//-- Espaces fonctionnels
fespace mP1(m,P1);
mP1<complex> A,Ap;

//-- Calcul 
solve Calc(A,Ap)=
  int2d(m,Ia,Ii1,Ii2,Ii3,Ii4)((x*dy(A)*dy(Ap)+x*dx(A)*dx(Ap)+dx(A)*Ap+dx(Ap)*A+1/x*Ap*A))
 +int2d(m,Ib)((x*dy(A)*dy(Ap)+x*dx(A)*dx(Ap)+dx(A)*Ap+dx(Ap)*A+1/x*Ap*A)/mur)
 +int2d(m,Ib)(x*A*Ap*cf)
 +int2d(m,Ii1)(x*Ap*cjs)
 +int2d(m,Ii2)(x*Ap*cjs)
 +int2d(m,Ii3)(x*Ap*cjs)
 +int2d(m,Ii4)(x*Ap*cjs)
 +on(2,5,A=czero);

//-- Exploitation
mP1 Ar=real(A), Ai=imag(A), rAr=x*Ar, rAi=x*Ai;
mP1<complex>  Ac=conj(A);
plot(rAr,cmm="r Ar");
plot(rAi,cmm="r Ai");

complex ctest=int2d(m,Ii1)(A);
cout << ctest << endl;

complex cS=int2d(m,Ii1,Ii2,Ii3,Ii4)(x*A)*cjs*2*pi/mu0*1i*omega*2;
complex cPj=int2d(m,Ib)(x*A*Ac)*2*pi*sigma*omega*omega*2;
complex cQ=
  int2d(m,Ib)((x*dy(A)*dy(Ac)+x*dx(A)*dx(Ac)+dx(A)*Ac+dx(Ac)*A+1/x*Ac*A)/mur)
+ int2d(m,Ia,Ii1,Ii2,Ii3,Ii4)((x*dy(A)*dy(Ac)+x*dx(A)*dx(Ac)+dx(A)*Ac+dx(Ac)*A+1/x*Ac*A));
cQ=1i*omega*cQ/mu0*2*pi*2;

cout << "S = " << real(cS) << " + " << imag(cS) << " i "  << endl ;
cout << "Pj = " << real(cPj)  << " ; Q = " << imag(cQ) << endl ;



cout << endl << "ok" << endl;