* fichier : pret3D2.dgibi ************************************************************************ * Section : Fluides Darcy ************************************************************************ ********************************************************************** **** 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 = 'SIGN' 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 = 'SIGN' 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' ;