* fichier : pret3D1.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 monoespece "thermally perfect" ****
**** Differents cas tests ****
**** ****
**** A. BECCANTINI DRN/DMT/SEMT/TTMF AOUT 1998 ****
**********************************************************************
'OPTION' 'DIME' 3 ;
'OPTION' 'ELEM' CUB8 ;
'OPTION' 'ECHO' 1 ;
'OPTION' 'TRAC' 'PS' ;
*
*** 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 ;
*
**** TEST1: si les limiteurs sont nuls ou les gradients son nuls
* le deuxieme ordre en espace 'et premier ordre
* en temps degenere en premier ordre en espace
*
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 1.0D0;
PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 2.0D0;
VIT = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (4.0 5.0 3.0);
GAMMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 1.5D0;
*
**** Limiteurs nuls, gradients non-nuls
*
GRADR = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P1DX' 'P1DY' 'P1DZ'
(1.0D0 1.0D0 1.0D0);
GRADP = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P1DX' 'P1DY' 'P1DZ'
(1.0D0 1.0D0 1.0D0);
GRADVX = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P1DX' 'P1DY' 'P1DZ'
(1.0D0 1.0D0 1.0D0);
GRADVY = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P2DX' 'P2DY' 'P2DZ'
(1.0D0 1.0D0 1.0D0);
GRADVZ = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P3DX' 'P3DY' 'P3DZ'
(1.0D0 1.0D0 1.0D0);
GRADV = GRADVX 'ET' GRADVY 'ET' GRADVZ;
ALR = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.0D0 ;
ALP = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.0D0 ;
ALV = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' 'P3'
(0.0 0.0D0 0.0D0) ;
*
ORDESP = 1;
ORDTEM = 1;
ROF1 VITF1 PF1 GAMF1 = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN
VIT
PN
GAMMA ;
'SI' GRAPH;
NN='VECT' (('DOMA' $DOMTOT 'TABLE').'MATROT')
1.0 'MX ' 'MY ' 'MZ ' ROUGE;
TT='VECT' (('DOMA' $DOMTOT 'TABLE').'MATROT')
1.0 'UX ' 'UY ' 'UZ ' VERT;
VV='VECT' (('DOMA' $DOMTOT 'TABLE').'MATROT')
1.0 'RX ' 'RY ' 'RZ ' BLEU;
'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET' ('DOMA' $DOMTOT 'FACEL')
'ET' P10) 'TITRE' 'Domaine et FACEL';
'TRACER' (NN 'ET' TT 'ET' VV) DOMTOT;
'FINSI' ;
ORDESP = 2;
ORDTEM = 1;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN GRADR ALR
VIT GRADV ALV
PN GRADP ALP
GAMMA ;
ERRO = 'MAXIMUM' ('PROG'
('MAXIMUM' (ROF '-' ROF1) 'ABS')
('MAXIMUM' (PF '-' PF1) 'ABS')
('MAXIMUM' (VITF '-' VITF1) 'ABS')
) 'ABS' ;
'SI' (ERRO > 1.0D-15);
'ERREUR' 5;
'FINSI' ;
*
**** Limiteurs non-nuls, gradients nuls
*
GRADR = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P1DX' 'P1DY' 'P1DZ'
(0.0D0 0.0D0 0.0D0);
GRADP = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P1DX' 'P1DY' 'P1DZ'
(0.0D0 0.0D0 0.0D0);
GRADVX = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P1DX' 'P1DY' 'P1DZ'
(0.0D0 0.0D0 0.0D0);
GRADVY = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P2DX' 'P2DY' 'P2DZ'
(0.0D0 0.0D0 0.0D0);
GRADVZ = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P3DX' 'P3DY' 'P3DZ'
(0.0D0 0.0D0 0.0D0);
GRADV = GRADVX 'ET' GRADVY 'ET' GRADVZ;
ALR = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.4D0 ;
ALP = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.8D0 ;
ALV = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
'COMP' 'P1' 'P2' 'P3'
(0.9 0.9D0 0.9D0) ;
ORDESP = 1;
ORDTEM = 1;
ROF1 VITF1 PF1 GAMF1 = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN
VIT
PN
GAMMA ;
ORDESP = 2;
ORDTEM = 1;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN GRADR ALR
VIT GRADV ALV
PN GRADP ALP
GAMMA ;
ERRO = 'MAXIMUM' ('PROG'
('MAXIMUM' (ROF '-' ROF1) 'ABS')
('MAXIMUM' (PF '-' PF1) 'ABS')
('MAXIMUM' (VITF '-' VITF1) 'ABS')
) 'ABS' ;
'SI' (ERRO > 1.0D-15);
'ERREUR' 5;
'FINSI' ;
*
**** TEST2 : Deuxieme ordre en espace, premiere ordre en
* temps: invariance par rotation
*
* Il faut remarquer que:
* densite, pression sont des tenseurs d'ordre 0,
* donc leurs gradients sont des tenseurs d'ordre 1
* la vitesse est un tenseur d'ordre 1,
* donc sont gradient est un tenseur d'ordre 2
*
**** On considere un repere solidal avec la surface, (n,t,v),
* et un repere (x,y,z) (en pratique on a toujours v=z)
*
*** Etats gauche et droite
*
gam = 1.4D0;
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 ;
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 ;
*
*** Etats a gauche et a droite du point face P10
*
* ro, un, ut, p
ETATG = 'PROG' 1.25000 3.625 2.1875;
ETATD = 'PROG' 1.12500 2.8125 2.28125;
*
*** 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;
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 ;
*
**** Redefinition de P1 dans $DOMTOT 'FACE'
*
P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1;
***********************
**** Les CHPOINTs ****
***********************
GAMMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam;
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);
*
**** 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.25D0) ;
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);
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 = 1;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN GRADR ALR
VIT GRADV ALV
PN GRADP ALP
GAMMA ;
'SI' GRAPH;
NN='VECT' (('DOMA' $DOMTOT 'TABLE').'MATROT')
1.0 'MX ' 'MY ' 'MZ ' ROUGE;
TT='VECT' (('DOMA' $DOMTOT 'TABLE').'MATROT')
1.0 'UX ' 'UY ' 'UZ ' VERT;
VV='VECT' (('DOMA' $DOMTOT 'TABLE').'MATROT')
1.0 'RX ' 'RY ' 'RZ ' BLEU;
'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET' ('DOMA' $DOMTOT 'FACEL')
'ET' P1) 'TITRE' 'Domaine et FACEL';
'TRACER' (NN 'ET' TT 'ET' VV) DOMTOT;
'FINSI' ;
*********************************************************
*** 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;
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;
uzfag = 'EXTRAIRE' VGEOP1 'UV ' 1 1 1;
uzfad = 'EXTRAIRE' VGEOP1 'UV ' 1 1 3;
pfag = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 1;
pfad = 'EXTRAIRE' PGEOP1 'SCAL' 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 ;
ERRLID = 'PROG' rofad (unfad '*' ORIENT) pfad ;
'SINON' ;
ERRLID = 'PROG' rofag (unfag '*' ORIENT) pfag ;
ERRLIG = 'PROG' rofad (unfad '*' ORIENT) pfad ;
'FINSI' ;
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 = 1';
'MESSAGE' ERRO;
'ERREUR' 5 ;
'FINSI' ;
'FIN' BL1;
*
**** TEST3: deuxieme ordre en espace et en temps
* Invariance par rotation
*
*
**** On considere un repere solidal avec la surface, (n,t),
* et un repere (x,y)
*
*** Etats gauche et droite
*
gam = 1.4D0;
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 = 0.0D0 ;
uzg = 4.0D0 ;
guzxg = 7.0D0 ;
guzyg = 8.0D0 ;
guzzg = 7.0D0 ;
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 ;
*
*** Etats a gauche et a droite du point face P10
*
* ro, un, ut, p
ETATG = 'PROG' 1.25000 3.625 2.1875;
ETATD = 'PROG' 1.12500 2.8125 2.28125;
*
*** 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;
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 ;
FACTOT = 'DOMA' $DOMTOT 'FACEL' ;
*
**** Redefinition de P1 dans $DOMTOT 'FACE'
*
P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1;
***********************
**** Les CHPOINTs ****
***********************
GAMMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam;
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);
*
**** 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.25D0) ;
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);
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 = (gam '*' pg) '/' rog ;
DTCFLG = DTCFLG '**' 0.5 ;
DTCFLG = DTCFLG '+' ung ;
DTCFLG = 1 '/' DTCFLG ;
DTCFLD = (gam '*' pd) '/' rod ;
DTCFLD = DTCFLD '**' 0.5 ;
DTCFLD = DTCFLD '+' und ;
DTCFLD = 1 '/' DTCFLD ;
DTCFLD = 'MINIMUM' ('PROG' DTCFLD DTCFLG);
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN GRADR ALR
VIT GRADV ALV
PN GRADP ALP
GAMMA (1D-3 '*' DTCFLD) ;
'SI' GRAPH;
NN='VECT' (('DOMA' $DOMTOT 'TABLE').'MATROT')
1.0 'MX ' 'MY ' 'MZ ' ROUGE;
TT='VECT' (('DOMA' $DOMTOT 'TABLE').'MATROT')
1.0 'UX ' 'UY ' 'UZ ' VERT;
VV='VECT' (('DOMA' $DOMTOT 'TABLE').'MATROT')
1.0 'RX ' 'RY ' 'RZ ' BLEU;
'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET' ('DOMA' $DOMTOT 'FACEL')
'ET' P1) 'TITRE' 'Domaine et FACEL';
'TRACER' (NN 'ET' TT 'ET' VV) DOMTOT;
'FINSI' ;
*********************************************************
*** 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;
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;
uzfag = 'EXTRAIRE' VGEOP1 'UV ' 1 1 1;
uzfad = 'EXTRAIRE' VGEOP1 'UV ' 1 1 3;
pfag = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 1;
pfad = 'EXTRAIRE' PGEOP1 'SCAL' 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 ;
ERRLID = 'PROG' rofad (unfad '*' ORIENT) pfad ;
'SINON' ;
ERRLID = 'PROG' rofag (unfag '*' ORIENT) pfag ;
ERRLIG = 'PROG' rofad (unfad '*' ORIENT) pfad ;
'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';
'MESSAGE' ERRO ;
'ERREUR' 5 ;
'FINSI' ;
'FINSI' ;
ETATG = ERRLIG ;
ETATD = ERRLID ;
'FIN' BL1;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer