* fichier : pret3D2.dgibi
**********************************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la solution des ****
**** Equations d'Euler pour un gaz parfait. ****
**** ****
**** OPERATEUR PRET ****
**** Operateur qui 'recontruit les variables primitives aux faces ****
**** Cas gaz multiespeces "thermally perfect" ****
**** Differents cas tests ****
**** ****
**** A. BECCANTINI DRN/DMT/SEMT/TTMF AOUT 1998 ****
**********************************************************************
'OPTION' 'DIME' 3 ;
'OPTION' 'ELEM' CUB8 ;
'OPTION' 'ECHO' 0 ;
'OPTION' 'TRAC' 'X' ;
*
*** GRAPH
*
*GRAPH = VRAI ;
GRAPH = FAUX ;
***************************
*** PROCEDURE RUO1 ****
***************************
*
*** Precedure pour la rotation d'un tenseur de premier
* ordre autour de l'axe Oz
*
'DEBPROC' RUO1 ;
'ARGUMENT' ALPHA*'FLOTTANT' UN*'FLOTTANT' UT*'FLOTTANT';
SINA = 'SIN' ALPHA ;
COSA = 'COS' ALPHA ;
UX = (UN * COSA ) '-' (UT * SINA) ;
UY = (UN * SINA ) '+' (UT * COSA) ;
'FINPROC' UX UY;
***************************
*** PROCEDURE RUO2 ****
***************************
*
*** Precedure pour la rotation d'un tenseur de deuxieme
* ordre autour de l'axe Oz
*
'DEBPROC' RUO2 ;
'ARGUMENT' ALPHA*'FLOTTANT'
UNN*'FLOTTANT' UNT*'FLOTTANT' UNZ*'FLOTTANT'
UTN*'FLOTTANT' UTT*'FLOTTANT' UTZ*'FLOTTANT';
*
**** (n,t) -> (x,y)
*
* n = CA x '+' SA y ;
* t = -SA x '+' CA y ;
*
* UNT = DUN/DT
*
SA = 'SIN' ALPHA ;
CA = 'COS' ALPHA ;
CA2 = CA * CA ;
CASA = CA * SA ;
SA2 = SA * SA;
UXX = (CA2 * UNN) '-' (CASA * (UNT '+' UTN)) '+' (SA2 *UTT) ;
UYX = (CASA * (UNN '-' UTT)) '-' (SA2 * UNT) '+' (CA2 * UTN ) ;
UXY = (CASA * (UNN '-' UTT)) '+' (CA2 * UNT) '-' (SA2 * UTN ) ;
UYY = (SA2 * UNN) '+' (CASA * (UNT '+' UTN)) '+' (CA2 *UTT) ;
UXZ = (CA * UNZ) '-' (SA * UTZ) ;
UYZ = (SA * UNZ) '+' (CA * UTZ) ;
'FINPROC' UXX UXY UXZ UYX UYY UYZ;
***************************
***** DOMAINE SPATIAL ****
***************************
*
**** Deux carre
*
A1 = 0.0D0 0.0D0 0.0D0;
A2 = 1.0D0 0.0D0 0.0D0;
A3 = 2.0D0 0.0D0 0.0D0;
A4 = 2.0D0 1.0D0 0.0D0;
A5 = 1.0D0 1.0D0 0.0D0;
A6 = 0.0D0 1.0D0 0.0D0;
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;
DOM100 = 'DALL' L12 L25 L56 L61
'PLANE';
DOM200 = 'DALL' L23 L34 L45 ('INVERSE' L25)
'PLANE';
DOM10 = DOM100 'VOLU' 1 'TRAN' (0.0 0.0 1.0);
DOM20 = DOM200 'VOLU' 1 'TRAN' (0.0 0.0 1.0);
*
*** Point face entre le deux carre, ou on fait les controlles
*
P10 = 1.0 0.5 0.5;
DOM1 = DOM10 ;
DOM2 = DOM20 ;
'ELIMINATION' (DOM1 ET DOM2) 1D-6;
DOMTOT = DOM1 ET DOM2;
$DOMTOT = 'MODELISER' DOMTOT 'EULER';
$DOM1 = 'MODELISER' DOM1 'EULER';
$DOM2 = 'MODELISER' DOM2 'EULER';
TDOMTOT = 'DOMA' $DOMTOT 'VF';
TDOM1 = 'DOMA' $DOM1 'VF';
TDOM2 = 'DOMA' $DOM2 'VF';
MDOM1 = TDOM1 . 'QUAF' ;
MDOM2 = TDOM2 . 'QUAF' ;
MDOMTOT = TDOMTOT . 'QUAF' ;
'ELIMINATION' (MDOMTOT ET MDOM1) 0.0001 ;
'ELIMINATION' (MDOMTOT ET MDOM2) 0.0001 ;
'SI' GRAPH;
'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET' ('DOMA' $DOMTOT 'FACEL')
'ET' P10) 'TITRE' 'Domaine et FACEL';
'FINSI' ;
*
**** Proprietes de gaz
*
*
*** GAZ: H_2, O_2
*
* CP, CV en J/Kg/K @ T = 3000
*
PGAZ = 'TABLE' ;
PGAZ . 'CP' = 'TABLE' ;
PGAZ . 'CP' . 'H2 ' = .18729066D+05 ;
PGAZ . 'CP' . 'O2 ' = .11886820D+04 ;
PGAZ . 'CV' = 'TABLE' ;
PGAZ . 'CV' . 'H2 ' = .14571861D+05 ;
PGAZ . 'CV' . 'O2 ' = .92885670D+03 ;
*
**** Especes qui sont dans les equations d'Euler
*
PGAZ . 'ESPEULE' = 'MOTS' 'H2 ';
*
**** Espece qui n'y est pas
*
PGAZ . 'ESPNEULE' = 'O2 ';
rog = 1.0 ;
rod = 2.0;
RN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rog ;
RN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rod ;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2);
pg = 3.0 ;
pd = 4.0 ;
PN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' pg ;
PN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' pd ;
PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2);
ung = 5.0 ;
utg = 6.0 ;
uzg = 6.5 ;
und = 7.0 ;
utd = 8.0 ;
uzd = 8.5 ;
VN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' (ung utg uzg) ;
VN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' (und utd uzd) ;
VN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (VN1 'ET' VN2) ;
* YN
YG = 'PROG' 0.2 ;
YD = 'PROG' 0.0 ;
'REPETER' BL1 ('DIME' (PGAZ . 'ESPEULE') );
CELCAR = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE');
YCELG = 'EXTRAIRE' &BL1 YG ;
YCELD = 'EXTRAIRE' &BL1 YD ;
YN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE'
'COMP' CELCAR ycelg ;
YN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE'
'COMP' CELCAR yceld ;
'SI' (&BL1 'EGA' 1);
YN = YN1 'ET' YN2;
'SINON' ;
YN = YN 'ET' YN1 'ET' YN2 ;
'FINSI' ;
'FIN' BL1 ;
* GAMMA
*
*** gamg
*
YTOT = 0.0 ;
CPTOT = 0.0 ;
CVTOT = 0.0 ;
'REPETER' BL1 ('DIME' (PGAZ . 'ESPEULE') );
CELCAR = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE');
CPCEL = PGAZ . 'CP' . CELCAR ;
CVCEL = PGAZ . 'CV' . CELCAR ;
YCEL = 'EXTRAIRE' &BL1 yg ;
YTOT = YTOT '+' YCEL ;
CPTOT = CPTOT '+' (YCEL '*' CPCEL) ;
CVTOT = CVTOT '+' (YCEL '*' CVCEL) ;
'FIN' BL1 ;
CELCAR = (PGAZ . 'ESPNEULE') ;
CPCEL = PGAZ . 'CP' . CELCAR ;
CVCEL = PGAZ . 'CV' . CELCAR ;
YCEL = 1.0D0 '-' YTOT;
CPTOT = CPTOT '+' (YCEL '*' CPCEL) ;
CVTOT = CVTOT '+' (YCEL '*' CVCEL) ;
gamg = CPTOT '/' CVTOT ;
*
*** gamd
*
YTOT = 0.0 ;
CPTOT = 0.0 ;
CVTOT = 0.0 ;
'REPETER' BL1 ('DIME' (PGAZ . 'ESPEULE') );
CELCAR = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE');
CPCEL = PGAZ . 'CP' . CELCAR ;
CVCEL = PGAZ . 'CV' . CELCAR ;
YCEL = 'EXTRAIRE' &BL1 yd ;
YTOT = YTOT '+' YCEL ;
CPTOT = CPTOT '+' (YCEL '*' CPCEL) ;
CVTOT = CVTOT '+' (YCEL '*' CVCEL) ;
'FIN' BL1 ;
CELCAR = (PGAZ . 'ESPNEULE') ;
CPCEL = PGAZ . 'CP' . CELCAR ;
CVCEL = PGAZ . 'CV' . CELCAR ;
YCEL = 1.0D0 '-' YTOT;
CPTOT = CPTOT '+' (YCEL '*' CPCEL) ;
CVTOT = CVTOT '+' (YCEL '*' CVCEL) ;
gamd = CPTOT '/' CVTOT ;
GAM1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' gamg ;
GAM2 = 'KCHT '$DOM2 'SCAL' 'CENTRE' gamd ;
GAMN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (GAM1 'ET' GAM2);
*
*** Etats gauche et droit
*
ETATG = YG 'ET' ('PROG' rog ung pg gamg) ;
ETATD = YD 'ET' ('PROG' rod und pd gamd) ;
*
*** TEST1: premier ordre en espace, premier ordre en temps
*
IE = 1;
IT = 1;
RNF VNF PNF YNF GAMF = 'PRET' 'PERFMULT' IE IT $DOMTOT
RN VN PN YN GAMN ;
*********************************************************
*** Control des etats sur la surface qui contient P1 ****
*********************************************************
P10 = 1.0 0.5 0.5;
P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P10;
GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
GEOP2 = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
ROGEOP1 = 'REDU' RNF GEOP1;
VGEOP1 = 'REDU' VNF GEOP1;
PGEOP1 = 'REDU' PNF GEOP1;
YGEOP1 = 'REDU' YNF GEOP1;
GAMGEOP1 = 'REDU' GAMF GEOP1;
REFGEOP1 = 'REDU' VNF GEOP2;
rofag = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 1;
rofad = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 3;
unfag = 'EXTRAIRE' VGEOP1 'UN ' 1 1 1;
unfad = 'EXTRAIRE' VGEOP1 'UN ' 1 1 3;
utfag = 'EXTRAIRE' VGEOP1 'UT ' 1 1 1;
utfad = 'EXTRAIRE' VGEOP1 'UT ' 1 1 3;
pfag = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 1;
pfad = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 3;
gamfag = 'EXTRAIRE' GAMGEOP1 'SCAL' 1 1 1;
gamfad = 'EXTRAIRE' GAMGEOP1 'SCAL' 1 1 3;
YFAG = 'PROG' ;
YFAD = 'PROG' ;
'REPETER' BL1 ('DIME' PGAZ . 'ESPEULE');
CELLCH = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE');
YCEL = 'EXTRAIRE' YGEOP1 CELLCH 1 1 1;
YFAG = YFAG 'ET' ('PROG' YCEL);
YCEL = 'EXTRAIRE' YGEOP1 CELLCH 1 1 3;
YFAD = YFAD 'ET' ('PROG' YCEL);
'FIN' BL1;
*
**** Orientation de la normal n de castem par raport a la
* notre; t est par consequence
*
NCOS = 'EXTRAIRE' REFGEOP1 'NX' 1 1 1;
NSIN = 'EXTRAIRE' REFGEOP1 'NY' 1 1 1;
ORIENT = 1.0 '/' NCOS;
ORIENT = 'SIGNE' 'FLOTTANT' ORIENT;
*
**** ORIENT = -1 -> Mon etat gauche est son etat droite
*
'SI' (ORIENT > 0);
ERRLIG = YFAG 'ET' ('PROG' rofag (unfag '*' ORIENT) pfag gamfag) ;
ERRLID = YFAD 'ET' ('PROG' rofad (unfad '*' ORIENT) pfad gamfad) ;
'SINON' ;
ERRLID = YFAG 'ET' ('PROG' rofag (unfag '*' ORIENT) pfag gamfag) ;
ERRLIG = YFAD 'ET' ('PROG' rofad (unfad '*' ORIENT) pfad gamfad) ;
'FINSI' ;
ERRO = 'MAXIMUM' ('PROG'
('MAXIMUM' (ETATG '-' ERRLIG) 'ABS')
('MAXIMUM' (ETATD '-' ERRLID) 'ABS')
);
'SI' (ERRO > 1.0D-14)
'ERREUR' 5 ;
'FINSI' ;
*
*** TEST2: deuxieme ordre en espace, deuxieme ordre en temps
*
* Invariance par rotation
*
* Les gradients
rog = 1.00 ;
grong = 1.0D0 ;
grotg = 2.0D0 ;
grozg = 1.5D0 ;
pg = 2.0D0 ;
gpng = 3.0D0 ;
gptg = 4.0D0 ;
gpzg = 5.0D0 ;
ung = 3.0D0 ;
gunng = 5.0D0 ;
guntg = 6.0D0 ;
gunzg = 7.0D0 ;
utg = 4.0D0 ;
gutng = 7.0D0 ;
guttg = 8.0D0 ;
gutzg = 7.0D0 ;
uzg = 4.0D0 ;
guzxg = 7.0D0 ;
guzyg = 8.0D0 ;
guzzg = 7.0D0 ;
gyng = 0.3 ;
gytg = 0.2 ;
gyzg = 0.1 ;
rod = 1.50 ;
grond = 1.5D0 ;
grotd = 2.5D0 ;
grozd = 2.0D0 ;
pd = 2.5D0 ;
gpnd = 3.5D0 ;
gptd = 4.5D0 ;
gpzd = 5.5D0 ;
und = 3.5D0 ;
gunnd = 5.5D0 ;
guntd = 6.5D0 ;
gunzd = 7.0D0 ;
utd = 4.5D0 ;
gutnd = 7.5D0 ;
guttd = 8.5D0 ;
gutzd = 7.0D0 ;
uzd = 4.5D0 ;
guzxd = 7.5D0 ;
guzyd = 8.5D0 ;
guzzd = 7.0D0 ;
gynd = 0.13 ;
gytd = 0.12 ;
gyzd = 0.11 ;
*
*** Rotation
*
ANGLE = 7.0D0;
DANGLE = 85;
ORIG = 0.0D0 0.0D0 0.0D0;
'REPETER' BL1 5;
ANGLE = ANGLE '+' DANGLE ;
uxg uyg = RUO1 ANGLE UNG UTG ;
uxd uyd = RUO1 ANGLE UND UTD ;
groxg groyg = RUO1 ANGLE grong grotg;
groxd groyd = RUO1 ANGLE grond grotd;
gpxg gpyg = RUO1 ANGLE gpng gptg;
gpxd gpyd = RUO1 ANGLE gpnd gptd;
gyxg gyyg = RUO1 ANGLE gyng gytg;
gyxd gyyd = RUO1 ANGLE gynd gytd;
guxxg guxyg guxzg guyxg guyyg guyzg=
RUO2 ANGLE gunng guntg gunzg gutng guttg gutzg;
guxxd guxyd guxzd guyxd guyyd guyzd=
RUO2 ANGLE gunnd guntd gunzd gutnd guttd gutzd;
'MESSAGE' ;
'MESSAGE' (CHAIN 'Angle de rotation= ' ANGLE);
'MESSAGE' ;
DOM1 = DOM10 'TOURNER' ANGLE ORIG (0.0 0.0 1.0);
DOM2 = DOM20 'TOURNER' ANGLE ORIG (0.0 0.0 1.0);
P1 = P10 'TOURNER' ANGLE ORIG (0.0 0.0 1.0);
'ELIMINATION' (DOM1 ET DOM2) 1D-6;
DOMTOT = DOM1 ET DOM2;
$DOMTOT = 'MODELISER' DOMTOT 'EULER';
$DOM1 = 'MODELISER' DOM1 'EULER';
$DOM2 = 'MODELISER' DOM2 'EULER';
TDOMTOT = 'DOMA' $DOMTOT 'VF';
TDOM1 = 'DOMA' $DOM1 'VF';
TDOM2 = 'DOMA' $DOM2 'VF';
MDOM1 = TDOM1 . 'QUAF' ;
MDOM2 = TDOM2 . 'QUAF' ;
MDOMTOT = TDOMTOT . 'QUAF' ;
'ELIMINATION' (MDOMTOT ET MDOM1) 0.0001 ;
'ELIMINATION' (MDOMTOT ET MDOM2) 0.0001 ;
'SI' GRAPH;
'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET' ('DOMA' $DOMTOT 'FACEL')
'ET' P1) 'TITRE' 'Domaine et FACEL';
'FINSI' ;
*
**** Redefinition de P1 dans $DOMTOT 'FACE'
*
P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1;
***********************
**** Les CHPOINTs ****
***********************
GAM1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' gamg;
GAM2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' gamd;
GAMMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (GAM1 'ET' GAM2);
RN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rog;
RN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rod;;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2);
VIT1 = 'KCHT' $DOM1 'VECT' 'CENTRE' (uxg uyg uzg);
VIT2 = 'KCHT' $DOM2 'VECT' 'CENTRE' (uxd uyd uzd);
VIT = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (VIT1 'ET' VIT2);
PN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' pg;
PN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' pd;
PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2);
YN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2 '
('EXTRAIRE' 1 yg) ;
YN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'H2 '
('EXTRAIRE' 1 yd) ;
YN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2 '
(YN1 'ET' YN2);
*
**** On impose les gradients et le limiteurs
*
ALR = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.5D0 ;
ALP = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.125D0 ;
ALV = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' 'P3'
(0.25 0.25D0 0.25) ;
ALY = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1'
0.1 ;
GRADR1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(groxg groyg grozg) ;
GRADR2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(groxd groyd grozd) ;
GRADR = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(GRADR1 'ET' GRADR2);
GRADP1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(gpxg gpyg gpzg) ;
GRADP2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(gpxd gpyd gpzd) ;
GRADP = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(GRADP1 'ET' GRADP2);
GRADY1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(gyxg gyyg gyzg) ;
GRADY2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(gyxd gyyd gyzd) ;
GRADY = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(GRADY1 'ET' GRADY2);
GRADVX1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(guxxg guxyg guxzg);
GRADVX2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' 'P1DZ'
(guxxd guxyd guxzd);
GRADVY1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' 'P2DZ'
(guyxg guyyg guyzg);
GRADVY2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' 'P2DZ'
(guyxd guyyd guyzd);
GRADVZ1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P3DX' 'P3DY' 'P3DZ'
(guzxg guzyg guzzg);
GRADVZ2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P3DX' 'P3DY' 'P3DZ'
(guzxd guzyd guzzd);
GRADVX = GRADVX1 'ET' GRADVX2 ;
GRADVY = GRADVY1 'ET' GRADVY2 ;
GRADVZ = GRADVZ1 'ET' GRADVZ2 ;
GRADV = GRADVX 'ET' GRADVY 'ET' GRADVZ ;
ORDESP = 2;
ORDTEM = 2;
*
*** L = 1
*
DTCFLG = (gamg '*' pg) '/' rog ;
DTCFLG = DTCFLG '**' 0.5 ;
DTCFLG = DTCFLG '+' ung ;
DTCFLG = 1 '/' DTCFLG ;
DTCFLD = (gamd '*' pd) '/' rod ;
DTCFLD = DTCFLD '**' 0.5 ;
DTCFLD = DTCFLD '+' und ;
DTCFLD = 1 '/' DTCFLD ;
DTCFLD = 'MINIMUM' ('PROG' DTCFLD DTCFLG);
ROF VITF PF YF GAMF = 'PRET' 'PERFMULT' ORDESP ORDTEM
$DOMTOT PGAZ RN GRADR ALR
VIT GRADV ALV
PN GRADP ALP
YN GRADY ALY
GAMMA (1D-3 '*' DTCFLD) ;
*********************************************************
*** Control des etats sur la surface qui contient P1 ****
*********************************************************
GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
GEOP2 = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P1;
ROGEOP1 = 'REDU' ROF GEOP1;
VGEOP1 = 'REDU' VITF GEOP1;
PGEOP1 = 'REDU' PF GEOP1;
YGEOP1 = 'REDU' YF GEOP1;
GAMGEOP1 = 'REDU' GAMF GEOP1;
REFGEOP1 = 'REDU' VITF GEOP2;
rofag = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 1;
rofad = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 3;
unfag = 'EXTRAIRE' VGEOP1 'UN ' 1 1 1;
unfad = 'EXTRAIRE' VGEOP1 'UN ' 1 1 3;
utfag = 'EXTRAIRE' VGEOP1 'UT ' 1 1 1;
utfad = 'EXTRAIRE' VGEOP1 'UT ' 1 1 3;
pfag = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 1;
pfad = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 3;
yfag = 'EXTRAIRE' YGEOP1 'H2 ' 1 1 1;
yfad = 'EXTRAIRE' YGEOP1 'H2 ' 1 1 3;
*
**** Orientation de la normal n de castem par raport a la
* notre; t est par consequence
*
NCOS = 'EXTRAIRE' REFGEOP1 'NX' 1 1 1;
NSIN = 'EXTRAIRE' REFGEOP1 'NY' 1 1 1;
'SI' (('ABS' NCOS) > ('ABS' NSIN));
ORIENT = ('COS' ANGLE) '/' NCOS;
'SINON';
ORIENT = ('SIN' ANGLE) '/' NSIN;
'FINSI' ;
ORIENT = 'SIGNE' 'FLOTTANT' ORIENT;
*
**** ORIENT = -1 -> Mon etat gauche est son etat droite
*
'SI' (ORIENT > 0);
ERRLIG = 'PROG' rofag (unfag '*' ORIENT) pfag yfag;
ERRLID = 'PROG' rofad (unfad '*' ORIENT) pfad yfad;
'SINON' ;
ERRLID = 'PROG' rofag (unfag '*' ORIENT) pfag yfag;
ERRLIG = 'PROG' rofad (unfad '*' ORIENT) pfad yfad;
'FINSI' ;
*
**** 'SI' &BL1 > 1, on controle qui rien n'a change!
*
'SI' (&BL1 > 1);
ERRO = 'MAXIMUM' ('PROG'
('MAXIMUM' (ETATG '-' ERRLIG) 'ABS')
('MAXIMUM' (ETATD '-' ERRLID) 'ABS')
);
'SI' (ERRO > 1.0D-14)
'MESSAGE' 'Ordre en espace = 2';
'MESSAGE' 'Ordre en temps = 2';
'ERREUR' 5 ;
'FINSI' ;
'FINSI' ;
ETATG = ERRLIG ;
ETATD = ERRLID ;
'FIN' BL1;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer