* fichier : pret_gfmp.dgibi
**********************************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la solution des ****
**** Equations d'Euler. ****
**** "Ghost Fluid Method for the Poor" ****
**** ****
**** OPERATEUR PRET ****
**** Operateur qui 'recontruit les variables primitives aux faces ****
**** "Stiffened gas" ****
**** Interieur et murs ****
**** ****
**** A. BECCANTINI DEN/DM2S/SFME/LTMF MAI 2011 ****
**********************************************************************
'OPTION' 'DIME' 2 'ELEM' QUA4 'ECHO' 1
'TRAC' 'X' ;
*
*** GRAPH
*
GRAPH = VRAI ;
GRAPH = FAUX ;
***************************
***** DOMAINE SPATIAL ****
***************************
*
**** Deux carre
*
A1 = 0.0D0 0.0D0;
A2 = 2.0D0 0.0D0;
A3 = 2.2D0 0.0D0;
A4 = 2.0D0 0.8D0;
A5 = 1.0D0 0.90D0;
A6 = 0.0D0 1.1D0;
L12 = A1 'DROIT' 1 A2;
L23 = A2 'DROIT' 1 A3;
L34 = A3 'DROIT' 1 A4;
L45 = A4 'DROIT' 1 A5;
L56 = A5 'DROIT' 1 A6;
L61 = A6 'DROIT' 1 A1;
L25 = A2 'DROIT' 1 A5;
DOM10 = 'DALL' L12 L25 L56 L61
'PLANE';
DOM20 = 'DALL' L23 L34 L45 ('INVERSE' L25)
'PLANE';
*
*** Etats a gauche et a droite
*
*
al1g = 0.9 ; limgal1 = 0.11 ;
rn1g = 1.10 ; limgrn1 = 0.21 ;
ux1g = 130. ; limgux1 = 0.31 ;
uy1g = 111. ; limguy1 = 0.41 ;
pn1g = 1.13E5 ; limgpn1 = 0.51 ;
xgragal1 = 11.0D-2 ; ygragal1 = 12.0D-2 ;
xgragrn1 = 13.0D-1 ; ygragrn1 = 11.0D-1 ;
xgragux1 = 15.0 ; ygragux1 = 12.0 ;
xgraguy1 = 17.0 ; ygraguy1 = 19.0 ;
xgragpn1 = 21.0D2 ; ygragpn1 = 21.2D2 ;
al1d = 0.8 ; limdal1 = 0.17 ;
rn1d = 1.31 ; limdrn1 = 0.19 ;
ux1d = 230. ; limdux1 = 0.22 ;
uy1d = 431. ; limduy1 = 0.32 ;
pn1d = 1.31E5 ; limdpn1 = 0.42 ;
xgradal1 = 11.0D-2 ; ygradal1 = 11.2D-2 ;
xgradrn1 = 21.0D-1 ; ygradrn1 = 21.2D-1 ;
xgradux1 = 31.0 ; ygradux1 = 31.2 ;
xgraduy1 = 41.0 ; ygraduy1 = 41.3 ;
xgradpn1 = 51.0D1 ; ygradpn1 = 51.2D1 ;
ETATG = 'PROG' al1g rn1g ux1g uy1g pn1g ;
ETATD = 'PROG' al1d rn1d ux1d uy1d pn1d ;
DOM1 = DOM10 ;
DOM2 = DOM20 ;
DOMTOT = DOM1 ET DOM2;
'ELIMINATION' (DOM1 ET DOM2) 1D-6;
$DOMTOT = 'MODELISER' DOMTOT 'EULER';
$DOM1 = 'MODELISER' DOM1 'EULER';
$DOM2 = 'MODELISER' DOM2 'EULER';
$L25 = 'MODELISER' L25 'EULER';
$L12 = 'MODELISER' L12 'EULER';
TDOMTOT = 'DOMA' $DOMTOT 'VF';
TDOM1 = 'DOMA' $DOM1 'VF';
TDOM2 = 'DOMA' $DOM2 'VF';
TL25 = 'DOMA' $L25 'VF';
TL12 = 'DOMA' $L12 'VF';
MDOM1 = TDOM1 . 'QUAF' ;
MDOM2 = TDOM2 . 'QUAF' ;
MDOMTOT = TDOMTOT . 'QUAF' ;
ML25 = TL25 . 'QUAF' ;
ML12 = TL12 . 'QUAF' ;
'ELIMINATION' (MDOMTOT ET MDOM1 ET ML25) 0.0001 ;
'ELIMINATION' (MDOMTOT ET MDOM2 ET ML12) 0.0001 ;
P1 = 'POIN' 1 (TL25 . 'CENTRE') ;
P2 = 'POIN' 1 (TL12 . 'CENTRE') ;
'SI' GRAPH;
'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET'
('MANUEL' 'POI1' P1 'COULEUR' 'ROUG') 'ET'
('MANUEL' 'POI1' P2 'COULEUR' 'ROUG'))
'TITRE' 'Domaine et points';
'FINSI' ;
**** La densité, la pression, la vitesse
AL1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL' al1g) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL' al1d) ;
GRADAL1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 2
'P1DX' xgragal1 'P1DY' ygragal1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 2
'P1DX' xgradal1 'P1DY' ygradal1)
;
LIMAL1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'P1' limgal1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'P1' limdal1) ;
RN1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL' rn1g) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL' rn1d) ;
GRADRN1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 2
'P1DX' xgragrn1 'P1DY' ygragrn1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 2
'P1DX' xgradrn1 'P1DY' ygradrn1)
;
LIMRN1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'P1' limgrn1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'P1' limdrn1) ;
PN1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL' pn1g) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL' pn1d) ;
GRADPN1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 2
'P1DX' xgragpn1 'P1DY' ygragpn1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 2
'P1DX' xgradpn1 'P1DY' ygradpn1)
;
LIMPN1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'P1' limgpn1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'P1' limdpn1) ;
UXN1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'UX' ux1g) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'UX' ux1d) ;
GRADUX1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 2
'P1DX' xgragux1 'P1DY' ygragux1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 2
'P1DX' xgradux1 'P1DY' ygradux1)
;
LIMUX1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'P1' limgux1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'P1' limdux1) ;
UYN1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'UY' uy1g) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'UY' uy1d) ;
GRADUY1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 2
'P2DX' xgraguy1 'P2DY' ygraguy1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 2
'P2DX' xgraduy1 'P2DY' ygraduy1)
;
LIMUY1 = ('MANU' 'CHPO' (TDOM1 . 'CENTRE') 1 'P2' limguy1) '+'
('MANU' 'CHPO' (TDOM2 . 'CENTRE') 1 'P2' limduy1) ;
VN1 = UXN1 '+' UYN1 ;
GRADVN1 = GRADUX1 '+' GRADUY1 ;
LIMVN1 = LIMUX1 '+' LIMUY1 ;
*************************************************************************
*** We check if the reconstruction degenerates to first order as the ****
*** gradients or the limiters are zero ****
*************************************************************************
NESP = 'MOTS' ;
AL1G0 RN1G0 VN1G0 PN1G0 =
'PRET' 'GFMP' NESP $DOMTOT
AL1 (0.0 * GRADAL1) LIMAL1
RN1 (0.0 * GRADRN1) LIMRN1
VN1 (0.0 * GRADVN1) LIMVN1
PN1 (0.0 * GRADPN1) LIMPN1
;
AL1L0 RN1L0 VN1L0 PN1L0 =
'PRET' 'GFMP' NESP $DOMTOT
AL1 GRADAL1 (0.0 * LIMAL1)
RN1 GRADRN1 (0.0 * LIMRN1)
VN1 GRADVN1 (0.0 * LIMVN1)
PN1 GRADPN1 (0.0 * LIMPN1)
;
ERRO = 'PROG' ;
ERRO = ERRO 'ET' ('PROG' ('MAXIMUM' (AL1G0 '-' AL1L0) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ('MAXIMUM' (RN1G0 '-' RN1L0) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ('MAXIMUM' (VN1G0 '-' VN1L0) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ('MAXIMUM' (PN1G0 '-' PN1L0) 'ABS'));
ERRO = 'MAXIMUM' ERRO 'ABS' ;
'SI' (ERRO > 1.0D-16) ;
'MESSAGE' 'Limiteur et gradients zero' ;
'MESSAGE' 'Erreur trop importante' ;
'MESSAGE' ('CHAINE' 'erro = ' ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
*
*** Control des etats sur la surface qui contient P1
*
GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
GEP1F = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
AL1GEOP1 = 'REDU' AL1G0 GEOP1;
RN1GEOP1 = 'REDU' RN1G0 GEOP1;
PN1GEOP1 = 'REDU' PN1G0 GEOP1;
VN1GEOP1 = 'REDU' VN1G0 GEOP1;
TN1GEP1F = 'REDU' VN1G0 GEP1F;
al1gp1 = 'EXTRAIRE' AL1GEOP1 'SCAL' 1 1 1 ;
al1dp1 = 'EXTRAIRE' AL1GEOP1 'SCAL' 1 1 3 ;
rn1gp1 = 'EXTRAIRE' RN1GEOP1 'SCAL' 1 1 1 ;
rn1dp1 = 'EXTRAIRE' RN1GEOP1 'SCAL' 1 1 3 ;
pn1gp1 = 'EXTRAIRE' PN1GEOP1 'SCAL' 1 1 1 ;
pn1dp1 = 'EXTRAIRE' PN1GEOP1 'SCAL' 1 1 3 ;
un1gp1 = 'EXTRAIRE' VN1GEOP1 'UN ' 1 1 1 ;
un1dp1 = 'EXTRAIRE' VN1GEOP1 'UN ' 1 1 3 ;
ut1gp1 = 'EXTRAIRE' VN1GEOP1 'UT ' 1 1 1 ;
ut1dp1 = 'EXTRAIRE' VN1GEOP1 'UT ' 1 1 3 ;
nxp1 = 'EXTRAIRE' TN1GEP1F 'NX ' 1 1 1 ;
nyp1 = 'EXTRAIRE' TN1GEP1F 'NY ' 1 1 1 ;
txp1 = 'EXTRAIRE' TN1GEP1F 'TX ' 1 1 1 ;
typ1 = 'EXTRAIRE' TN1GEP1F 'TY ' 1 1 1 ;
ux1gp1 = (un1gp1 * nxp1) '+' (ut1gp1 * txp1) ;
uy1gp1 = (un1gp1 * nyp1) '+' (ut1gp1 * typ1) ;
ux1dp1 = (un1dp1 * nxp1) '+' (ut1dp1 * txp1) ;
uy1dp1 = (un1dp1 * nyp1) '+' (ut1dp1 * typ1) ;
ERRO = 'PROG' ;
ERRO = ERRO 'ET' ('PROG' ((al1gp1 '-' al1g) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((al1dp1 '-' al1d) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((rn1gp1 '-' rn1g) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((rn1dp1 '-' rn1d) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((pn1gp1 '-' pn1g) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((pn1dp1 '-' pn1d) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((ux1gp1 '-' ux1g) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((ux1dp1 '-' ux1d) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((uy1gp1 '-' uy1g) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((uy1dp1 '-' uy1d) 'ABS'));
ERRO = 'MAXIMUM' ERRO 'ABS' ;
'SI' (ERRO > 1.0D-12) ;
'MESSAGE' 'Limiteur et gradients zero' ;
'MESSAGE' 'Les etats' ;
'MESSAGE' 'Erreur trop importante' ;
'MESSAGE' ('CHAINE' 'erro = ' ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
*
**** Wall
*
GEOP2 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P2;
GEP2F = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P2;
AL1GEOP2 = 'REDU' AL1G0 GEOP2;
RN1GEOP2 = 'REDU' RN1G0 GEOP2;
PN1GEOP2 = 'REDU' PN1G0 GEOP2;
VN1GEOP2 = 'REDU' VN1G0 GEOP2;
TN1GEP2F = 'REDU' VN1G0 GEP2F;
al1gp2 = 'EXTRAIRE' AL1GEOP2 'SCAL' 1 1 1 ;
al1dp2 = 'EXTRAIRE' AL1GEOP2 'SCAL' 1 1 3 ;
rn1gp2 = 'EXTRAIRE' RN1GEOP2 'SCAL' 1 1 1 ;
rn1dp2 = 'EXTRAIRE' RN1GEOP2 'SCAL' 1 1 3 ;
pn1gp2 = 'EXTRAIRE' PN1GEOP2 'SCAL' 1 1 1 ;
pn1dp2 = 'EXTRAIRE' PN1GEOP2 'SCAL' 1 1 3 ;
un1gp2 = 'EXTRAIRE' VN1GEOP2 'UN ' 1 1 1 ;
un1dp2 = 'EXTRAIRE' VN1GEOP2 'UN ' 1 1 3 ;
ut1gp2 = 'EXTRAIRE' VN1GEOP2 'UT ' 1 1 1 ;
ut1dp2 = 'EXTRAIRE' VN1GEOP2 'UT ' 1 1 3 ;
ERRO = 'PROG' ;
ERRO = ERRO 'ET' ('PROG' ((al1gp2 '-' al1dp2) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((rn1gp2 '-' rn1dp2) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((pn1gp2 '-' pn1dp2) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((ut1gp2 '-' ut1dp2) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((un1gp2 '+' un1dp2) 'ABS'));
ERRO = 'MAXIMUM' ERRO 'ABS' ;
'SI' (ERRO > 1.0D-16) ;
'MESSAGE' 'Limiteur et gradients zero' ;
'MESSAGE' 'Les etats (mur)' ;
'MESSAGE' 'Erreur trop importante' ;
'MESSAGE' ('CHAINE' 'erro = ' ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
*****************************************************
*** We check the second order reconstruction. *******
*****************************************************
AL1F RN1F VN1F PN1F =
'PRET' 'GFMP' NESP $DOMTOT
AL1 GRADAL1 LIMAL1
RN1 GRADRN1 LIMRN1
VN1 GRADVN1 LIMVN1
PN1 GRADPN1 LIMPN1
;
*** Internal point
GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
GEP1F = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
AL1GEOP1 = 'REDU' AL1F GEOP1;
RN1GEOP1 = 'REDU' RN1F GEOP1;
PN1GEOP1 = 'REDU' PN1F GEOP1;
VN1GEOP1 = 'REDU' VN1F GEOP1;
TN1GEP1F = 'REDU' VN1F GEP1F;
al1gp1 = 'EXTRAIRE' AL1GEOP1 'SCAL' 1 1 1 ;
al1dp1 = 'EXTRAIRE' AL1GEOP1 'SCAL' 1 1 3 ;
rn1gp1 = 'EXTRAIRE' RN1GEOP1 'SCAL' 1 1 1 ;
rn1dp1 = 'EXTRAIRE' RN1GEOP1 'SCAL' 1 1 3 ;
pn1gp1 = 'EXTRAIRE' PN1GEOP1 'SCAL' 1 1 1 ;
pn1dp1 = 'EXTRAIRE' PN1GEOP1 'SCAL' 1 1 3 ;
un1gp1 = 'EXTRAIRE' VN1GEOP1 'UN ' 1 1 1 ;
un1dp1 = 'EXTRAIRE' VN1GEOP1 'UN ' 1 1 3 ;
ut1gp1 = 'EXTRAIRE' VN1GEOP1 'UT ' 1 1 1 ;
ut1dp1 = 'EXTRAIRE' VN1GEOP1 'UT ' 1 1 3 ;
nxp1 = 'EXTRAIRE' TN1GEP1F 'NX ' 1 1 1 ;
nyp1 = 'EXTRAIRE' TN1GEP1F 'NY ' 1 1 1 ;
txp1 = 'EXTRAIRE' TN1GEP1F 'TX ' 1 1 1 ;
typ1 = 'EXTRAIRE' TN1GEP1F 'TY ' 1 1 1 ;
ux1gp1 = (un1gp1 * nxp1) '+' (ut1gp1 * txp1) ;
uy1gp1 = (un1gp1 * nyp1) '+' (ut1gp1 * typ1) ;
ux1dp1 = (un1dp1 * nxp1) '+' (ut1dp1 * txp1) ;
uy1dp1 = (un1dp1 * nyp1) '+' (ut1dp1 * typ1) ;
* Selection of the left and the right state
MAICEN = 'DIFF' ('CHANGER' 'POI1' GEOP1) GEP1F ;
'SI' (('NBNO' MAICEN) 'NEG' 2) ;
'MESSAGE' 'Probleme en maicen' ;
'ERREUR' 5 ;
'FINSI' ;
PC1 = 'POIN' 1 MAICEN ;
PC2 = 'POIN' 2 MAICEN ;
DXPC1 = ('COORDONNEE' 1 P1) '-' ('COORDONNEE' 1 PC1) ;
DYPC1 = ('COORDONNEE' 2 P1) '-' ('COORDONNEE' 2 PC1) ;
DXPC2 = ('COORDONNEE' 1 P1) '-' ('COORDONNEE' 1 PC2) ;
DYPC2 = ('COORDONNEE' 2 P1) '-' ('COORDONNEE' 2 PC2) ;
PSPC1 = (DXPC1 * nxp1) '+' (DYPC1 * nyp1) ;
PSPC2 = (DXPC2 * nxp1) '+' (DYPC2 * nyp1) ;
LOG1 = (PSPC1 > 0.0) 'ET' (PSPC2 < 0.0) ;
'SI' LOG1 ;
PG = PC1 ;
DXG = DXPC1 ;
DYG = DYPC1 ;
PD = PC2 ;
DXD = DXPC2 ;
DYD = DYPC2 ;
'FINSI' ;
LOG2 = (PSPC1 < 0.0) 'ET' (PSPC2 > 0.0) ;
'SI' LOG2 ;
PG = PC2 ;
DXG = DXPC2 ;
DYG = DYPC2 ;
PD = PC1 ;
DXD = DXPC1 ;
DYD = DYPC1 ;
'FINSI' ;
'SI' (('NON' LOG1) 'ET' ('NON' LOG2)) ;
'MESSAGE' 'Probleme en facel' ;
'ERREUR' 5 ;
'FINSI' ;
al1gn = ('EXTRAIRE' AL1 'SCAL' PG) '+' (
(
(('EXTRAIRE' GRADAL1 'P1DX' PG) * DXG) '+'
(('EXTRAIRE' GRADAL1 'P1DY' PG) * DYG)
) * ('EXTRAIRE' LIMAL1 'P1' PG)) ;
rn1gn = ('EXTRAIRE' RN1 'SCAL' PG) '+' (
(
(('EXTRAIRE' GRADRN1 'P1DX' PG) * DXG) '+'
(('EXTRAIRE' GRADRN1 'P1DY' PG) * DYG)
) * ('EXTRAIRE' LIMRN1 'P1' PG)) ;
pn1gn = ('EXTRAIRE' PN1 'SCAL' PG) '+' (
(
(('EXTRAIRE' GRADPN1 'P1DX' PG) * DXG) '+'
(('EXTRAIRE' GRADPN1 'P1DY' PG) * DYG)
) * ('EXTRAIRE' LIMPN1 'P1' PG)) ;
ux1gn = ('EXTRAIRE' VN1 'UX ' PG) '+' (
(
(('EXTRAIRE' GRADUX1 'P1DX' PG) * DXG) '+'
(('EXTRAIRE' GRADUX1 'P1DY' PG) * DYG)
) * ('EXTRAIRE' LIMUX1 'P1' PG)) ;
uy1gn = ('EXTRAIRE' VN1 'UY ' PG) '+' (
(
(('EXTRAIRE' GRADUY1 'P2DX' PG) * DXG) '+'
(('EXTRAIRE' GRADUY1 'P2DY' PG) * DYG)
) * ('EXTRAIRE' LIMUY1 'P2' PG)) ;
ERRO = 'PROG' ;
ERRO = ERRO 'ET' ('PROG' ((al1gp1 '-' al1gn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((rn1gp1 '-' rn1gn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((pn1gp1 '-' pn1gn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((ux1gp1 '-' ux1gn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((uy1gp1 '-' uy1gn) 'ABS'));
ERRO = 'MAXIMUM' ERRO 'ABS' ;
'SI' (ERRO > 1.0D-13) ;
'MESSAGE' 'Deuxieme ordre' ;
'MESSAGE' 'Les etats gauches' ;
'MESSAGE' 'Erreur trop importante' ;
'MESSAGE' ('CHAINE' 'erro = ' ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
al1dn = ('EXTRAIRE' AL1 'SCAL' PD) '+' (
(
(('EXTRAIRE' GRADAL1 'P1DX' PD) * DXD) '+'
(('EXTRAIRE' GRADAL1 'P1DY' PD) * DYD)
) * ('EXTRAIRE' LIMAL1 'P1' PD)) ;
rn1dn = ('EXTRAIRE' RN1 'SCAL' PD) '+' (
(
(('EXTRAIRE' GRADRN1 'P1DX' PD) * DXD) '+'
(('EXTRAIRE' GRADRN1 'P1DY' PD) * DYD)
) * ('EXTRAIRE' LIMRN1 'P1' PD)) ;
pn1dn = ('EXTRAIRE' PN1 'SCAL' PD) '+' (
(
(('EXTRAIRE' GRADPN1 'P1DX' PD) * DXD) '+'
(('EXTRAIRE' GRADPN1 'P1DY' PD) * DYD)
) * ('EXTRAIRE' LIMPN1 'P1' PD)) ;
ux1dn = ('EXTRAIRE' VN1 'UX ' PD) '+' (
(
(('EXTRAIRE' GRADUX1 'P1DX' PD) * DXD) '+'
(('EXTRAIRE' GRADUX1 'P1DY' PD) * DYD)
) * ('EXTRAIRE' LIMUX1 'P1' PD)) ;
uy1dn = ('EXTRAIRE' VN1 'UY ' PD) '+' (
(
(('EXTRAIRE' GRADUY1 'P2DX' PD) * DXD) '+'
(('EXTRAIRE' GRADUY1 'P2DY' PD) * DYD)
) * ('EXTRAIRE' LIMUY1 'P2' PD)) ;
ERRO = 'PROG' ;
ERRO = ERRO 'ET' ('PROG' ((al1dp1 '-' al1dn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((rn1dp1 '-' rn1dn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((pn1dp1 '-' pn1dn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((ux1dp1 '-' ux1dn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((uy1dp1 '-' uy1dn) 'ABS'));
ERRO = 'MAXIMUM' ERRO 'ABS' ;
'SI' (ERRO > 1.0D-13) ;
'MESSAGE' 'Deuxieme ordre' ;
'MESSAGE' 'Les etats droites' ;
'MESSAGE' 'Erreur trop importante' ;
'MESSAGE' ('CHAINE' 'erro = ' ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
*** Wall point
GEOP2 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P2;
GEP2F = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P2;
'SI' GRAPH ;
'TRACER' (DOMTOT 'ET' ('MANUEL' 'POI1' P2 'COULEUR' 'ROUG'))
'TITRE' 'Wall' ;
'FINSI' ;
AL1GEOP2 = 'REDU' AL1F GEOP2;
RN1GEOP2 = 'REDU' RN1F GEOP2;
PN1GEOP2 = 'REDU' PN1F GEOP2;
VN1GEOP2 = 'REDU' VN1F GEOP2;
TN1GEP2F = 'REDU' VN1F GEP2F;
al1gp2 = 'EXTRAIRE' AL1GEOP2 'SCAL' 1 1 1 ;
al1dp2 = 'EXTRAIRE' AL1GEOP2 'SCAL' 1 1 3 ;
rn1gp2 = 'EXTRAIRE' RN1GEOP2 'SCAL' 1 1 1 ;
rn1dp2 = 'EXTRAIRE' RN1GEOP2 'SCAL' 1 1 3 ;
pn1gp2 = 'EXTRAIRE' PN1GEOP2 'SCAL' 1 1 1 ;
pn1dp2 = 'EXTRAIRE' PN1GEOP2 'SCAL' 1 1 3 ;
un1gp2 = 'EXTRAIRE' VN1GEOP2 'UN ' 1 1 1 ;
un1dp2 = 'EXTRAIRE' VN1GEOP2 'UN ' 1 1 3 ;
ut1gp2 = 'EXTRAIRE' VN1GEOP2 'UT ' 1 1 1 ;
ut1dp2 = 'EXTRAIRE' VN1GEOP2 'UT ' 1 1 3 ;
nxp2 = 'EXTRAIRE' TN1GEP2F 'NX ' 1 1 1 ;
nyp2 = 'EXTRAIRE' TN1GEP2F 'NY ' 1 1 1 ;
txp2 = 'EXTRAIRE' TN1GEP2F 'TX ' 1 1 1 ;
typ2 = 'EXTRAIRE' TN1GEP2F 'TY ' 1 1 1 ;
ux1gp2 = (un1gp2 * nxp2) '+' (ut1gp2 * txp2) ;
uy1gp2 = (un1gp2 * nyp2) '+' (ut1gp2 * typ2) ;
* Selection of the left and the right state
MAICEN = 'DIFF' ('CHANGER' 'POI1' GEOP2) GEP2F ;
'SI' (('NBNO' MAICEN) 'NEG' 1) ;
'MESSAGE' 'Probleme en maicen' ;
'ERREUR' 5 ;
'FINSI' ;
PC1 = 'POIN' 1 MAICEN ;
DXPC1 = ('COORDONNEE' 1 P2) '-' ('COORDONNEE' 1 PC1) ;
DYPC1 = ('COORDONNEE' 2 P2) '-' ('COORDONNEE' 2 PC1) ;
PSPC1 = (DXPC1 * nxp2) '+' (DYPC1 * nyp2) ;
LOG1 = (PSPC1 > 0.0) ;
'SI' LOG1 ;
PG = PC1 ;
DXG = DXPC1 ;
DYG = DYPC1 ;
'FINSI' ;
'SI' (('NON' LOG1)) ;
'MESSAGE' 'Probleme en facel' ;
'ERREUR' 5 ;
'FINSI' ;
al1gn = ('EXTRAIRE' AL1 'SCAL' PG) '+' (
(
(('EXTRAIRE' GRADAL1 'P1DX' PG) * DXG) '+'
(('EXTRAIRE' GRADAL1 'P1DY' PG) * DYG)
) * ('EXTRAIRE' LIMAL1 'P1' PG)) ;
rn1gn = ('EXTRAIRE' RN1 'SCAL' PG) '+' (
(
(('EXTRAIRE' GRADRN1 'P1DX' PG) * DXG) '+'
(('EXTRAIRE' GRADRN1 'P1DY' PG) * DYG)
) * ('EXTRAIRE' LIMRN1 'P1' PG)) ;
pn1gn = ('EXTRAIRE' PN1 'SCAL' PG) '+' (
(
(('EXTRAIRE' GRADPN1 'P1DX' PG) * DXG) '+'
(('EXTRAIRE' GRADPN1 'P1DY' PG) * DYG)
) * ('EXTRAIRE' LIMPN1 'P1' PG)) ;
ux1gn = ('EXTRAIRE' VN1 'UX ' PG) '+' (
(
(('EXTRAIRE' GRADUX1 'P1DX' PG) * DXG) '+'
(('EXTRAIRE' GRADUX1 'P1DY' PG) * DYG)
) * ('EXTRAIRE' LIMUX1 'P1' PG)) ;
uy1gn = ('EXTRAIRE' VN1 'UY ' PG) '+' (
(
(('EXTRAIRE' GRADUY1 'P2DX' PG) * DXG) '+'
(('EXTRAIRE' GRADUY1 'P2DY' PG) * DYG)
) * ('EXTRAIRE' LIMUY1 'P2' PG)) ;
ERRO = 'PROG' ;
ERRO = ERRO 'ET' ('PROG' ((al1gp2 '-' al1gn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((rn1gp2 '-' rn1gn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((pn1gp2 '-' pn1gn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((ux1gp2 '-' ux1gn) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((uy1gp2 '-' uy1gn) 'ABS'));
ERRO = 'MAXIMUM' ERRO 'ABS' ;
'SI' (ERRO > 1.0D-13) ;
'MESSAGE' 'Deuxieme ordre' ;
'MESSAGE' 'Les etats gauches au mur' ;
'MESSAGE' 'Erreur trop importante' ;
'MESSAGE' ('CHAINE' 'erro = ' ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
ERRO = 'PROG' ;
ERRO = ERRO 'ET' ('PROG' ((al1gp2 '-' al1dp2) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((rn1gp2 '-' rn1dp2) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((pn1gp2 '-' pn1dp2) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((un1gp2 '+' un1dp2) 'ABS'));
ERRO = ERRO 'ET' ('PROG' ((ut1gp2 '-' ut1dp2) 'ABS'));
ERRO = 'MAXIMUM' ERRO 'ABS' ;
'SI' (ERRO > 1.0D-13) ;
'MESSAGE' 'Deuxieme ordre' ;
'MESSAGE' 'Les etats gauches au mur' ;
'MESSAGE' 'Erreur trop importante' ;
'MESSAGE' ('CHAINE' 'erro = ' ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
* One species involved
* For the sake of simplicity, we put Y = alpha = phi
* and we test with respect to the previous case.
*
PH1F RN1F VN1F PN1F =
'PRET' 'GFMP' NESP $DOMTOT
AL1 GRADAL1 LIMAL1
RN1 GRADRN1 LIMRN1
VN1 GRADVN1 LIMVN1
PN1 GRADPN1 LIMPN1
;
PH0F RN1F VN1F PN1F =
'PRET' 'GFMP' NESP $DOMTOT
AL1 GRADAL1 (LIMAL1 * 0.0D0)
RN1 GRADRN1 LIMRN1
VN1 GRADVN1 LIMVN1
PN1 GRADPN1 LIMPN1
;
NESP1 = 'MOTS' 'ESP2' ;
PHFM RN1F VN1F PN1F YNF ALF =
'PRET' 'GFMP' NESP1 $DOMTOT
AL1 GRADAL1 (LIMAL1 * 0.0D0)
RN1 GRADRN1 LIMRN1
VN1 GRADVN1 LIMVN1
PN1 GRADPN1 LIMPN1
('NOMC' AL1 'ESP2') (0.0 * GRADAL1) LIMAL1
('NOMC' AL1 'ESP2') GRADAL1 LIMAL1
;
PHFM = 'NOMC' 'ESP2' PHFM ;
PH0F = 'NOMC' 'ESP2' PH0F ;
PH1F = 'NOMC' 'ESP2' PH1F ;
ERRO = 'MAXIMUM' (PHFM '-' YNF) 'ABS' ;
ERRO = ERRO '+' (
'MAXIMUM' (PH0F '-' YNF) 'ABS'
) ;
ERRO = ERRO '+' (
'MAXIMUM' (ALF '-' PH1F) 'ABS'
) ;
PHFM RN1F VN1F PN1F YNF1 ALF1 =
'PRET' 'GFMP' NESP1 $DOMTOT
AL1 GRADAL1 (LIMAL1 * 0.0D0)
RN1 GRADRN1 LIMRN1
VN1 GRADVN1 LIMVN1
PN1 GRADPN1 LIMPN1
('NOMC' AL1 'ESP2') (1.0 * GRADAL1) LIMAL1
('NOMC' AL1 'ESP2') (0.0 * GRADAL1) LIMAL1
;
ERRO = ERRO '+' (
'MAXIMUM' (ALF1 '-' YNF) 'ABS'
) ;
ERRO = ERRO '+' (
'MAXIMUM' (ALF '-' YNF1) 'ABS'
) ;
'SI' (ERRO > 1.0D-13) ;
'MESSAGE' 'Multi-espece' ;
'MESSAGE' 'Erreur trop importante' ;
'MESSAGE' ('CHAINE' 'erro = ' ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer