Télécharger pret3D2.dgibi

Retour à la liste

* 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 = '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' ;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales