OPTION 'TRAC' PSC;
nomfic = 'AILE00.ps';
OPTION 'FTRAC' nomfic;
*  OK janvier 2003
*
*                  Z
*
*
*
*   A---------------0------------------B
*
*
*
*
*                  N

TITR 'Aile';
OPTI DIME 2 ELEM TRI6 MODE PLAN;
* parametres
Ld=5.0;

*points
pA=(-1.0*Ld) 0.0;
pB= Ld 0.0;
pZ= 0.0 Ld;
pN= 0.0 (-1.0*Ld);
p0= 0.0 0.0;

*nbr de points
nL=35;
nR=25;
nA=2;
nE=25;

* cercle au loin
LOINhd= CERC nL pB p0 pZ;
LOINhg= CERC nL pZ p0 pA;
LOINbg= CERC nL pA p0 pN;
LOINbd= CERC nL pN p0 pB;
CNT=LOINhd et LOINhg et LOINbg et LOINbd;

* Aile
bd1=0.0 -0.1;
bd2=0. 0.1;
bd0=-0.1 0.0;
bdf=1.0 0.0;

dh = droite pZ  nR bd2;
db = droite bd1 nR pN;
bdah= CERC nA bd1 p0 bd0;
bdab= CERC nA bd2 p0 bd0;
bda= bdah et bdab;

intrd=droi bdf bd1 nE;
vintrd=bdah et intrd;
extrd=droi bd2 bdf nE;
vextrd=bdab et extrd; 
Aile=   vextrd et vintrd;

* surface Maillage  
adroite=LOINbd et LOINhd et dh et extrd et intrd et db;
agauche=LOINhg  et LOINbg et db et dh et bda;
SURFd= SURF(adroite);
SURFg= SURF(agauche);
su0=SURFd et SURFg;    
TRAC su0 'TITR' 'Aile';
 

************************FIN DU MAILLAGE

*Modele
MOD1=MODE su0 THERMIQUE ISOTROPE ;
* DONNEES MATERIAU ET ELEMENTS
MAT1=MATR MOD1 K 1.0 ;

* MATRICES DE CONDUCTIVITE
COND1=CONDUCTIVITE MOD1 MAT1 ;
*
* CONDITIONS AUX LIMITES
*BLOQue/ DEPlacement Impose!
CLS =BLOQ CNT T ;
DCLS=DEPI CLS 0. ;
CLL =BLOQ Aile T ;
DCLL=DEPI CLL 1. ;

***************************

***************************
* RESOLUTION  0 pbm homogene: 

CONDTOT1=COND1 ET CLS ET CLL ;
FLUTOT1= DCLL ET DCLS;
 
* RESOLUTION
PSI1=RESO CONDTOT1 FLUTOT1 ;
*
* POST-TRAITEMENT
TITR 'aile PSI';
TRAC PSI1 su0 ;
*
* CONDITIONS AUX LIMITES
*BLOQue/ DEPlacement Impose!
yy=coord 2 CNT;
xx=coord 1 CNT;
psiCNT=yy +(-0.0*xx);
CLS =BLOQ CNT T ;
DCLS=DEPI CLS psiCNT ;
CLL =BLOQ Aile T ;
DCLL=DEPI CLL 0. ;

***************************
* RESOLUTION
CONDTOT0=COND1 ET CLS ET CLL ;
FLUTOT0= DCLL ET DCLS;
* RESOLUTION
PSI0=RESO CONDTOT0 FLUTOT0 ;

*
* POST-TRAITEMENT
TITR 'aile PSI';
TRAC PSI0 su0 CNT  ;
** TRACE des vitesses du probleme 1
GRADPSI1 = GRAD MOD1 PSI1 MAT1;
 TITR 'aile GRADPSI1';
*TRAC  MOD1 GRADPSI1   su0;
totox = exco T,X GRADPSI1 ;
totoy = exco T,Y GRADPSI1;
totoyx = (nomc ux totoy) et (nomc uy (-1. * totox) ) ; 
totoyxn = chan chpo mod1 totoyx ; 
vec1 = vect totoyxn UX UY .1 bleu;    
trac vec1 su0 Aile; 

 
** TRACE des vitesses
GRADPSI0 = GRAD MOD1 PSI0 MAT1;
 TITR 'aile GRADPSI0';
*TRAC  MOD1 GRADPSI0   su0;
totox = exco T,X GRADPSI0 ;
totoy = exco T,Y GRADPSI0;
totoyx = (nomc ux totoy) et (nomc uy (-1. * totox) ) ; 
totoyxn = chan chpo mod1 totoyx ; 
vec0 = vect totoyxn UX UY .25 jaune;    
trac vec0 su0 ; 

 
************KUTTA*******
eps=0.025;
sbfp=su0 poin proc  ((1.-eps) eps);
 zzp=extr PSI0 T sbfp;
 zp= extr PSI1 T sbfp;
 list zzp;
