OPTI DIME 2 ELEM TRI3;
* example pour MPMA + JPMA
* calcul d'inductance d'un parallélépipède dans CAST3M
* et comparaison avec la formule analytique

* l'intérêt de MPMA + JPMA : on peut calculer l'énergie magnétique
* (et donc l'inductance) associée à une densité de courant J(x,y) quelconque

* IMPORTANT: J(x,y) est donné en 2D (avec une épaisseur selon Z)
* mais le potentiel magnétique est calculé en 3D



OPTI 'ECHO' 0;

* afficher ou pas
AFFICH = FAUX;

* longueur, largeur, épaisseur
LX = 100.E-3;
LY = 5.E-3;
EPAI = 12.E-3;



CT0 = POIN 0. 0.;

N1 = 6;

TAILLE1 = LY/N1;

* maillage
DENS TAILLE1;


* on prend la densité de courant pour que le courant soit égal à 1 A

JA = (1./LY)/EPAI;


X1 = (-1.)*LX/2;
Y1 = (-1.)*LY/2;

X2 = (-1.)*LX/2;
Y2 = LY/2;

X3 = LX/2;
Y3 = LY/2;

X4 = LX/2;
Y4 = (-1.)*LY/2;



P1 = POIN X1 Y1;
P2 = POIN X2 Y2;
P3 = POIN X3 Y3;
P4 = POIN X4 Y4;

LD1 = D P1 P2;
LD2 = D P2 P3;
LD3 = D P3 P4;
LD4 = D P4 P1;


CT = LD1 ET LD2 ET LD3 ET LD4;


SU1 = SURF CT;


ELIM SU1 1.E-7;

* il faut un modèle
MTB = 'MODE' SU1 'THERMIQUE' 'ISOTROPE';

* remarque: ici on fait le plus simple, on impose une densité de courant;
*           il est également possible de résoudre un problème "thermique", 
*           afin de trouver une densité de courant inconnue,
*           et ensuite on peut calculer le potentiel vecteur, l'énergie magnétique, etc.
* dans ce cas on le ferait de mainère suivante :
*SIGMA = 1;
*
*MAT1 = MATE MTB 'K' SIGMA;
*
*
*COND1 = COND MAT1 MTB;
*V0 = 0;
*V1 = 1;
*
*B0 = BLOQ 'T' LD3;
*
*CL0 = DEPI B0 V0;
*CL1 = FLUX MTB LD1 JA;
*
*V = RESO (COND1 ET B0) (CL0 ET CL1);
*
*J = GRAD V MTB;
*
*J = J * -1.;
*
*JX = EXCO 'T,X' J;
*JY = EXCO 'T,Y' J;
*
*JX = CHAN 'COMP' 'JX' JX;
*
*JY = CHAN 'COMP' 'JY' JY;
*
*J = JX ET JY;

*TRAC J MTB SU1 TITR 'J obtenu par RESO';

JXC = MANU 'CHPO'  SU1  'JX' JA;
JYC = MANU 'CHPO' SU1 'JY' 0.;


JX = CHAN 'CHAM' JXC SU1;

JY = CHAN 'CHAM' JYC SU1;



* il faut une composante 'JX' et une autre 'JY'
JX = CHAN 'COMP' 'JX' JX;

JY = CHAN 'COMP' 'JY' JY;



J = JX ET JY;


SI AFFICH;
TRAC J MTB SU1 TITR 'DENSITE DE COURANT';
FINSI;

* on calcule la martice M (une seule fois)
M = MPMA MTB EPAI;

* et ensuite on calcule A
* si J change, on peut recalculer A avec la même matrice M
A = JPMA J MTB M;


SI AFFICH;
TRAC A MTB SU1 TITR 'POT MAG VECT MOYENNE SUR L''EPAISSEUR';
FINSI;

* on a le potentiel vecteur, on peut calculer l'énergie magnétique

AX = EXCO 'AX' A;
AY = EXCO 'AY' A;

CJX = CHAN 'CHPO' MTB JX;

CAX = CHAN 'CHPO' MTB AX;

CJX = CHAN 'COMP' 'SCAL' CJX;
CAX = CHAN 'COMP' 'SCAL' CAX;


WX        = CAX * CJX ;

CJY = CHAN 'CHPO' MTB JY;

CAY = CHAN 'CHPO' MTB AY;


CJY = CHAN 'COMP' 'SCAL' CJY;
CAY = CHAN 'COMP' 'SCAL' CAY;

WY        = CAY * CJY ;

CW = WX + WY;

W = CHAN 'CHAM' SU1 CW;

* intégrale de la densité d'énergie magnétique (x2)

L1 = INTG W MTB;

* on multiplie par l'épaisseur car l'intégration est faite en 3D
L = (L1*EPAI);

MOT1 =  CHAI 'INDUCTANCE TROUVEE PAR LES ELEMENTS FINIS:' L;


* on applique la formule analytique
* Rosa, Grover - FORMULAS AND TABLES FOR THE CALCULATION OF
* MUTUAL AND SELF-INDUCTANCE
* Bulletin of Bureau of Standards vol 8 no 1, p.153

AJOUT1 = 0.5+(0.2235*((LY+EPAI)/LX));

LREF = ((2*LX)*100)*(LOG((2*LX)/(LY+EPAI))+AJOUT1)*(1.E-9);



MOT2 =  CHAI 'INDUCTANCE TROUVEE PAR LA FORMULE ANALYTIQUE:' LREF;

MESS MOT1;
MESS MOT2;

*  test de validite
 EC = (L - LREF) ABS;
SI (EC > 1d-10); MESS 'resultat incorrect' ' ' ec; erreur 5; finsi;
fin;