sbfm=su0 poin proc ((1.-eps) (-1.*eps));;
 zzm=extr PSI0 T sbfm;
 zm= extr PSI1 T sbfm;

Gamma = -1.0* (zzp+zzm)/(zm+zp-2.);
list  Gamma;


**********recomposition
psi= Gamma*psi1 + psi0;
TITR 'aile PSI';
*TRAC PSI su0 ;


********** gradient pour vitesse
GRADPSI = GRAD MOD1 PSI MAT1;
*TITR 'aile GRADPSI';
*TRAC  MOD1 GRADPSI   su0;


totox = exco T,X gradpsi;
totoy = exco T,Y gradpsi;

****PRESSION..... BERNOULLI
gr1=chan CHPO mod1 gradpsi;
c1=nomc SCAL (exco T,X gr1);
c2=nomc SCAL (exco T,Y gr1);

 pr0 = -0.5*((c1 * c1) + (c2 * c2)) ;     
  TITR 'aile Press';
 trac mod1 pr0 su0 cnt 20 ;    
  
*** PRESSION SUR L'EXTRADOS
pextr=INT_COMP su0 pr0 vextrd;
*DESS(evol CHPO pextr vextrd) 'TITRE' 'p extrados';
*** PRESSION SUR L'INTRADOS
pintr=INT_COMP su0 pr0 vintrd;
*DESS(evol CHPO pintr vintrd)  'TITRE' 'p intrados';

 
*****vitesse finale
totoyx = (nomc UX totoy) et (nomc UY (-1. * totox) ) ; 
totoyxn = chan chpo mod1 totoyx ; 
vec1 = vect totoyxn UX UY .05 bleu;    
trac vec1 su0 CNT;         
  
*******
 trac psi su0 Aile 20;     
***
XAE= evol CHPO (coord 1 vextrd) vextrd;
YAE= evol CHPO (coord 2 vextrd) vextrd;
*DESS (EVOL MANU (EXTR XAE ORDO)(EXTR YAE ORDO)) titre 'aile extrados';
XAI= evol CHPO (coord 1 vintrd) vintrd;
YAI= evol CHPO (coord 2 vintrd) vintrd;
*DESS (EVOL MANU (EXTR XAI ORDO)(EXTR YAI ORDO)) titre 'aile intrados';

evaile= (EVOL bleu MANU (EXTR XAI ORDO)(EXTR YAI ORDO)) 
   et (EVOL rouge MANU (EXTR XAE ORDO)(EXTR YAE ORDO));
*DESS(evaile) titre 'aile';
****************** trace pressions
EVOPE=evol CHPO pextr vextrd;
EVOPI=evol CHPO pintr vintrd;
evpail=(EVOL rouge MANU (EXTR XAE ORDO)(EXTR EVOPE ORDO)  )  et 
     (EVOL bleu MANU (EXTR XAI ORDO)(EXTR EVOPI ORDO));
DESS(evaile et evpail) titre 'profil et Pe (rouge) et Pint (bleu)';
     


******GNUPLOT
*enlever a la main les $$ dans les fichier  .OUT
* lancer gnuplot:
* set key bottom
* p'CPextr.OUT't'extr'w l,'CPextr.OUT'u 1:3t'Pe'w l 1,'CPintr.OUT't'intr'w l 3,'CPintr.OUT'u 1:3t'Pi'w l 3
np= DIME (EXTR XAE ORDO);
I=0;
OPTI IMPR 10 IMPR 'CPextr.OUT';
REPETER BCL1 (np-1);
 I=I+1;
 Xi = EXTR (EXTR XAE ORDO) I;
 Yi = EXTR (EXTR YAE ORDO) I;
 Pyi= EXTR (EXTR EVOPE ORDO) I;
 MESSAGE Xi ' ' Yi ' ' Pyi;
FIN BCL1;

******GNUPLOT
np= DIME (EXTR XAI ORDO);
I=0;
OPTI IMPR 10 IMPR 'CPintr.OUT';
REPETER BCL1 (np);
 I=I+1;
 Xi = EXTR (EXTR XAI ORDO) I;
 Yi = EXTR (EXTR YAI ORDO) I;
 Pyi= EXTR (EXTR EVOPI ORDO) I;
 MESSAGE Xi ' ' Yi ' ' Pyi;
FIN BCL1;

OPTI IMPR 6;