Télécharger couplage_TH1D_Th3D.dgibi

Retour à la liste

* fichier :  couplage_TH1D_Th3D.dgibi
************************************************************************
************************************************************************
'OPTION' 'ECHO' 0 ;
* From ~/nlin/sources_dev_new/util_proc :
*BEGINPROCEDUR append
************************************************************************
* NOM         : APPEND
* DESCRIPTION : Rajoute :
*               - un entier à un listentier
*               - un réel   à un listreel
*               - un objet (liste, evolution, matrice ou chpoint)
*                 à un indice de table ('MOT' ou 'ENTIER')
*                    * si l'indice n'existe pas
*                    * 'ET' si l'indice existe 
*                
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 10/09/2004, version initiale
* HISTORIQUE : v1, 10/09/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' APPEND ;
'ARGUMENT' tab/'TABLE' ;
'SI' ('EXISTE' tab) ;
   'ARGUMENT' itab/'MOT' ;
   'SI' ('NON' ('EXISTE' itab)) ;
      'ARGUMENT' itab*'ENTIER' ;
   'FINSI' ;      
   lobj = FAUX ;
   'SI' ('NON' lobj) ;
      'ARGUMENT' lr/'LISTREEL' ;
      'SI' ('EXISTE' lr) ;
         obj = lr ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' le/'LISTENTI' ;
      'SI' ('EXISTE' le) ;
         obj = le ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' lev/'EVOLUTION' ;
      'SI' ('EXISTE' lev) ;
         obj = lev ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' lm/'MAILLAGE' ;
      'SI' ('EXISTE' lm) ;
         obj = lm ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' chpo/'CHPOINT' ;
      'SI' ('EXISTE' chpo) ;
         obj = chpo ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' rig/'RIGIDITE' ;
      'SI' ('EXISTE' rig) ;
         obj = rig ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' matk/'MATRIK' ;
      'SI' ('EXISTE' matk) ;
         obj = matk ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;
   'SI' ('NON' lobj) ;
      cherr = 'CHAINE'
       'Il faut fournir un objet liste, evolution, matrice ou chpoint.'
        ; 
      'ERREUR' cherr ;
   'FINSI' ;      
   'SI' ('EXISTE' tab itab) ;
      tab . itab = 'ET' (tab . itab) obj ;
   'SINON' ;
      tab . itab = obj ;
   'FINSI' ;
   'RESPRO' tab ;      
'FINSI' ;   
'ARGUMENT' lenti/'LISTENTI' ;
'ARGUMENT' lreel/'LISTREEL' ;
'SI' ('EXISTE' lenti) ;
   'ARGUMENT' enti*'ENTIER' ;
   lenti = 'ET' lenti ('LECT' enti) ;
   'RESPRO' lenti ;
'FINSI' ;   
'SI' ('EXISTE' lreel) ;
   'ARGUMENT' reel*'FLOTTANT' ;
   lreel = 'ET' lreel ('PROG' reel) ;
   'RESPRO' lreel ;
'FINSI' ;   
*
* End of procedure file APPEND
*
'FINPROC' ;
*ENDPROCEDUR append
*BEGINPROCEDUR errrel
************************************************************************
* NOM         : ERRREL
* DESCRIPTION : Calcul d'une erreur relative
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/04/2003, version initiale
* HISTORIQUE : v1, 23/04/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' ERRREL ;
'ARGUMENT' val*'FLOTTANT' ;
'ARGUMENT' valref*'FLOTTANT' ;
*
'SI' ('<' ('ABS' valref) 1.D-10) ;
   echref = 1.D0 ;
'SINON' ;
   echref = valref ;   
'FINSI' ;
*
errabs = 'ABS' ('/' ('-' val valref) echref);
*
'RESPRO' errabs ;
*
* End of procedure file ERRREL
*
'FINPROC' ;
*ENDPROCEDUR errrel
*BEGINPROCEDUR exmomod
************************************************************************
* NOM         : EXMOMOD
* DESCRIPTION : Extraction d'un mot d'un listmots
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/06/2003, version initiale
* HISTORIQUE : v1, 23/06/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' EXMOMOD ;
'ARGUMENT' lm*'LISTMOTS' i*'ENTIER' ;
j = 'DIME' lm ;
k = '+' (MODULO ('-' i 1) j) 1 ;
lemot = 'EXTRAIRE' lm k ;
* Usage de l'opérateur text pour éviter que lemot
* ne soit interprété comme un opérateur
'RESPRO' 'TEXTE' lemot ;
*
* End of procedure file EXMOMOD
*
'FINPROC' ;
*ENDPROCEDUR exmomod
*BEGINPROCEDUR formar
************************************************************************
* NOM         : FORMAR
* DESCRIPTION : formate un réel de facon courte
*               pratique pour les noms de
*               sauvegarde
*               Exemples :
* 'MESSAGE' ('CHAINE' (formar 2.9e5   1)) ;
* 2.9E5
* 'MESSAGE' ('CHAINE' (formar -2.9e5  1)) ;
* -2.9E5
* 'MESSAGE' ('CHAINE' (formar 2.9e-5  1)) ;
* 2.9E-5
* 'MESSAGE' ('CHAINE' (formar -2.9e-5 1)) ;
* -2.9E-5
* 'MESSAGE' ('CHAINE' (formar 2.9     1)) ;
* 2.9
* 'MESSAGE' ('CHAINE' (formar -2.9    1)) ;
* -2.9
* 'MESSAGE' ('CHAINE' (formar 0       1)) ;
* 0
* 'MESSAGE' ('CHAINE' (formar 0       1)) ;
* 0
* 'MESSAGE' ('CHAINE' (formar 2.9e5   0)) ;
* 3E5
* 'MESSAGE' ('CHAINE' (formar -2.9e5  0)) ;
* -3E5
* 'MESSAGE' ('CHAINE' (formar 2.9e-5  0)) ;
* 3E-5
* 'MESSAGE' ('CHAINE' (formar -2.9e-5 0)) ;
* -3E-5
* 'MESSAGE' ('CHAINE' (formar 2.9     0)) ;
* 3
* 'MESSAGE' ('CHAINE' (formar -2.9    0)) ;
* -3
* 'MESSAGE' ('CHAINE' (formar 0       0)) ;
* 0
* 'MESSAGE' ('CHAINE' (formar 0       0)) ;
* 0
*               
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 18/02/2003, version initiale
* HISTORIQUE : v1, 18/02/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' FORMAR ;
'ARGUMENT' fl*'FLOTTANT' ;
'ARGUMENT' vir/'ENTIER  ' ;
'SI' ('NON' ('EXISTE' vir)) ;
   vir = 1 ;
'SINON' ;   
   'SI' ('<' vir 0) ;
      'ERREUR' 'fournir un entier positif' ;
   'FINSI' ;
'FINSI' ;         
'SI' ('<' ('ABS' fl) 10.D-100) ;
   chfl = 'CHAINE' '0' ;
'SINON' ;
*! sans le 1.D-10, ca ne fonctionne pas
*! qd on entre pile poil une puissance de 10
   lfl = LOG10 ('ABS' fl) ;
*   lfl = '+' (LOG10 ('ABS' fl)) 1.D-10  ;
   slfl = 'SIGNE' ('ENTIER' lfl) ;
   'SI' ('EGA' slfl 1) ;
      elfl = 'ENTIER' lfl ;
   'SINON' ;
      elfl = '-' ('ENTIER' lfl) 1 ;
   'FINSI' ;
   man = '/' fl ('**' 10.D0 elfl) ;
*
* Une verrue pour des histoires de précision...
*   
   'SI' ('EGA' man 10.D0 ('**' 10.D0 ('*' vir -1.D0))) ;
      man  = '/' man 10.D0 ;
      elfl = '+' elfl 1 ;
   'FINSI' ;
*   
   sman = 'SIGNE' man ;
   'SI' ('EGA' sman 1) ;   
      fman = 'CHAINE' '(F' ('+' vir 2) '.0' vir ')' ;
   'SINON' ;
      fman = 'CHAINE' '(F' ('+' vir 3) '.0' vir ')' ;
   'FINSI' ;
   'SI' ('NEG' vir 0) ;   
      'SI' ('NEG' elfl 0) ;
         chfl = 'CHAINE' 'FORMAT' fman man 'E' elfl ;
      'SINON' ;
         chfl = 'CHAINE' 'FORMAT' fman man ;
      'FINSI' ;
   'SINON' ;
      man2 = 'ENTIER' ('+' man ('*' 0.5D0 sman)) ;
      'SI' ('NEG' elfl 0) ;
         chfl = 'CHAINE' man2 'E' elfl ;
      'SINON' ;
         chfl = 'CHAINE' man2 ;
      'FINSI' ;
   'FINSI' ;            
'FINSI' ;   
'RESPRO' chfl ;
*
* End of procedure file FORMAR
*
'FINPROC' ;
*ENDPROCEDUR formar
*BEGINPROCEDUR log10
************************************************************************
* NOM         : LOG10
* DESCRIPTION : Log_10
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 18/02/2003, version initiale
* HISTORIQUE : v1, 18/02/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' LOG10 ;
'ARGUMENT' fl/'FLOTTANT' ;
'ARGUMENT' lr/'LISTREEL' ;
'ARGUMENT' cp/'CHPOINT ' ;
'ARGUMENT' cm/'MCHAML  ' ;
'SI' ('EXISTE' fl) ;
   'RESPRO' ('/' ('LOG' fl) ('LOG' 10.D0)) ;
'FINSI' ;
'SI' ('EXISTE' lr) ;
   'RESPRO' ('/' ('LOG' lr) ('LOG' 10.D0)) ;
'FINSI' ;
'SI' ('EXISTE' cp) ;
   'RESPRO' ('/' ('LOG' cp) ('LOG' 10.D0)) ;
'FINSI' ;
'SI' ('EXISTE' cm) ;
   'RESPRO' ('/' ('LOG' cm) ('LOG' 10.D0)) ;
'FINSI' ;
*
* End of procedure file LOG10
*
'FINPROC' ;
*ENDPROCEDUR log10
*BEGINPROCEDUR modulo
************************************************************************
* NOM         : MODULO
* DESCRIPTION : Calcule un entier modulo un autre...
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 15/10/2002, version initiale
* HISTORIQUE : v1, 15/10/2002, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' MODULO ;
'ARGUMENT' i*'ENTIER' j*'ENTIER' ;
'SI' ('EGA' j 0) ;
  'MESSAGE' 'Impossible de faire modulo 0' ;
  'ERREUR' 5 ;
'SINON' ;
  k=i '/' j ;
  mod=i '-' ( k '*'j ) ;
  'RESPRO' mod ;
'FINSI' ;
*
* End of procedure file MODULO
*
'FINPROC' ;
*ENDPROCEDUR modulo
* From . :
*BEGINPROCEDUR calcoeil
************************************************************************
* NOM         : CALCOEIL
* DESCRIPTION : Calcul un oeil absolu convenable pour les tracés
*               en donnant un oeil relatif et un maillage
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 16/11/2004, version initiale
* HISTORIQUE : v1, 16/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' CALCOEIL ;
'ARGUMENT' mt*'MAILLAGE' ;
'ARGUMENT' foeil*'POINT' ;
   foeilx foeily foeilz = 'COORDONNEE' foeil ;
   xmt ymt zmt = 'COORDONNEE' mt ;
   maxmt = 'MAXIMUM' xmt ;
   maymt = 'MAXIMUM' ymt ;
   mazmt = 'MAXIMUM' zmt ;
   dx = '-' ('MAXIMUM' xmt) ('MINIMUM' xmt) ;
   dy = '-' ('MAXIMUM' ymt) ('MINIMUM' ymt) ;
   dxy = 'MAXIMUM' ('PROG' dx dy) ;
   dz = '-' ('MAXIMUM' zmt) ('MINIMUM' zmt) ;
   xoeil = '+' maxmt ('*' dxy foeilx) ;
   yoeil = '+' maymt ('*' dxy foeily) ;
   zoeil = '+' mazmt ('*' dz foeilz) ;
   oeil = xoeil yoeil zoeil ;
'RESPRO' oeil ;
*
* End of procedure file CALCOEIL
*
'FINPROC' ;
*ENDPROCEDUR calcoeil
*BEGINPROCEDUR compute
************************************************************************
* NOM         : COMPUTE
* DESCRIPTION : Do the computations of the stationnary fields 
*               in the reactor
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 09/11/2004, version initiale
* HISTORIQUE : v1, 09/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' COMPUTE ;
*
'ARGUMENT' tabdat*'TABLE' ;
*
* First store all the parameters to make tabres self-contained
*
tabres = 'TABLE' 'INCO' ;
tabres . 'tabdat' = tabdat ;
*
* Second, retrieve parameters from tables
*
tabgeo = tabdat . 'tabgeo' ;
tabphy = tabdat . 'tabphy' ;
tabnum = tabdat . 'tabnum' ;
*
H       = tabgeo . 'H'       ;
nh      = tabgeo . 'nh'      ;
$mtcbst = tabgeo . '$mtcbst' ;
$mtcb   = tabgeo . '$mtcb'   ;
$mthe   = tabgeo . '$mthe'   ;
$sthe   = tabgeo . '$sthe'   ;
$mtif   = tabgeo . '$mtif'   ;
*'SI' logb ;
*   $mtifb   = tabgeo . '$mtifb'   ;
*'FINSI' ;   
mcf     = tabgeo . 'mcf'     ;
mfc     = tabgeo . 'mfc'     ;
*
Dh     = tabphy . 'Dh'    ;
chblock= tabphy . 'chblock' ;
*vhe    = tabphy . 'vhe'   ;
qhe    = tabphy . 'qhe'   ;
prhe   = tabphy . 'prhe'  ;
theek  = tabphy . 'theek' ;
pvol   = tabphy . 'pvol'  ;
pvar   = tabphy . 'pvar'  ;
gamma  = tabphy . 'gamma' ;
irrad  = tabphy . 'irrad' ;
pcb    = tabphy . 'pcb'   ;
*
maxit  = tabnum . 'maxit'  ;
rescvg = tabnum . 'rescvg' ;
solvit = tabnum . 'solvit' ;
omcbst = tabnum . 'omcbst' ;
omhe   = tabnum . 'omhe'   ;
*
* Initialization of the tables for the thermal
* computation in the fuel+structure (rvcbst)
* and in the helium (rvhe)
*
rvcbst = 'EQEX' 'ITMA' 1 'NITER' 1 'OMEGA' omcbst 'FIDT' 1000 ;
rvcbst = 'EQEX' rvcbst
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtcbst 'OPER' 'LAPN' 'LCBST' 'INCO' 'TCBS' 
         ;
rvcbst = 'EQEX' rvcbst
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtcb 'OPER' 'FIMP' 'PCB' 'INCO' 'TCBS' 
         ;
rvcbst = 'EQEX' rvcbst
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtif 'OPER' 'ECHI' 'H' 'T0' 'INCO' 'TCBS' 'TCBS' 
         ;
*
* Iterative solver for the fuel+structure part
*
* 1 : direct solver, more robust, slower
* 3 : iterative solver, less robust, faster
'SI' solvit ;
   rvcbst . 'METHINV' . 'TYPINV' = 3 ;
'SINON' ;   
   rvcbst . 'METHINV' . 'TYPINV' = 1 ;
'FINSI' ;   
*         
rvhe = 'EQEX' 'ITMA' 1 'NITER' 1 'OMEGA' omhe 'FIDT' 1000 ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL' 'CENTREE' 'NOCONS'
         'ZONE' $mthe 'OPER' 'KONV' 'RCPHE' 'UHE' 'LHE' 'INCO' 'THE'
         ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mthe 'OPER' 'LAPN' 'LHE' 'INCO' 'THE'
         ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mthe 'OPER' 'FIMP' 'PHE' 'INCO' 'THE' 
         ;
rvhe = 'EQEX' rvhe
         'CLIM' 'THE' 'TIMP' ('DOMA' $sthe 'MAILLAGE') theek ;
rvcbst . 'INCO' = tabres ;
rvhe . 'INCO' = tabres ;
tabres . 'TCBS' = 'KCHT'  $mtcbst 'SCAL' 'SOMMET' theek ;
lupu lsic lcb okcb = PRPCB irrad theek ;
tabres . 'LCBST' = 'KCHT' $mtcbst 'SCAL' 'CENTRE' lsic
                    ('KCHT' $mtcb 'SCAL' 'CENTRE' lcb) ;
tabres . 'PCB' = pcb ;
tabres . 'T0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
tabres . 'F0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
tabres . 'H'  = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
tabres . 'THE' = 'KCHT' $mthe 'SCAL' 'SOMMET' theek ;
tabres . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' 0.D0 ;
thep  = 'COPIER' (tabres . 'THE') ;
tcbsp = 'COPIER' (tabres . 'TCBS') ;
*
* Initialization of the physical properties of helium
*
rhe muhe lhe cphe okhe = PRPHE prhe theek ;
*qhe = '*' rhe vhe ;
qheb  = '*' qhe chblock ;
qhebi = '*' mcf ('NOEL' $mthe qheb) ;
* Ce KCHT est indispensable car '*' supprimme les valeurs nulles...
qhebi = 'KCHT' $mtif 'SCAL' 'CENTRE' qhebi ;
uhe  = 'NOMC' 'UZ' ('/' qheb rhe) 'NATURE' 'DISCRET' ;
tabres . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhe cphe) ;
tabres . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhe ;
*tabres . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' (0.D0 0.D0 uhe) ;
*tabres . 'QHE' = qhe ;
tabres . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' uhe ;
tabres . 'QHE' = qheb ;
***************************************
*
* COMPUTATIONAL LOOP
*
***************************************
* Number of non-linearity iterations
maxiter = vrai ;
*
'REPETER' bouc maxit ;
*   1) Compute a mean temperature in helium
   tmoyhe = 'NOEL' $mthe (tabres . 'THE') ;
      'MESSAGE' 'Maximum mean helium temperature (degrees C) : '
             ('-' ('MAXIMUM' tmoyhe) 273.15D0) ;
*   2) Transfer this temperature on the interfaces
   tmtif  = '*' mcf tmoyhe ;
* Ce KCHT est indispensable car '*' supprimme les valeurs nulles...
   tabres . 'T0' = 'KCHT' $mtif 'SCAL' 'CENTRE' tmtif ;
*   3) Compute the exchange coefficient at the interface
*      and the physical properties of the fuel and the
*      structure with the temperature from the
*      previous iteration
* Mixing temperature
   tm = tabres . 'T0' ;
* Wall temperature
   tp = 'NOEL' $mtif (tabres . 'TCBS') ;
   t0 = '*' ('+' tm tp) 0.5D0 ;
*   t0 = tabres . 'T0' ;
   rho0 mu0 l0 cp0 okhe = PRPHE prhe t0 ;
   pran0 = '*' ('*' mu0 cp0) ('INVERSE' l0) ;
*   reyn0 = '*' ('*' qhe Dh) ('INVERSE' mu0) ;
   reyn0 = '*' ('*' qhebi Dh) ('INVERSE' mu0) ;
*  Colburn correlation
   nuss0 = 0.023 '*' ('**' reyn0 0.8D0)
                 '*' ('**' pran0 ('/' 1.D0 3.D0)) ;
   h0 = '/' ('*' nuss0 l0) Dh ;
* Tests
*   asp = '/' H Dh ;
   minasp  = 60.  ;
   minreyn = 1.D4 ;
   minpran = 0.65 ;
   maxpran = 160. ;
   asp = '/' H ('MAXIMUM' Dh) ;
   mireyn = 'MINIMUM' reyn0 ;
   mareyn = 'MAXIMUM' reyn0 ;
   okreyn = '>' mireyn minreyn ;
   mipran = 'MINIMUM' pran0 ;
   mapran = 'MAXIMUM' pran0 ;
*   okpr   = 'ET' ('>' mipran 0.7) ('<' mapran 160.) ;
   okipr =  '>' mipran minpran ;
   okapr =  '<' mapran maxpran ; 
   okpr   = 'ET' okipr okapr ;
   okasp  = '>' asp minasp ;
   oknuss = okreyn 'ET' okpr 'ET' okasp ;
   'SI' ('NON' oknuss) ;
      cherr = 'CHAINE' '!! Correlation de Colburn non valide' ;
      'MESSAGE' cherr ;
      'SI' ('NON' okreyn) ;
         cherr = 'CHAINE' ' min Re=' (formar mireyn 2)
                          ' (should be '  (formar minreyn 2) ')' ;
         'MESSAGE' cherr ;
         cherr = 'CHAINE' '  (Note : max Re=' (formar mareyn 2) ')' ;
        'MESSAGE' cherr ;
      'FINSI' ;          
      'SI' ('NON' okipr) ;         
         cherr = 'CHAINE' ' min Pr=' (formar mipran 2)
                          ' (should be '  (formar minpran 2) ')' ;
         'MESSAGE' cherr ;
      'FINSI' ;
      'SI' ('NON' okapr) ;                                   
         cherr = 'CHAINE' ' max Pr=' (formar mapran 2)
                          ' (should be '  (formar maxpran 2) ')' ;
         'MESSAGE' cherr ;
      'FINSI' ;
      'SI' ('NON' okasp) ;                                   
         cherr = 'CHAINE' ' min H/Dh=' (formar asp 2)
                          ' (should be '  (formar minasp 2) ')' ;
         'MESSAGE' cherr ;
      'FINSI' ;
   'FINSI' ;      
*   tabres . 'H' =  'KCHT' $mtif 'SCAL' 'CENTRE' ('*' h0 chblock) ;
   tabres . 'H' =  'KCHT' $mtif 'SCAL' 'CENTRE' h0 ;
   tabres . 'nuss0' = 'KCHT' $mtif 'SCAL' 'CENTRE' nuss0 ;
   tabres . 'reyn0' = 'KCHT' $mtif 'SCAL' 'CENTRE' reyn0 ;
   tabres . 'pran0' = 'KCHT' $mtif 'SCAL' 'CENTRE' pran0 ;
   tabres . 'l0' = 'KCHT' $mtif 'SCAL' 'CENTRE' l0 ;
   tabres . 'Dh' = 'KCHT' $mtif 'SCAL' 'CENTRE' Dh ;
* Physical properties of fuel and structure
   tcbck = 'NOEL' $mtcbst (tabres . 'TCBS') ;
   lupu lsic lcb okcb = PRPCB irrad tcbck ;
   tabres . 'LCBST' = 'KCHT' $mtcbst 'SCAL' 'CENTRE' lsic
                       ('KCHT' $mtcb 'SCAL' 'CENTRE' lcb) ;
*   4) Thermal computation in fuel + structure
   EXEC rvcbst ;
*   5) Flux computation on the surface of the fuel
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif 'OPER' 'MDIA' 'H' 'INCO' 'TCBS' 'F0' ;
   rvmdia1 . 'INCO' = tabres ;                   
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif 'OPER' 'MDIA' 'H' 'INCO' 'T0'  'F0' ;
   rvmdia2 . 'INCO' = tabres ;                   
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   f0 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCBS' (tabres . 'TCBS')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'T0'  (tabres . 'T0'))))
                ;
      'MESSAGE' 'Integrated flux on the surface of the fuel (W) : '
             ('MAXIMUM' ('RESULT' f0)) ;
*   6) Transfer it to a volume power in the helium
   tabres . 'F0' = f0 ;
   phe = '*' mfc f0 ;
   tabres . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' phe ;
*
*   7) Compute the physical properties of helium with temperature and
*      pressure from the previous iteration
*
   theck = 'NOEL' $mthe (tabres . 'THE') ;
   thesk = tabres . 'THE' ;
   rhec muhec lhec cphec okhe = PRPHE prhe theck ;
   rhes muhes lhes cphes okhe = PRPHE prhe thesk ;
*   uhes = 'NOMC' 'UZ' ('*' qhe ('INVERSE' rhes)) ;
   uhes = 'NOMC' 'UZ' ('*' qheb ('INVERSE' rhes)) ;
   tabres . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhec cphec) ;
   tabres . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhec ;
   tabres . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' uhes ;
*   8) Thermal computation in the helium
   EXEC rvhe ;
*
* Convergence check
*
   then  = 'COPIER' (tabres . 'THE') ;
   tcbsn = 'COPIER' (tabres . 'TCBS') ;
   errhe  = '/' ('MAXIMUM'  ('-' then thep) 'ABS')
               ('MAXIMUM'  thep 'ABS') ;
   errcbs = '/' ('MAXIMUM'  ('-' tcbsn tcbsp) 'ABS')
               ('MAXIMUM'  tcbsp 'ABS') ;
   'MESSAGE' ('CHAINE' 'Iteration ' &bouc) ;
   'MESSAGE' ('CHAINE' '   errhe  = ' errhe) ;
   'MESSAGE' ('CHAINE' '   errcbs = ' errcbs) ;
   thep = then ;
   tcbsp = tcbsn ;
   testcvg = 'ET' (errhe '<' rescvg)  (errcbs '<' rescvg) ;
*
* Plots 
*
*   explores tabres ;
   'SI' testcvg ;
      maxiter = faux ;
      'QUITTER' bouc ;
    'FINSI' ;      
'FIN' bouc ;
'SI' ('EGA' maxiter vrai) ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'The computation has not converged') ;
   'ERREUR' 'Check compute.procedur !' ;
'FINSI' ;
'RESPRO' tabres ;
*
* End of procedure file COMPUTE
*
'FINPROC' ;
*ENDPROCEDUR compute
*BEGINPROCEDUR cpowdis
************************************************************************
* NOM         : CPOWDIS
* DESCRIPTION : Compute cos (z) power distribution
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 09/11/2004, version initiale
* HISTORIQUE : v1, 09/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' CPOWDIS ;
'ARGUMENT' tabgeo*'TABLE' ;
'ARGUMENT' pvol*'FLOTTANT' ;
'ARGUMENT' gamma*'FLOTTANT' ;
'ARGUMENT' pvar*'LOGIQUE' ;
*
$mtcbst = tabgeo . '$mtcbst' ;
$mtcb   = tabgeo . '$mtcb'   ;
$mthe   = tabgeo . '$mthe'   ;
*
vcbst = 'MAXIMUM' ('RESULT' ('DOMA' $mtcbst 'VOLUME')) ;
vcb   = 'MAXIMUM' ('RESULT' ('DOMA' $mtcb 'VOLUME'))   ;
vhe   = 'MAXIMUM' ('RESULT' ('DOMA' $mthe 'VOLUME'))   ;
vcoeur = '+' vcbst vhe ;
*
pmoy = pvol '*' ('/' vcoeur vcb) ;
'SI' ('NON' pvar) ;
   pcb = 'KCHT' $mtcb 'SCAL' 'CENTRE' pmoy ;
'SINON' ;   
   raddeg = '/' 180.D0 PI ;
   zcb = 'COORDONNEE' 3 ('DOMA' $mtcb 'CENTRE') ;
   zmoy = '/' H 2.D0 ;
   zcb = '+' zcb ('*' zmoy -1.D0) ;
   pmax = '*' pmoy ('/' gamma ('SIN' ('*' raddeg gamma))) ;
   'MESSAGE' ('CHAINE' 'gamma=' gamma) ;
   'MESSAGE' ('CHAINE' 'pmax=' pmax) ;
   fccos = '/' (2.D0 '*' gamma '*' raddeg) H ;
   pcb1 = '*' ('COS' ('*' zcb fccos)) pmax ;
   chvol = 'DOMA' $mtcb 'VOLUME' ;
   puis1 = 'MAXIMUM' ('RESULT' ('*' chvol pmoy)) ;
   puis2 = 'MAXIMUM' ('RESULT' ('*' chvol pcb1)) ;
  'SI' ('NON' ('EGA' puis1 puis2 ('*' 0.1 ('ABS' puis2)))) ;
      'MESSAGE' ('CHAINE' 'puis1 = ' puis1 ) ;
      'MESSAGE' ('CHAINE' 'puis2 = ' puis2 ) ;
      'ERREUR' 'Les deux puissances devraient etre a peu pres egales' ;
  'FINSI' ;         
* Correction
   facor = '/' puis1 puis2 ;
   pcb2 = '*' pcb1 facor ;
   pcb = 'KCHT' $mtcb 'SCAL' 'CENTRE' pcb2 ;
'FINSI' ;
'RESPRO' pcb ;
*
* End of procedure file CPOWDIS
*
'FINPROC' ;
*ENDPROCEDUR cpowdis
*BEGINPROCEDUR dessevol
************************************************************************
* NOM         : DESSEVOL
* DESCRIPTION : Dessine des évolutions : choisit automatiquement
*               les options, marqueurs, couleurs...
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 16/11/2004, version initiale
* HISTORIQUE : v1, 16/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' DESSEVOL ;
'ARGUMENT' evtot*'EVOLUTION' ;
'ARGUMENT' tabt*'TABLE' ;
'ARGUMENT' tit*'MOT' ;
'ARGUMENT' tix*'MOT' ;
'ARGUMENT' tiy*'MOT' ;
'ARGUMENT' lnclk*'LOGIQUE' ;
'ARGUMENT' nb/'LOGIQUE' ;
*
'SI' ('NON' ('EXISTE' nb)) ;
   nb = FAUX ;
'FINSI' ;   
*
nt  = 'DIME' tabt ;
nev = 'DIME' evtot ;
*
* Attention, dans evtot, il y a une évolution avec des noms de points ?
*
*'SI' ('NEG' nev nt) ;
*   cherr = 'CHAINE' 'Evolution and title table : not same dim.' ;
*   'ERREUR' cherr ;
*'FINSI' ;
*
tev = 'TABLE' ;
tev . 'TITRE' = tabt ;
*
toto = 'TABLE' ;
*
lcoul = 'MOTS' 'TURQ' 'VERT' 'JAUN' 'ROSE' 'ROUG' 'BLEU' ;
lmarq = 'MOTS' 'TRIB' 'TRIA' 'LOSA' 'CARR' 'ETOI' 'PLUS' 'CROI' ;
ltirr = 'MOTS' 'TIRR' 'TIRC' 'TIRL' 'TIRM' ;
*
'SI' nb ;
   ev2 = evtot ;
'SINON' ;
   icou = 0 ;
   'REPETER' iev nev ;
      ii = &iev ;
      evi = 'EXTRAIRE' evtot 'COUR' ii ;
      'SI' ('NEG' ('TYPE' ('EXTRAIRE' evi 'ORDO')) 'LISTMOTS') ;
         icou = '+' icou 1 ;
      'FINSI' ;         
*      ii2 = '/' ('+' ii 1) 2 ;
*      ci  = EXMOMOD lcoul ii2 ;
*      ci  = EXMOMOD lcoul ii ;
      ci  = EXMOMOD lcoul icou ;
      APPEND toto 'EVOLUTION' ('COULEUR' evi ci) ;
   'FIN' iev ;
   ev2 = toto . 'EVOLUTION' ;
'FINSI' ;
*      
'REPETER' iev nev ;
   ii = &iev ;
   mi = EXMOMOD lmarq ii ;
   ti = EXMOMOD ltirr ii ;
   'SI' nb ;
      tev . ii = 'CHAINE' 'MARQ ' mi ' ' ti ;
   'SINON' ;      
      tev . ii = 'CHAINE' 'MARQ ' mi ;
   'FINSI' ;      
'FIN' iev ;   
*
'SI' lnclk ;
   'DESSIN' ev2 'TITR' tit 'TITX' tix 'TITY' tiy 'LEGE' 'MIMA' tev
   'NCLK' ;
'SINON' ;   
   'DESSIN' ev2 'TITR' tit 'TITX' tix 'TITY' tiy 'LEGE' 'MIMA' tev ;
'FINSI' ;   
*
* End of procedure file DESSEVOL
*
'FINPROC' ;
*ENDPROCEDUR dessevol
*BEGINPROCEDUR explodat
************************************************************************
* NOM         : EXPLODAT
* DESCRIPTION : Explore the data interactively or not.
*
* - first argument (mandatory) : the tables with data
* - second argument (optional) : VRAI (default) interactive use
*                                     (mouse clicks)
*                                FAUX batch use
* - third argument (mandatory if 2ndarg = FAUX) : a list of integer or
*                                                  an integer  
*              an integer corresponds to a particular view
*   View 1 : the base mesh
*   View 2 : the mesh
*   View 3 : the blocked channels
*   View 4 : the power distribution
*   View 5 : fuel and structure thermal conductivity (T)
*   View 6 : Helium density
*   View 7 : Helium thermal conductivity
*   View 8 : Helium heat capacity
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 17/11/2004, version initiale
* HISTORIQUE : v1, 17/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' EXPLODAT ;
'ARGUMENT' tabdat*'TABLE';
'ARGUMENT' interact/'LOGIQUE' ;
'SI' ('NON' ('EXISTE' interact)) ;
   interact = VRAI ;
'FINSI' ;   
*
'SI' ('NON' interact) ;
   'ARGUMENT' liview/'LISTENTI' ;
   'SI' ('NON' ('EXISTE' liview)) ;
      'ARGUMENT' iview*'ENTIER' ;
      liview = 'LECT' iview ;
   'FINSI' ;      
   niview = 'DIME' liview ;
   'SI' ('<' niview 1) ;
      'QUITTER' EXPLODAT ;
   'FINSI' ;       
   iiview = 1 ;
'FINSI' ;
*
* Retrieve parameters from tables
*
tabgeo = tabdat . 'tabgeo' ;
tabphy = tabdat . 'tabphy' ;
tabin  = tabgeo . 'tabin' ;
logst  = 'EXISTE' tabin 'STRUCTURE' ;
*
H       = tabgeo . 'H'       ;
nh      = tabgeo . 'nh'      ;
$mtcbst = tabgeo . '$mtcbst' ;
$mtcb   = tabgeo . '$mtcb'   ;
$mthe   = tabgeo . '$mthe'   ;
$sthe   = tabgeo . '$sthe'   ;
$mtif   = tabgeo . '$mtif'   ;
cfc     = tabgeo . 'cfc'     ;
amtcb   = tabgeo . 'amtcb' ;
amthe   = tabgeo . 'amthe'   ;
* QUAF meshes
_mtcbst = 'DOMA' $mtcbst 'QUAF' ;
_mthe   = 'DOMA' $mthe   'QUAF' ;
mtcbst = 'DOMA' $mtcbst 'MAILLAGE' ;
mtcb   = 'DOMA' $mtcb   'MAILLAGE' ;
mthe   = 'DOMA' $mthe   'MAILLAGE' ;
'SI' logst ;
   mtst   = 'DIFF' mtcbst mtcb ;
'FINSI' ;   
mthe   = 'DOMA' $mthe   'MAILLAGE' ;
mtif   = 'DOMA' $mtif   'MAILLAGE' ;
*
Dh     = tabphy . 'Dh'    ;
chblock= tabphy . 'chblock' ;
*vhe    = tabphy . 'vhe'   ;
qhe    = tabphy . 'qhe'   ;
prhe   = tabphy . 'prhe'  ;
theek  = tabphy . 'theek' ;
pvol   = tabphy . 'pvol'  ;
pvar   = tabphy . 'pvar'  ;
gamma  = tabphy . 'gamma' ;
irrad  = tabphy . 'irrad' ;
pcb    = tabphy . 'pcb'   ;
* compute bottom He speed
rhohe muhe lhe cphe dum = PRPHE prhe theek ;
vhe = '/' qhe rhohe ;
*
* Precomputations
*
* One wants the height (z) not to be more than four times x or y dimension
*   targetz = 4.D0 ;
   targetz = 1.5D0 ;
*   targetz = 0.5D0 ;
   ftour = 0.1D0 ;
   _mtot = 'ET' _mtcbst _mthe ;
   xmt ymt zmt = 'COORDONNEE' _mtot ;
   dx = '-' ('MAXIMUM' xmt) ('MINIMUM' xmt) ;
   dy = '-' ('MAXIMUM' ymt) ('MINIMUM' ymt) ;
   dxy = 'MAXIMUM' ('PROG' dx dy) ;
   dz = '-' ('MAXIMUM' zmt) ('MINIMUM' zmt) ;
   dzdx = '/' dz dx ;
   chzscale = '>' dzdx targetz ;
   'SI' chzscale ;
      zfactor = '/' dzdx targetz ;
      factech =  ('-' ('/' 1.D0 zfactor) 1.D0) ;
      dzmt = 'NOMC' 'UZ' ('*' zmt factech) 'NATURE' 'DISCRET' ;
      orig = 'FORME' ;
      defrm = 'FORME' dzmt ;
   'FINSI' ;      
*
   'SI' interact ;                             
      mail = ('COULEUR' mtcb 'ROUG')
        'ET' ('COULEUR' mthe 'BLEU') ;
      'SI' (logst) ;
         mail = 'ET' mail ('COULEUR' mtst 'BLANC') ;
      'FINSI' ;                 
      npomail = 'NBNO' mail ;
      nelmail = 'NBEL' mail ;
      titmail = 'CHAINE' 'Mesh ' 'NBPO=' npomail
                                ' NBELEM=' nelmail ;
      TRACMESH mail titmail VRAI ;
   'FINSI' ;      
*
* Precomp of physical properties
*
*  Temperature evolutions (in °C !!!)
   Tmin = 500.D0 ;
   Tmax = 1600.D0 ;
   npasT = 10.D0 ;
   pasT = '/' ('-' Tmax Tmin) npasT ;
   prT = 'PROG' Tmin 'PAS' pasT Tmax ;
   prlupu = 'PROG' ;
   prlsic = 'PROG' ;
   prlcb  = 'PROG' ;
   nbc = 'DIME' prT ;
   'REPETER' ibc nbc ;
      T = 'EXTRAIRE' pRT &ibc ;
      T = '+' T 273.15D0 ;
      lupu lsic lcb dum = PRPCB FAUX  T ;
      prlupu  = 'ET' prlupu  ('PROG' lupu) ;
      prlsic  = 'ET' prlsic  ('PROG' lsic) ;
      prlcb   = 'ET' prlcb   ('PROG' lcb) ; 
   'FIN' ibc ;   
   evlupu = 'EVOL' 'MANU' 'T' prT 'Lambda' prlupu ;
   evlsic = 'EVOL' 'MANU' 'T' prT 'Lambda' prlsic ;
   evlcb =  'EVOL' 'MANU' 'T' prT 'Lambda' prlcb ;
* Helium prop
*  Temperature evolutions (in °C !!!)
Tmin = 450.D0 ;
Tmax = 1150.D0 ;
npasT = 10.D0 ;
pasT = '/' ('-' Tmax Tmin) npasT ;
prT    = 'PROG' Tmin 'PAS' pasT Tmax ;
prrho  = 'PROG' ;
prlamb = 'PROG' ;
prcp   = 'PROG' ;
nbc = 'DIME' prT ;
'REPETER' ibc nbc ;
   T = 'EXTRAIRE' pRT &ibc ;
   T = '+' T 273.15D0 ;
   rho mu lamb cp dum = PRPHE prhe T ;
   prrho = 'ET' prrho ('PROG' rho) ;
   prlamb = 'ET' prlamb ('PROG' lamb) ;
   prcp = 'ET' prcp ('PROG' cp) ;
'FIN' ibc ;
evrho  = 'EVOL' 'MANU' 'T' prT 'Rho' prrho ;
evlamb = 'EVOL' 'MANU' 'T' prT 'Lambda' prlamb ;
evcp   = 'EVOL' 'MANU' 'T' prT 'Cp' prcp ;
*
* Arete precomp
*   
*amtcb = 'ARETE' mtcb ;
*
   tcha = 'TABLE' ;
   tcha . 1 = 'BMesh'  ;
   tcha . 2 = 'Mesh'   ;
   tcha . 3 = 'BChan'  ;
   tcha . 4 = 'FuelT2' ;
   tcha . 5 = 'LFuel'  ;
   tcha . 6 = 'RhoHe'  ;
   tcha . 7 = 'LamHe'  ;
   tcha . 8 = 'CpHe'   ;
   ntcha = 'DIME' tcha ;
   'SI' interact ;
      itcha = 1 ;
      noclic = VRAI ;
   'SINON' ;
      itcha = 'EXTRAIRE' liview iiview ;
      noclic = FAUX ;
   'FINSI' ;            
*
'REPETER' bouc2 ;
'SI' interact ;
   'MESSAGE' 'Noclic allows not to click between each view' ;   
   'MESSAGE'
'However, one will deactivate it in order to zoom or save' ;
   noclic = 'CHOI' 'Choose desired options' noclic ;
'FINSI' ;   
*
'REPETER' bouc1;
*
'SI' ('EGA' itcha 1) ;
   bmesh = ('COULEUR' (tabin . 'FUEL') 'ROUG') 'ET' 
           ('COULEUR' (@STBL (tabin . 'CHANNEL')) 'BLEU') ;
   'SI' logst ;
      bmesh = 'ET' bmesh (tabin . 'STRUCTURE') ;
   'FINSI' ;
   oeil = (0. 0. 1000.) ;
   npomail = 'NBNO' bmesh ;
   nelmail = 'NBEL' bmesh ;
   tit = 'CHAINE' 'Base Mesh ' 'NBPO=' npomail
                               ' NBELEM=' nelmail ;
   'SI' noclic ;
      TRAC oeil bmesh 'TITRE' tit 'NCLK' ;
   'SINON' ;      
      TRAC oeil bmesh 'TITRE' tit  ;
   'FINSI' ;      
'FINSI' ;
'SI' ('EGA' itcha 2) ;
   mail = ('COULEUR' mtcb 'ROUG')
     'ET' ('COULEUR' mthe 'BLEU') ;
   'SI' (logst) ;
      mail = 'ET' mail ('COULEUR' mtst 'BLANC') ;
   'FINSI' ;                 
   npomail = 'NBNO' mail ;
   nelmail = 'NBEL' mail ;
   tit = 'CHAINE' 'Mesh ' 'NBPO=' npomail
                          ' NBELEM=' nelmail ;
   TRACMESH mail tit noclic ;
'FINSI' ;
'SI' ('EGA' itcha 3) ;
*   tit = 'CHAINE' 'Blocked channels'  ;
*   tracchml $mtif chblock tit noclic ;
   tit = 'CHAINE' 'Blocked channels'  ;
   tracchpo chblock mthe amthe tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 4) ;
   tit = 'CHAINE' 'Power Distribution (W.m-3)' ;
   tracchml $mtcb pcb amtcb tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 5) ;
   evt = evlupu 'ET' evlsic 'ET' evlcb ;
   tt = 'TABLE' ;
   'SI' irrad ;
      tt . 1 = 'CHAINE' '(U,Pu)C irr.' ;
      tt . 2 = 'CHAINE' 'SiC irr.'     ;
      tt . 3 = 'CHAINE' 'CerCer irr.'  ;
   'SINON' ;      
      tt . 1 = 'CHAINE' '(U,Pu)C' ;
      tt . 2 = 'CHAINE' 'SiC' ;
      tt . 3 = 'CHAINE' 'CerCer' ;
   'FINSI' ;      
   tit = 'CHAINE' 'Lambda(Temperature)' ;
   titx = 'CHAINE' 'T (°C)' ;
   tity = 'CHAINE' 'Lambda (W.m-1.K-1)' ;
   DESSEVOL evt tt tit titx tity noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 6) ;
   tt = 'TABLE' ;
   tt . 1 = 'CHAINE' 'Helium' ;
   tit = 'CHAINE' 'Density(Temperature)'
         ' pres= ' (formar prhe 2) ' bar' ;
   titx = 'CHAINE' 'T (°C)' ;
   tity = 'CHAINE' 'Rho (kg.m-3)' ;
   DESSEVOL evrho tt tit titx tity noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 7) ;
   tt = 'TABLE' ;
   tt . 1 = 'CHAINE' 'Helium' ;
   tit = 'CHAINE' 'Thermal Conductivity (Temperature)'
         ' pres= ' (formar prhe 2) ' bar' ;
   titx = 'CHAINE' 'T (°C)' ;
   tity = 'CHAINE' 'Lambda (W.m-1.K-1)' ;
   DESSEVOL evlamb tt tit titx tity noclic ;
'FINSI' ;   
*
'SI' ('EGA' itcha 8) ;
   tt = 'TABLE' ;
   tt . 1 = 'CHAINE' 'Helium' ;
   tit = 'CHAINE' 'Specific Heat (Temperature)'
         ' pres= ' (formar prhe 2) ' bar' ;
   titx = 'CHAINE' 'T (°C)' ;
   tity = 'CHAINE' 'Cp (J.kg-1.K-1)' ;
   DESSEVOL evcp tt tit titx tity noclic ;
'FINSI' ;   
*
'SI' interact ;
   cha = tit ;
   ret = 'MENU' cha 'Next' 'Previous' 'Options' ;
   'SI' ('EGA' ret 'Options') ; 'QUITTER' bouc1  ; 'FINSI';
   'SI' ('EGA' ret 'Quitter') ; 'QUITTER' bouc2; 'FINSI';
   'SI' ('EGA' ret 'Next') ; itcha = '+' itcha 1 ; 'FINSI' ; 
   'SI' ('EGA' ret 'Previous') ; itcha = '-' itcha 1 ; 'FINSI' ; 
   'SI' (itcha > ntcha); itcha = 1      ; 'FINSI';
   'SI' (itcha < 1); itcha = ntcha      ; 'FINSI';
'SINON' ;
   iiview = '+' iiview 1 ;
   'SI' ('>' iiview niview) ;
      'QUITTER' bouc2 ;
   'SINON' ;
      itcha = 'EXTRAIRE' liview iiview ;
   'FINSI' ;
'FINSI' ;            
*
'FIN' bouc1;
'FIN' bouc2 ;
*
'SI' chzscale ;
   'FORME' orig ;
'FINSI' ;      
*
* End of procedure file EXPLODAT
*
'FINPROC' ;
*ENDPROCEDUR explodat
*BEGINPROCEDUR explores
************************************************************************
* NOM         : EXPLORES
* DESCRIPTION : Explore the Results interactively or not.
*
* EXPLORES does graphical output of results
* - first argument (mandatory) : the table with the fields
* - second argument (optional) : VRAI (default) interactive use
*                                     (mouse clicks)
*                                FAUX batch use
* - third argument (mandatory if 2ndarg = FAUX) : a list of integer or
*                                                 an integer  
*              an integer corresponds to a particular view
*   View 1 : temperature (fuel + He)
*   View 2 : temperature (fuel)
*   View 3 : temperature (fuel : 2 viewpoints)
*   View 4 : exchange coefficient
*   View 5 : temperature (helium)
*   View 6 : vertical speed (helium)
*   View 7 : temperature profiles per channel (helium)
*   View 8 : vertical speed profiles per channel (helium)
*   View 9  : Nusselt number
*   View 10 : Reynolds number
*   View 11 : Prandtl number
*   View 12 : Lambda Helium
*   View 13 : Hydraulic diameter
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 17/11/2004, version initiale
* HISTORIQUE : v1, 17/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' EXPLORES ;
'ARGUMENT' tabres*'TABLE' ;
'ARGUMENT' interact/'LOGIQUE' ;
'SI' ('NON' ('EXISTE' interact)) ;
   interact = VRAI ;
'FINSI' ;   
*
'SI' ('NON' interact) ;
   'ARGUMENT' liview/'LISTENTI' ;
   'SI' ('NON' ('EXISTE' liview)) ;
      'ARGUMENT' iview*'ENTIER' ;
      liview = 'LECT' iview ;
   'FINSI' ;      
   niview = 'DIME' liview ;
   'SI' ('<' niview 1) ;
      'QUITTER' EXPLORES ;
   'FINSI' ;       
   iiview = 1 ;
'FINSI' ;
*
* Retrieve parameters from tables
*
tabdat = tabres . 'tabdat' ;
tabgeo = tabdat . 'tabgeo' ;
tabin  = tabgeo . 'tabin' ;
logst  = 'EXISTE' tabin 'STRUCTURE' ;
*
H       = tabgeo . 'H'       ;
nh      = tabgeo . 'nh'      ;
$mtcbst = tabgeo . '$mtcbst' ;
$mtcb   = tabgeo . '$mtcb'   ;
$mthe   = tabgeo . '$mthe'   ;
$sthe   = tabgeo . '$sthe'   ;
$mtif   = tabgeo . '$mtif'   ;
cfc     = tabgeo . 'cfc'     ;
amtcbst = tabgeo . 'amtcbst' ;
amthe   = tabgeo . 'amthe'   ;
* QUAF meshes
_mtcbst = 'DOMA' $mtcbst 'QUAF' ;
_mthe   = 'DOMA' $mthe   'QUAF' ;
mtcbst = 'DOMA' $mtcbst 'MAILLAGE' ;
mtcb   = 'DOMA' $mtcb   'MAILLAGE' ;
'SI' logst ;
   mtst   = 'DIFF' mtcbst mtcb ;
'FINSI' ;   
mtst   = 'DIFF' mtcbst mtcb ;
mthe   = 'DOMA' $mthe   'MAILLAGE' ;
mtif   = 'DOMA' $mtif   'MAILLAGE' ;
*
* Precomputations
*
* One wants the height (z) not to be more than four times x or y dimension
*   targetz = 4.D0 ;
   targetz = 1.5D0 ;
   ftour = 0.1D0 ;
   _mtot = 'ET' _mtcbst _mthe ;
   xmt ymt zmt = 'COORDONNEE' _mtot ;
   dx = '-' ('MAXIMUM' xmt) ('MINIMUM' xmt) ;
   dy = '-' ('MAXIMUM' ymt) ('MINIMUM' ymt) ;
   dxy = 'MAXIMUM' ('PROG' dx dy) ;
   dz = '-' ('MAXIMUM' zmt) ('MINIMUM' zmt) ;
   dzdx = '/' dz dx ;
   chzscale = '>' dzdx targetz ;
   'SI' chzscale ;
      zfactor = '/' dzdx targetz ;
      factech =  ('-' ('/' 1.D0 zfactor) 1.D0) ;
      dzmt = 'NOMC' 'UZ' ('*' zmt factech) 'NATURE' 'DISCRET' ;
      orig = 'FORME' ;
      defrm = 'FORME' dzmt ;
   'FINSI' ;      
*
* manipulation maillage et champ de temperature
   tcbsc = ('-' (tabres . 'TCBS') 273.15D0) ;
   thec  = ('-' (tabres . 'THE')  273.15D0) ;
   cech  = tabres . 'H' ;
   nuss0 = tabres . 'nuss0' ;
   reyn0 = tabres . 'reyn0' ;
   pran0 = tabres . 'pran0' ;
   l0 = tabres . 'l0' ;
   Dh = tabres . 'Dh' ;
   uz = 'EXCO' 'UZ' (tabres . 'UHE') ;
* tournage 
   mtot = 'ET' mtcbst mthe ;
   xmt ymt zmt = 'COORDONNEE' mtot ;
   mixmt = 'MINIMUM' xmt ;
   miymt = 'MINIMUM' ymt ;
   mizmt = 'MINIMUM' zmt ;
   dx = '-' ('MAXIMUM' xmt) ('MINIMUM' xmt) ;
   dy = '-' ('MAXIMUM' ymt) ('MINIMUM' ymt) ;
   dz = '-' ('MAXIMUM' zmt) ('MINIMUM' zmt) ;
   p1 = ('-' mixmt ('*' dx ftour))
        ('-' miymt ('*' dy ftour))
        mizmt ;
   p2 = ('-' mixmt ('*' dx ftour))
        ('-' miymt ('*' dy ftour))
        ('+' mizmt 1.D0) ;
*   amtcbst  = 'ARETE' mtcbst  ;
*   amthe    = 'ARETE' mthe    ;
   tcbsc2 thec2 mtcbst2 mthe2 amtcbst2 amthe2 =
    'TOURNER' tcbsc thec mtcbst mthe amtcbst amthe 180. p1 p2 ;
   tcbsct  = 'ET' tcbsc tcbsc2 ;
   thect   = 'ET' thec  thec2 ;
   mtcbstt = 'ET' mtcbst mtcbst2 ;
   mthet   = 'ET' mthe   mthe2 ;
   amtcbstt = 'ET' amtcbst amtcbst2 ;
   amthet   = 'ET' amthe   amthe2 ;
*
   'SI' interact ;                             
      mail = ('COULEUR' mtcb 'ROUG')
        'ET' ('COULEUR' mthe 'BLEU') ;
      'SI' (logst) ;
         mail = 'ET' mail ('COULEUR' mtst 'BLANC') ;
      'FINSI' ;                 
      npomail = 'NBNO' mail ;
      nelmail = 'NBEL' mail ;
      titmail = 'CHAINE' 'Mesh ' 'NBPO=' npomail
                                ' NBELEM=' nelmail ;
      TRACMESH mail titmail VRAI ;
   'FINSI' ;      
*
* Evol precalc
*
   'SI' chzscale ;
      'FORME' orig ;
   'FINSI' ;      
   nchan = 'DIME' (tabin . 'CHANNEL') ;
   tt = 'TABLE' ;
   tev = 'TABLE' ;
   'REPETER' ichan nchan ;
      lighe = tabgeo . 'tlchan' . &ichan ;
      evthec = 'EVOL' 'CHPO' thec lighe ;
      APPEND tev 'TEMPERATURE' evthec ;
      evuhe = 'EVOL' 'CHPO' uz lighe ;
      APPEND tev 'VITESSE' evuhe ;
      tt . &ichan = 'CHAINE' 'Canal ' &ichan ;
   'FIN' ichan ;
   'SI' chzscale ;
      'FORME' defrm ;
   'FINSI' ;      
*
   tcha = 'TABLE' ;
   tcha . 1 = 'T' ;
   tcha . 2 = 'FuelT' ;
   tcha . 3 = 'FuelT2' ;
   tcha . 4 = 'Excoef' ;
   tcha . 5 = 'HeT' ;
   tcha . 6 = 'HeUz' ;
   tcha . 7 = 'HeevT' ;
   tcha . 8 = 'HeevUz' ;
   tcha . 9 = 'Nuss' ;
   tcha . 10 = 'Reyn' ;
   tcha . 11 = 'Prand' ;
   tcha . 12 = 'LamHe' ;
   tcha . 13 = 'Dh' ;
   ntcha = 'DIME' tcha ;
   'SI' interact ;
      itcha = 1 ;
      noclic = VRAI ;
   'SINON' ;
      itcha = 'EXTRAIRE' liview iiview ;
      noclic = FAUX ;
   'FINSI' ;            
*
'REPETER' bouc2 ;
'SI' interact ;
   'MESSAGE' 'Noclic allows not to click between each view' ;   
   'MESSAGE'
'However, one will deactivate it in order to zoom or save' ;
   noclic = 'CHOI' 'Choose desired options' noclic ;
'FINSI' ;   
*
'REPETER' bouc1;
*
'SI' ('EGA' itcha 1) ;
   tit = 'CHAINE' 'Temp. (°C)' ;
   tracchpo  (tcbsct 'ET' thect)
             (mtcbstt 'ET' mthet)
             (amtcbstt 'ET' amthet)
             tit noclic ;
'FINSI' ;
'SI' ('EGA' itcha 2) ;
   tit = 'CHAINE' 'Fuel Temp. (°C)' ;
   tracchpo tcbsc mtcbst amtcbst tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 3) ;
   tit = 'CHAINE' 'Fuel Temp. (°C)' ;
   tracchpo tcbsct mtcbstt amtcbstt tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 4) ;
   tit = 'CHAINE' 'Exchange coeff. (W.m-2.K-1)'  ;
   tracchml $mtif cech ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 5) ;
   tit = 'CHAINE' 'Temp. He (°C)' ;
   tracchpo thect mthe amthe tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 6) ;
   tit = 'CHAINE' 'Uz He (m.s-1)' ;
   tracchpo uz  mthe amthe tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 7) ;
    tit = 'CHAINE' 'Helium Temperature' ;
    titxthe = 'CHAINE' 's (m)' ;
    titythe = 'CHAINE' 'T (°C)' ;
    DESSEVOL (tev . 'TEMPERATURE') tt tit titxthe titythe noclic ;
'FINSI' ;
'SI' ('EGA' itcha 8) ;
    tit = 'CHAINE' 'Helium Vertical Speed' ;
    titxuhe = 'CHAINE' 's (m)' ;
    tityuhe = 'CHAINE' 'uz (m.s-1)' ;
    DESSEVOL (tev . 'VITESSE') tt tit titxuhe tityuhe noclic ;
'FINSI' ;       
'SI' ('EGA' itcha 9) ;
   tit = 'CHAINE' 'Nusselt number'  ;
   tracchml $mtif nuss0 ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 10) ;
   tit = 'CHAINE' 'Reynolds number'  ;
   tracchml $mtif reyn0 ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 11) ;
   tit = 'CHAINE' 'Prandtl number'  ;
   tracchml $mtif pran0 ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 12) ;
   tit = 'CHAINE' 'Lambda (He)'  ;
   tracchml $mtif l0 ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 13) ;
   tit = 'CHAINE' 'Hydraulic diameter (m)'  ;
   tracchml $mtif Dh ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
*
'SI' interact ;
   cha = tit ;
   ret = 'MENU' cha 'Next' 'Previous' 'Options' ;
   'SI' ('EGA' ret 'Options') ; 'QUITTER' bouc1  ; 'FINSI';
   'SI' ('EGA' ret 'Quitter') ; 'QUITTER' bouc2; 'FINSI';
   'SI' ('EGA' ret 'Next') ; itcha = '+' itcha 1 ; 'FINSI' ; 
   'SI' ('EGA' ret 'Previous') ; itcha = '-' itcha 1 ; 'FINSI' ; 
   'SI' (itcha > ntcha); itcha = 1      ; 'FINSI';
   'SI' (itcha < 1); itcha = ntcha      ; 'FINSI';
'SINON' ;
   iiview = '+' iiview 1 ;
   'SI' ('>' iiview niview) ;
      'QUITTER' bouc2 ;
   'SINON' ;
      itcha = 'EXTRAIRE' liview iiview ;
   'FINSI' ;
'FINSI' ;            
*
'FIN' bouc1;
'FIN' bouc2 ;
*
'SI' chzscale ;
   'FORME' orig ;
'FINSI' ;      
*
* End of procedure file EXPLORES
*
'FINPROC' ;
*ENDPROCEDUR explores
*BEGINPROCEDUR genmema
************************************************************************
* NOM         : GENMEMA
* DESCRIPTION : Generate meshes and coupling matrix given base surface
*               mesh
*
* Input for GENMEMA procedure
* tabin . 'INTERFACE'  . i : i-th interface       (line mesh)
* tabin . 'CHANNEL'    . i : i-th helium channel  (surface mesh)
* tabin . 'FUEL'           : fuel mesh (surface mesh)
* tabin . 'STRUCTURE'      : structure mesh (surface mesh) (optional)
* H : height (m) and nh : number of elements in the height
* (regular meshing in z direction)
*
* Temporary table for GENMEMA procedure
* tabout . 'INTERFACE'  . i : i-th interface       (surface mesh)
* tabout . 'CHANNEL'    . i : i-th helium channel  (volume mesh)
* tabout . 'CCONINCH'   . i : i-th coupling connectivity intreface <-> channel
* tabout . 'CMATCHIN'   . i : i-th coupling matrix channel -> interface
* tabout . 'CMATINCH'   . i : i-th coupling matrix interface -> channel
* tabin . 'FUEL'           : fuel mesh (volume mesh)
* tabin . 'STRUCTURE'      : structure mesh (volume mesh) (optional)
*
* Output table for GENMEMA procedure
*
* output tabgeo should contain all the necessary geometric information
* needed further
* tabgeo . 'tabin'   = tabin
* tabgeo . 'H'       = H   : height
* tabgeo . 'nh'      = nh  : number of elements in the height
* tabgeo . '$mtcbst' = NS model objet (fuel+structure)
* tabgeo . '$mtcb'   = NS model objet (fuel only)
* tabgeo . '$mthe'   = NS model objet (helium channels)
* tabgeo . '$sthe'   = NS model objet (bottom of the helium channels)
* tabgeo . '$mtif'   = NS model objet (interfaces)
* tabgeo . '$mtifb'  = NS model objet (blocked interfaces)
* tabgeo . 'mcf'     = coupling matrix helium -> fuel
* tabgeo . 'mfc'     = coupling matrix fuel -> helium
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 09/11/2004, version initiale
* HISTORIQUE : v1, 09/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GENMEMA ;
'ARGUMENT' tabin*'TABLE' ;
'ARGUMENT' H*'FLOTTANT' ;
'ARGUMENT' nh*'ENTIER' ;
*
tabout = 'TABLE' ;
tabout . 'INTERFACE' = 'TABLE' ;
tabout . 'CHANNEL'   = 'TABLE' ;
tabout . 'CCONINCH'  = 'TABLE' ;
tabout . 'CMATCHIN'  = 'TABLE' ;
tabout . 'CMATINCH'  = 'TABLE' ;
*
hnh    = '/' H ('FLOTTANT' nh) ;
vechnh = (0. 0. hnh) ;
*
* Treatment of everything except 'OTHER'
*
ninterf = 'DIME' (tabin . 'INTERFACE') ;
nchanne = 'DIME' (tabin . 'CHANNEL') ;
'SI' ('NEG' ninterf nchanne) ;
   cherr = 'CHAINE' 'NEG (DIME (tabin . INTERFACE)) '
                    '(DIME (tabin . CHANNEL))' ;
   'ERREUR' cherr ;                    
'FINSI' ;
*
* Loop on interfaces
*
*   Préconditionnement
*
   mott = 'MOTS' 'T' ;
   pr1  = 'PROG' 0.D0 0.D0 ;
   pr3  = 'PROG' 0.D0 1.D0 ;
*
'REPETER' iinterf ninterf ;
   lpreums = vrai ;
   linterf = tabin  . 'INTERFACE' . &iinterf ;
   scan    = tabin  . 'CHANNEL'   . &iinterf ;
* Loop on height   
   'REPETER' ih nh ;
      vectran = (0. 0. ('*' hnh ('FLOTTANT' ('-' &ih 1)))) ;
      tli = 'PLUS' linterf vectran ;
      tsc = 'PLUS' scan    vectran ;
      sli = 'TRANSLATION' tli 1 vechnh ;
      vsc  = 'VOLUME' tsc 'TRAN' 1 vechnh ;
      cvsc = 'BARYCENTRE' vsc ;
      vvsc = 'MESURE' vsc ;
      ivvsc = '/' 1.D0 vvsc ;
      pr2i  = 'PROG' ivvsc 0.D0 ;
*      'LISTE' cvsc ; 'LISTE' cvsc2 ;
      APPEND (tabout . 'CHANNEL')  &iinterf vsc ;
* Loop on interface elements
      nsli = 'NBEL' sli ;
      'REPETER' isli nsli ;
         esli  = 'ELEM' sli &isli ;
         cesli = 'BARYCENTRE' esli ;
         cco = 'MANUEL' 'SEG2' cesli cvsc ;
         mci = 'MANUEL' 'RIGIDITE' cco mott 'QUEL' pr3 pr1 ;
*                ('PROG' 0.D0 1.D0) ('PROG' 0.D0 0.D0) ;
         mic = 'MANUEL' 'RIGIDITE' cco mott 'QUEL' pr1 pr2i ;
         'SI' lpreums ;
            tabout . 'INTERFACE' . &iinterf = esli ;
            tabout . 'CCONINCH'  . &iinterf = cco ;
            tabout . 'CMATCHIN'  . &iinterf = mci ;
            tabout . 'CMATINCH'  . &iinterf = mic ;
            lpreums = faux ;
         'SINON' ;
            tabout . 'INTERFACE' . &iinterf = 'ET'
               (tabout . 'INTERFACE' . &iinterf) esli ;
            tabout . 'CCONINCH'  . &iinterf = 'ET'
               (tabout . 'CCONINCH'  . &iinterf) cco ;
            tabout . 'CMATCHIN'  . &iinterf = 'ET' 
               (tabout . 'CMATCHIN'  . &iinterf) mci ;
            tabout . 'CMATINCH'  . &iinterf = 'ET'
               (tabout . 'CMATINCH' . &iinterf) mic ;
         'FINSI' ;
      'FIN' isli ;
   'FIN' ih ;
'FIN' iinterf ;                            
*
* Treatment of 'FUEL' and 'STRUCTURE'
*
sother = tabin . 'FUEL' ;
vother = 'VOLUME' sother 'TRAN' nh (0. 0. H) ;
tabout . 'FUEL' = vother ;
tabout . 'FUEL2' = 'VOLUME' sother 'TRAN' 1 (0. 0. H) ;
'SI' ('EXISTE' tabin 'STRUCTURE') ;
   sother = tabin . 'STRUCTURE' ;
   vother = 'VOLUME' sother 'TRAN' nh (0. 0. H) ;
   tabout . 'STRUCTURE' = vother ;
   tabout . 'STRUCT2'   = 'VOLUME' sother 'TRAN' 1 (0. 0. H) ;
'FINSI' ;   
*
* Some eliminations
*
tinter = @STBL (tabout . 'INTERFACE') ;
tother = tabout . 'FUEL' ;
'SI' ('EXISTE' tabout 'STRUCTURE') ;
   tother = 'ET' tother (tabout . 'STRUCTURE') ;
'FINSI' ;   
'ELIMINATION' ('ET' tinter tother) 1.D-6 ;
*
* Second part : tabgeo
*
tabgeo = 'TABLE' ;
*
* Build the global meshes
*
*  Structure & Combustible
mtcb   = tabout . 'FUEL' ;
amtcb  = 'ARETE' (tabout . 'FUEL2') ;
'SI' ('EXISTE' tabout 'STRUCTURE') ;
   mtcbst  = 'ET' mtcb (tabout . 'STRUCTURE') ;
   amtcbst = 'ET' amtcb ('ARETE' (tabout . 'STRUCT2')) ;
'SINON' ;
   mtcbst  = mtcb ;
   amtcbst = amtcb ;
'FINSI' ;   
*  Channels
*mthe = @STBL (tabout . 'CHANNEL') ;
tmthe = tabout . 'CHANNEL' ;
sthe =  @STBL (tabin  . 'CHANNEL') ;
amthe = 'TRANSLATION' ('CONTOUR' sthe) 1 (0. 0. H) ;
*  Line in channels
nchan = 'DIME' (tabin . 'CHANNEL') ;
tlchan = 'TABLE' ;
tpechan = 'TABLE' ;
tpschan = 'TABLE' ;
'REPETER' ichan nchan ;
   selem1 = 'ELEM' (tabin . 'CHANNEL' . &ichan) 1 ;
   po1    = 'POIN' ('CONTOUR' selem1) 1 ;
   po2    =  'PLUS' po1 (0. 0. H) ;
   lighe  = 'DROIT' nh po1 po2 ;
   tlchan . &ichan  = lighe ;
   tpechan . &ichan = po1 ;
   tpschan . &ichan = po2 ;
'FIN' ichan ;   
*  Interfaces
tmtif = tabout . 'INTERFACE' ;
mtif = @STBL (tabout . 'INTERFACE') ;
* Blocked interfaces
lblock = tabin . 'BLOCKED' ;
dblock = 'DIME' lblock ;
logb =  ('>' dblock 0) ;
'SI' logb ;
   tchan = 'TABLE' ;
   'REPETER' iblock dblock ;
      iib = 'EXTRAIRE' lblock &iblock ;
      APPEND tchan 'mtifb' (tabout . 'INTERFACE' . iib) ;
   'FIN' iblock ;
   mtifb = tchan . 'mtifb' ;      
'FINSI' ;
*  Coupling
cfc = @STBL (tabout . 'CCONINCH') ;
mcf = @STBL (tabout . 'CMATCHIN') ;
mfc = @STBL (tabout . 'CMATINCH') ;
mcf = 'CHANGER' mcf 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
mfc = 'CHANGER' mfc 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
*  Generate the QUAF meshes
*   Channels
nhe = 'DIME' tmthe ;
_tmthe = 'TABLE' ;
'REPETER' ihe nhe ;
   _tmthe . &ihe = QUAFME (tmthe . &ihe) ;
'FIN' ihe ;
_mthe = @STBL _tmthe ;
*   Interface
nif = 'DIME' tmtif ;
_tmtif = 'TABLE' ;
'REPETER' iif nif ;
   _tmtif . &iif = QUAFME (tmtif . &iif) ;
'FIN' iif ;
_mtif = @STBL _tmtif ;
*   Other
 _mtcbst _mtcb _sthe =
   QUAFME mtcbst mtcb sthe ;
mcbt = _mtcbst 'ET' _mtcb 'ET' _mtif 'ET' cfc
       'ET' amtcb 'ET' amtcbst ;
'SI' logb ;
   _mtifb = QUAFME mtifb ;
   mcbt = 'ET' mcbt _mtifb ;
'FINSI' ;   
*  Elimination of points with same coordinates
'ELIMINATION' mcbt 1.D-6 ;
ttlchan = @STBL tlchan ;
'ELIMINATION' (_mthe 'ET' _sthe 'ET' cfc 'ET' amthe 'ET' ttlchan)
              1.D-6 ;
* Create the Navier-Stokes model objects
*  Channels
$tmthe = 'TABLE' ;
'REPETER' ihe nhe ;
   $tmthe . &ihe = MODNAV 'LINE' (_tmthe . &ihe) ;
'FIN' ihe ;
*  Interface
$tmtif = 'TABLE' ;
'REPETER' iif nif ;
   $tmtif . &iif = MODNAV 'LINE' (_tmtif . &iif) ;
'FIN' iif ;
$mtcbst $mtcb $mthe $mtif $sthe =
   MODNAV 'LINE' _mtcbst _mtcb _mthe _mtif _sthe ;
* Create the hydraulic diameter table 
tdhif = 'TABLE' ;
'REPETER' iif nif ;
   tdhif . &iif = 'KCHT' ($tmtif . &iif) 'SCAL' 'CENTRE'
                          (tabin . 'DH' .&iif) ;
'FIN' iif ;
dh = @STBL tdhif ;
*
* Store results
*
tabgeo . 'tabin'   = tabin   ;
tabgeo . 'H'       = H       ;
tabgeo . 'nh'      = nh      ;
tabgeo . '$mtcbst' = $mtcbst ;
tabgeo . '$mtcb'   = $mtcb   ;
tabgeo . '$tmthe'  = $tmthe   ;
tabgeo . '$mthe'   = $mthe   ;
tabgeo . '$sthe'   = $sthe   ;
tabgeo . '$mtif'   = $mtif   ;
tabgeo . 'amtcb'   = amtcb   ;
tabgeo . 'amtcbst' = amtcbst ;
tabgeo . 'amthe'   = amthe   ;
tabgeo . 'mcf'     = mcf     ;
tabgeo . 'mfc'     = mfc     ;
tabgeo . 'cfc'     = cfc     ;
'SI' logb ;
   $mtifb = MODNAV 'LINE' _mtifb ;
   tabgeo . '$mtifb' = $mtifb ;
'FINSI' ;
tabgeo . 'Dh' = dh ;
tabgeo . 'tlchan' = tlchan ;
tabgeo . 'tpechan' = tpechan ;
tabgeo . 'tpschan' = tpschan ;
'RESPRO' tabgeo ;
*
* End of procedure file GENMEMA
*
'FINPROC' ;
*ENDPROCEDUR genmema
*BEGINPROCEDUR geomhex1
************************************************************************
* NOM         : GEOMHEX1
* DESCRIPTION : 30° symmetry cell of an hexagonal geometry
*               in the x-y plane
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, ??/??/2004, version initiale
* HISTORIQUE : v1, ??/??/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMHEX1 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' pas*'FLOTTANT' ;
*
* Geometric quantities deduced from the one given above
*
q = '/' pas 2.D0 ;
R = '/' Dh 2.D0 ;
*
*  Geometric discretization parameters 
*  (number of elements)
*  nb (base) : number of elements for the line in (x-y) plane
*              with smallest length
*
nb = 7 ;
s30 = 'SIN' 30.D0 ;
c30 = 'COS' 30.D0 ;
t30 = '/' s30 c30 ;
* length of each line
lcer = ('*' ('/' PI 6.D0) R) ;
lpdia = '-' q R ;
lgdia = '-' ('/' q c30) R ;
lhex  = ('*' q t30) ;
* scaling based on shortest line
lbase = 'MINIMUM' ('PROG' lcer lpdia lgdia lhex) ;
* factors in order to have an almost regular mesh
fcer  = '/' lcer lbase ;
fpdia = '/' lpdia lbase ;
fgdia = '/' lgdia lbase ;
fhex  = '/' lhex lbase ;
* mesh
rbase = 'FLOTTANT' nb ;
ncer  = '*' rbase fcer ;
nhex  = '*' rbase fhex ;
npdia  = '*' rbase fpdia ;
ngdia  = '*' rbase fgdia ;
*
ncer = 'ENTIER' ('+' ncer 0.99D0) ;
nhex = 'ENTIER' ('+' nhex 0.99D0) ;
npdia = 'ENTIER' ('+' npdia 0.99D0) ;
ngdia = 'ENTIER' ('+' ngdia 0.99D0) ;
* Points definition for fuel
p0 = 0.D0 0.D0 0.D0 ;
p1 = R    0.D0 0.D0 ;
p2 = q    0.D0 0.D0 ;
p3 = q    ('*' q t30) 0.D0 ;
p4 = ('*' R c30) ('*' R s30) 0.D0 ;
* Lines definition for fuel
l1 = 'DROIT' npdia p1 p2 ;
l2 = 'DROIT' nhex  p2 p3 ;
l3 = 'DROIT' ngdia p3 p4 ;
l4 = 'CERCLE' ncer p4 p0 p1 ;
* Surfaces definition for fuel
ctn  = l1 'ET' l2 'ET' l3 'ET' l4 ;
stcb = 'SURFACE' ctn 'PLANE' ;
* Points and surface definitions for helium (an isocele triangle with
* surface equal to the corresponding circular sector
Rp = '*' R ('**' ('/' PI 3.D0) 0.5D0) ;
p5 = Rp   0.D0 0.D0 ;
p6 = ('*' Rp c30) ('*' Rp s30) 0.D0 ;
sthe = 'MANUEL' 'TRI3' p0 p5 p6 ;
*
* Output results in a table suitable for GENMEMA
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
*
tabin . 'INTERFACE' . 1 = l4 ;
tabin . 'CHANNEL'   . 1 = sthe ;
tabin . 'FUEL'          = stcb ;
'RESPRO' tabin ;
*
* End of procedure file GEOMHEX1
*
'FINPROC' ;
*ENDPROCEDUR geomhex1
*BEGINPROCEDUR geomhex2
************************************************************************
* NOM         : GEOMHEX2
* DESCRIPTION : Same as GEOMHEX1 but larger symmetry cell
*               in the x-y plane
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 18/11/2004, version initiale
* HISTORIQUE : v1, 18/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMHEX2 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' pas*'FLOTTANT' ;
*
* Geometric quantities deduced from the one given above
*
q = '/' pas 2.D0 ;
R = '/' Dh 2.D0 ;
*
*  Geometric discretization parameters 
*  (number of elements)
*  nb (base) : number of elements for the line in (x-y) plane
*              with smallest length
*
nb = 7 ;
s30 = 'SIN' 30.D0 ;
c30 = 'COS' 30.D0 ;
t30 = '/' s30 c30 ;
* length of each line
lcera = ('*' ('/' PI 6.D0) R) ;
lcerb = ('*' ('/' PI 2.D0) R) ;
lpdia = '-' pas Dh ;
lgdia = '-' ('/' pas c30) R ;
lhex  = '-' ('*' pas t30) R ;
* scaling based on shortest line
lbase = 'MINIMUM' ('PROG' lcera lcerb lpdia lgdia lhex) ;
* factors in order to have an almost regular mesh
fcera = '/' lcera lbase ;
fcerb = '/' lcerb lbase ;
fpdia = '/' lpdia lbase ;
fgdia = '/' lgdia lbase ;
fhex  = '/' lhex lbase ;
* mesh
rbase = 'FLOTTANT' nb ;
ncera  = '*' rbase fcera ;
ncerb  = '*' rbase fcerb ;
nhex  = '*' rbase fhex ;
npdia  = '*' rbase fpdia ;
ngdia  = '*' rbase fgdia ;
*
ncera = 'ENTIER' ('+' ncera 0.99D0) ;
ncerb = 'ENTIER' ('+' ncerb 0.99D0) ;
nhex = 'ENTIER' ('+' nhex 0.99D0) ;
npdia = 'ENTIER' ('+' npdia 0.99D0) ;
ngdia = 'ENTIER' ('+' ngdia 0.99D0) ;
* Points definition for fuel
p0 = 0.D0 0.D0 0.D0 ;
p1 = R    0.D0 0.D0 ;
p2 = ('-' pas R) 0.D0 0.D0 ;
p3 = pas  0.D0 0.D0 ;
p4 = pas  R 0.D0    ;
p5 = pas  ('*' pas t30) 0.D0 ;
p6 = ('*' R c30) ('*' R s30) 0.D0 ;
* Lines definition for fuel
l1 = 'DROIT' npdia p1 p2 ;
l2 = 'CERCLE' ncerb p2 p3 p4 ;
l3 = 'DROIT' nhex  p4 p5 ;
l4 = 'DROIT' ngdia p5 p6 ;
l5 = 'CERCLE' ncera p6 p0 p1 ;
* Surfaces definition for fuel
ctn  = l1 'ET' l2 'ET' l3 'ET' l4 'ET' l5 ;
stcb = 'SURFACE' ctn 'PLANE' ;
* Points and surface definitions for helium (two isocele triangles with
* surface equal to the corresponding circular sector
Rpa = '*' R ('**' ('/' PI 3.D0) 0.5D0) ;
p7 = Rpa   0.D0  0.D0 ;
p8 = ('*' Rpa c30) ('*' Rpa s30) 0.D0 ;
sthea = 'MANUEL' 'TRI3' p0 p7 p8 ;
Rpb = '*' R ('**' ('/' PI 2.D0) 0.5D0) ;
p9  = ('-' pas Rpb) 0.D0 0.D0 ;
p10 = pas Rpb 0.D0 ;
stheb = 'MANUEL' 'TRI3' p9 p3 p10 ;
*
* Output results in a table suitable for GENMEMA
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
*
tabin . 'INTERFACE' . 1 = l5 ;
tabin . 'CHANNEL'   . 1 = sthea ;
tabin . 'INTERFACE' . 2 = l2 ;
tabin . 'CHANNEL'   . 2 = stheb ;
tabin . 'FUEL'          = stcb ;
'RESPRO' tabin ;
*
* End of procedure file GEOMHEX2
*
'FINPROC' ;
*ENDPROCEDUR geomhex2
*BEGINPROCEDUR geompl21
************************************************************************
* NOM         : GEOMPL21
* DESCRIPTION : Type 1 symmetry cell for GEOMPLA2
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/11/2004, version initiale
* HISTORIQUE : v1, 23/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPL21 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
'ARGUMENT' nc2*'ENTIER' ;
'ARGUMENT' nf*'ENTIER' ;
'ARGUMENT' np*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
L4 = '/' L2 2.D0 ;
*
* Discretisation
* 
nplate2 = '/' nplate 2 ;
nc = '*' nc2 2 ;
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* Part1 : Lines
*
p0 = 0. 0. 0. ; 
p1 = 0. ('*' ef -1.D0) 0. ;
lcor = 'DROIT' nf p0 p1 ;
tablig = 'TABLE' ;
bp = p1 ;
'REPETER' ipl nplate2 ;
   p2 = 'PLUS' bp (0. ('*' ep -1.D0) 0.) ;
   APPEND tablig 'LIG' ('DROIT' bp p2 np) ;
   bp = p2 ;
* The last one is an half channel
   'SI' ('NEG' &ipl nplate2) ;   
      p2 = 'PLUS' bp (0. ('*' ec -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc) ;
      bp = p2 ;
   'SINON' ;      
      p2 = 'PLUS' bp (0. ('*' ec2 -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc2) ;
      bp = p2 ;
   'FINSI' ;
'FIN' ipl ;
lcb = tablig . 'LIG' ;
p2  = 'PLUS' p0 bp ;
lch = 'DROIT' p1 p2 1 ;
*
* Part1 : left structure
*
bp = 0. 0. 0. ;
vt = ef 0. 0. ;
APPEND tabin 'STRUCTURE' ('TRANSLATION' (lcor 'ET' lcb) nf vt) ;
bp = 'PLUS' bp vt ;
*
* Loop on plates
*
'REPETER' iplate nplate ;
*  Structure part
   vt = ep 0. 0. ;
   lcort = 'PLUS' lcor bp ;
   APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort np vt) ;
* Fuel part
   lcbt = 'PLUS' lcb bp ;
   APPEND tabin 'FUEL'  ('TRANSLATION' lcbt np vt) ;
   lcbtt = 'PLUS' lcbt vt ;
   'SI' ('>' &iplate 1) ;
      APPEND (tabin . 'INTERFACE') ('-' iitf 1) lcbt ;
   'FINSI' ;      
   'SI' ('<' &iplate nplate) ;
      APPEND (tabin . 'INTERFACE') iitf lcbtt ;
   'FINSI' ;      
   bp = 'PLUS' bp vt ;
   'SI' ('<' &iplate nplate) ;
* Structure part 
      vt = ec 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p2t nc) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
      bp = 'PLUS' bp vt ;
   'FINSI' ;      
'FIN' iplate ;    
* Part1 : right structure
vt = ef 0. 0. ;
lcct = 'PLUS' ('ET' lcor lcb) bp ;
APPEND tabin 'STRUCTURE' ('TRANSLATION' lcct nf vt) ;
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPL21
*
'FINPROC' ;
*ENDPROCEDUR geompl21
*BEGINPROCEDUR geompl22
************************************************************************
* NOM         : GEOMPL22
* DESCRIPTION : Type 2 symmetry cell for GEOMPLA2
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/11/2004, version initiale
* HISTORIQUE : v1, 23/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPL22 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
'ARGUMENT' nc2*'ENTIER' ;
'ARGUMENT' nf*'ENTIER' ;
'ARGUMENT' np*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
L4 = '/' L2 2.D0 ;
*
* Discretisation
* 
nplate2 = '/' nplate 2 ;
nc = '*' nc2 2 ;
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* Part1 : Lines
*
p0 = 0. 0. 0. ; 
p1 = 0. ('*' ef -1.D0) 0. ;
lcor = 'DROIT' nf p0 p1 ;
tablig = 'TABLE' ;
bp = p1 ;
'REPETER' ipl nplate2 ;
   p2 = 'PLUS' bp (0. ('*' ep -1.D0) 0.) ;
   APPEND tablig 'LIG' ('DROIT' bp p2 np) ;
   bp = p2 ;
* The last one is an half channel
   'SI' ('NEG' &ipl nplate2) ;   
      p2 = 'PLUS' bp (0. ('*' ec -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc) ;
      bp = p2 ;
   'SINON' ;      
      p2 = 'PLUS' bp (0. ('*' ec2 -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc2) ;
      bp = p2 ;
   'FINSI' ;
'FIN' ipl ;
lcb = tablig . 'LIG' ;
p2  = 'PLUS' p0 bp ;
lch = 'DROIT' p1 p2 1 ;
*
* Part1 : left structure
*
bp = 0. 0. 0. ;
vt = ef 0. 0. ;
APPEND tabin 'STRUCTURE' ('TRANSLATION' (lcor 'ET' lcb) nf vt) ;
bp = 'PLUS' bp vt ;
*
* Loop on plates
*
'REPETER' iplate nplate2 ;
*  Structure part
   vt = ep 0. 0. ;
   lcort = 'PLUS' lcor bp ;
   APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort np vt) ;
* Fuel part
   lcbt = 'PLUS' lcb bp ;
   APPEND tabin 'FUEL'  ('TRANSLATION' lcbt np vt) ;
   lcbtt = 'PLUS' lcbt vt ;
   'SI' ('>' &iplate 1) ;
      APPEND (tabin . 'INTERFACE') ('-' iitf 1) lcbt ;
   'FINSI' ;      
   APPEND (tabin . 'INTERFACE') iitf lcbtt ;
   bp = 'PLUS' bp vt ;
   'SI' ('<' &iplate nplate2) ;
* Structure part 
      vt = ec 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p2t nc) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
      bp = 'PLUS' bp vt ;
* The last one is a half channel      
   'SINON' ;      
* Structure part 
      vt = ec2 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc2 vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p2t nc2) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
   'FINSI' ;      
'FIN' iplate ;
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPL22
*
'FINPROC' ;
*ENDPROCEDUR geompl22
*BEGINPROCEDUR geompl23
************************************************************************
* NOM         : GEOMPL23
* DESCRIPTION : Type 3 symmetry cell for GEOMPLA2
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/11/2004, version initiale
* HISTORIQUE : v1, 23/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPL23 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
'ARGUMENT' nc2*'ENTIER' ;
'ARGUMENT' nf*'ENTIER' ;
'ARGUMENT' np*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
L4 = '/' L2 2.D0 ;
*
* Discretisation
* 
nplate2 = '/' nplate 2 ;
nc = '*' nc2 2 ;
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* Part1 : Lines
*
p0 = 0. 0. 0. ; 
p1 = ef 0. 0. ;
lcorg = 'DROIT' nf p0 p1 ;
tablig = 'TABLE' ;
bp = p1 ;
'REPETER' iplate nplate ;
   p2 = 'PLUS' bp (ep 0. 0.) ;
   APPEND tablig 'LIG' ('DROIT' bp p2 np) ;
   bp = p2 ;
   'SI' ('<' &iplate nplate) ;   
      p2 = 'PLUS' bp (ec 0. 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc) ;
      bp = p2 ;
   'FINSI' ;
'FIN' iplate ;
lcb = tablig . 'LIG' ;
p2  = 'PLUS' p0 bp ;
lch = 'DROIT' p1 p2 1 ;
bp = p2 ;
p3 = 'PLUS' bp (ef 0. 0.) ;
lcord = 'DROIT' nf p2 p3 ;
lcor = 'ET' lcorg lcord ;
*
* Part1 : bottom structure
*
bp = 0. 0. 0. ;
vt = 0. ef 0. ;
APPEND tabin 'STRUCTURE' ('TRANSLATION' (lcor 'ET' lcb) nf vt) ;
bp = 'PLUS' bp vt ;
*
* Loop on plates
*
'REPETER' iplate nplate2 ;
*  Structure part
   vt = 0. ep 0. ;
   lcort = 'PLUS' lcor bp ;
   APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort np vt) ;
* Fuel part
   lcbt = 'PLUS' lcb bp ;
   APPEND tabin 'FUEL'  ('TRANSLATION' lcbt np vt) ;
   lcbtt = 'PLUS' lcbt vt ;
   'SI' ('>' &iplate 1) ;
      APPEND (tabin . 'INTERFACE') ('-' iitf 1) lcbt ;
   'FINSI' ;      
   APPEND (tabin . 'INTERFACE') iitf lcbtt ;
   bp = 'PLUS' bp vt ;
   'SI' ('<' &iplate nplate2) ;
* Structure part 
      vt = 0. ec 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t    = 'PLUS' p1   bp ;
      p1tt   = 'PLUS' p1t  vt ;
      p2t    = 'PLUS' p2   bp ;
      p2tt   = 'PLUS' p2t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p1tt nc) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p2t p2tt nc) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
      bp = 'PLUS' bp vt ;
* The last one is a half channel      
   'SINON' ;      
* Structure part 
      vt = 0. ec2 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t    = 'PLUS' p1   bp ;
      p1tt   = 'PLUS' p1t  vt ;
      p2t    = 'PLUS' p2   bp ;
      p2tt   = 'PLUS' p2t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc2 vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p1tt nc2) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p2t p2tt nc2) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
   'FINSI' ;      
'FIN' iplate ;
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPL23
*
'FINPROC' ;
*ENDPROCEDUR geompl23
*BEGINPROCEDUR geompla1
************************************************************************
* NOM         : GEOMPLA1
* DESCRIPTION : Mesh of a plate geometry : small symmetry cell
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 22/11/2004, version initiale
* HISTORIQUE : v1, 22/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPLA1 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
L4 = '/' L2 2.D0 ;
*
* Discretisation
* 
nplate2 = '/' nplate 2 ;
nc2 = 1 ;
nc = '*' nc2 2 ;
nf = 2 ;
np = 2 ;
*
* Tables
*
tabin1 = 'TABLE' ;
tabin1 . 'INTERFACE' = 'TABLE' ;
tabin1 . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* Part1 : Lines
*
p0 = 0. 0. 0. ; 
p1 = 0. ('*' ef -1.D0) 0. ;
lcor = 'DROIT' nf p0 p1 ;
tablig = 'TABLE' ;
bp = p1 ;
'REPETER' ipl nplate2 ;
   p2 = 'PLUS' bp (0. ('*' ep -1.D0) 0.) ;
   APPEND tablig 'LIG' ('DROIT' bp p2 np) ;
   bp = p2 ;
* The last one is an half channel
   'SI' ('NEG' &ipl nplate2) ;   
      p2 = 'PLUS' bp (0. ('*' ec -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc) ;
      bp = p2 ;
   'SINON' ;      
      p2 = 'PLUS' bp (0. ('*' ec2 -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc2) ;
      bp = p2 ;
   'FINSI' ;
'FIN' ipl ;
lcb = tablig . 'LIG' ;
p2  = 'PLUS' p0 bp ;
lch = 'DROIT' p1 p2 1 ;
*
* Part1 : left structure
*
bp = 0. 0. 0. ;
vt = ef 0. 0. ;
APPEND tabin1 'STRUCTURE' ('TRANSLATION' (lcor 'ET' lcb) nf vt) ;
bp = 'PLUS' bp vt ;
*
* Loop on plates
*
'REPETER' iplate nplate2 ;
*  Structure part
   vt = ep 0. 0. ;
   lcort = 'PLUS' lcor bp ;
   APPEND tabin1 'STRUCTURE' ('TRANSLATION' lcort np vt) ;
* Fuel part
   lcbt = 'PLUS' lcb bp ;
   APPEND tabin1 'FUEL'  ('TRANSLATION' lcbt np vt) ;
   lcbtt = 'PLUS' lcbt vt ;
   'SI' ('>' &iplate 1) ;
      APPEND (tabin1 . 'INTERFACE') ('-' iitf 1) lcbt ;
   'FINSI' ;      
   APPEND (tabin1 . 'INTERFACE') iitf lcbtt ;
   bp = 'PLUS' bp vt ;
   'SI' ('<' &iplate nplate2) ;
* Structure part 
      vt = ec 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin1 'STRUCTURE' ('TRANSLATION' lcort nc vt) ;
      APPEND (tabin1 . 'INTERFACE') iitf ('DROIT' p1t p2t nc) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin1 . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
      bp = 'PLUS' bp vt ;
* The last one is a half channel      
   'SINON' ;      
* Structure part 
      vt = ec2 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin1 'STRUCTURE' ('TRANSLATION' lcort nc2 vt) ;
      APPEND (tabin1 . 'INTERFACE') iitf ('DROIT' p1t p2t nc2) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin1 . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
   'FINSI' ;      
'FIN' iplate ;
nitf = ('-' iitf 1) ;
*
* Perform rotations
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
pr = 0. 0. 1. ;
'REPETER' bclrot 4 ;
   arot = '*' 90. ('-' &bclrot 1) ;
   append tabin 'STRUCTURE'
      ('TOURNER' (tabin1 . 'STRUCTURE') arot p0 pr) ;
   append tabin 'FUEL'
      ('TOURNER' (tabin1 . 'FUEL') arot p0 pr) ;
   'REPETER' bitf nitf ;
      append (tabin . 'INTERFACE') iitf
         ('TOURNER' (tabin1 . 'INTERFACE' . &bitf) arot p0 pr) ;
      append (tabin . 'CHANNEL') iitf
         ('TOURNER' (tabin1 . 'CHANNEL' . &bitf) arot p0 pr) ;
      iitf = '+' iitf 1 ;
   'FIN' bitf ;
'FIN' bclrot ;            
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPLA1
*
'FINPROC' ;
*ENDPROCEDUR geompla1
*BEGINPROCEDUR geompla2
************************************************************************
* NOM         : GEOMPLA2
* DESCRIPTION : Mesh of a plate geometry
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 22/11/2004, version initiale
* HISTORIQUE : v1, 22/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPLA2 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
mL2 = '*' L2 -1.D0 ;
L4 = '/' L2 2.D0 ;
*
* Discretization
*
nc2 = 1 ;
nf  = 2 ;
np  = 2 ;
*
* Type 1 cells 
*
tabin1 = GEOMPL21 Dh ef ep nchan nc2 nf np ;
*
* Type 2 cells
*
tabin2 = GEOMPL22 Dh ef ep nchan nc2 nf np ;
*
* Type 3 cells
*
tabin3 = GEOMPL23 Dh ef ep nchan nc2 nf np ;
*
* Perform rotations and translations to generate total mesh
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* part1 : Type 1 cells
*
nitf = 'DIME' (tabin1 . 'INTERFACE') ;
append tabin 'STRUCTURE' (tabin1 . 'STRUCTURE') ;
append tabin 'FUEL'      (tabin1 . 'FUEL') ;
'REPETER' bitf nitf ;
   append (tabin . 'INTERFACE') iitf
          (tabin1 . 'INTERFACE' . &bitf) ;
   append (tabin . 'CHANNEL') iitf
          (tabin1 . 'CHANNEL' . &bitf) ;
   iitf = '+' iitf 1 ;
'FIN' bitf ;
*
* part2 : two Type2 cells
*
'REPETER' btype2 2 ;
   'SI' ('EGA' &btype2 1) ;
      arot = 90.D0 ;
      vtra = (L2 0.D0 0.D0) ;
      lsym = VRAI ;
   'FINSI' ;    
   'SI' ('EGA' &btype2 2) ;
      arot = 180.D0 ;
      vtra = (L2 0.D0 0.D0) ;
      lsym = VRAI ;
   'FINSI' ;    
   nitf = 'DIME' (tabin2 . 'INTERFACE') ;
   append tabin 'STRUCTURE'
      (SYTOPL (tabin2 . 'STRUCTURE') arot vtra lsym) ;
   append tabin 'FUEL'
      (SYTOPL (tabin2 . 'FUEL') arot vtra lsym) ;
   'REPETER' bitf nitf ;
      append (tabin . 'INTERFACE') iitf
         (SYTOPL (tabin2 . 'INTERFACE' . &bitf) arot vtra lsym) ;
      append (tabin . 'CHANNEL') iitf
         (SYTOPL (tabin2 . 'CHANNEL' . &bitf) arot vtra lsym) ;
      iitf = '+' iitf 1 ;
   'FIN' bitf ;
'FIN' btype2 ;   
*
* part3 : Type 1 cells
*
nitf = 'DIME' (tabin3 . 'INTERFACE') ;
append tabin 'STRUCTURE' (tabin3 . 'STRUCTURE') ;
append tabin 'FUEL'      (tabin3 . 'FUEL') ;
'REPETER' bitf nitf ;
   append (tabin . 'INTERFACE') iitf
          (tabin3 . 'INTERFACE' . &bitf) ;
   append (tabin . 'CHANNEL') iitf
          (tabin3 . 'CHANNEL' . &bitf) ;
   iitf = '+' iitf 1 ;
'FIN' bitf ;
*
* part4 : two Type2 cells
*
'REPETER' btype2 2 ;
   'SI' ('EGA' &btype2 1) ;
      arot = 180.D0 ;
      vtra = (0.D0 0.D0 0.D0) ;
      lsym = FAUX ;
   'FINSI' ;    
   'SI' ('EGA' &btype2 2) ;
      arot = 270.D0 ;
      vtra = (0.D0 0.D0 0.D0) ;
      lsym = FAUX ;
   'FINSI' ;    
   nitf = 'DIME' (tabin2 . 'INTERFACE') ;
   append tabin 'STRUCTURE'
      (SYTOPL (tabin2 . 'STRUCTURE') arot vtra lsym) ;
   append tabin 'FUEL'
      (SYTOPL (tabin2 . 'FUEL') arot vtra lsym) ;
   'REPETER' bitf nitf ;
      append (tabin . 'INTERFACE') iitf
         (SYTOPL (tabin2 . 'INTERFACE' . &bitf) arot vtra lsym) ;
      append (tabin . 'CHANNEL') iitf
         (SYTOPL (tabin2 . 'CHANNEL' . &bitf) arot vtra lsym) ;
      iitf = '+' iitf 1 ;
   'FIN' bitf ;
'FIN' btype2 ;   
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPLA2
*
'FINPROC' ;
*ENDPROCEDUR geompla2
*BEGINPROCEDUR geomplh1
************************************************************************
* NOM         : GEOMPLH1
* DESCRIPTION : Géométrie hexagonale à plaques
*               Technical report DM2S/SFME/LTMF/RT/05-050/A
*               typgeo = 1 small losange cell
*               typgeo = 2 normal hexagonal cell
*               typgeo = 3 large two hexagonal cells
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 21/03/2005, version initiale
* HISTORIQUE : v1, 21/03/2005, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPLH1 ;
'ARGUMENT' ept*'FLOTTANT' ;
'ARGUMENT' eps*'FLOTTANT' ;
'ARGUMENT' epp*'FLOTTANT' ;
'ARGUMENT' epg*'FLOTTANT' ;
'ARGUMENT' npp*'ENTIER' ;
'ARGUMENT' typgeo*'ENTIER' ;
*
* On commence par projeter toutes les quantités sur une génératrice
*
ff = '/' 2.D0 ('**' 3.D0 0.5D0) ;
l   = '*' ('/' ept 2.D0) ff ;
eps = '*' eps ff ;
epp = '*' epp ff ;
epg = '*' epg ff ;
eps2 = '/' eps 2.D0 ;
*
res = '-' l (eps2 '+' ('*' epp npp) '+' ('*' epg ('-' npp 1))
                  '+' eps) ;
epg2 = '/' res 2.D0 ;
*
* Discrétisation 
*
ns2 = 1 ;
*ns2 = 2 ;
ns  = '*' ns2 2 ;
ng  = 2 ;
ng2 = 'ENTIER' ('+' ('*' ng ('/' epg2 epg)) 0.5D0) ;
*'MESSAGE' ('CHAINE' 'ng2 = ' ng2) ;
*'MESSAGE' ('CHAINE' 'epg = ' epg) ;
*'MESSAGE' ('CHAINE' 'epg2 = ' epg2) ;
*'MESSAGE' ('CHAINE' 'epg '-' epg2  = ' ('-' epg epg2)) ;
*ng2 = 1 ;
*ng2 = 2 ;
np  = 3 ;
*ns2 = 3 ;
*ns  = 3 ;
*ng  = 3 ;
*ng2 = 3 ;
*np  = 3 ;
*
* Build the horizontal generator
*
itg = 0 ;
tabgen = 'TABLE' ;
bp  = 0.D0 0.D0 0.D0 ;
*  left structure
bpp = 'PLUS'  bp (eps2 0. 0.) ;
ss  = 'DROIT' bp ns2 bpp ;
bp  = bpp ;
itg = '+' itg 1 ;
tabgen . itg = 'TABLE' ;
tabgen . itg . 'TYPE' = 'CHAINE' 'ST' ;
tabgen . itg . 'SEG'  = ss ;
*  left gas channel
bpp = 'PLUS' bp (epg2 0. 0.) ;
ss  = 'DROIT' bp ng2 bpp ;
bp  = bpp ;
itg = '+' itg 1 ;
tabgen . itg = 'TABLE' ;
tabgen . itg . 'TYPE' = 'CHAINE' 'CH' ;
tabgen . itg . 'SEG'  = ss ;
*  plates and channel
'REPETER' ipp npp ;
   bpp = 'PLUS' bp (epp 0. 0.) ;
   ss  = 'DROIT' bp np bpp ;
   bp  = bpp ;
   itg = '+' itg 1 ;
   tabgen . itg = 'TABLE' ;
   tabgen . itg . 'TYPE' = 'CHAINE' 'FU' ;
   tabgen . itg . 'SEG'  = ss ;
   'SI' ('NEG' &ipp npp) ;
      bpp = 'PLUS' bp (epg 0. 0.) ;
      ss  = 'DROIT' bp ng bpp ;
      bp  = bpp ;
      itg = '+' itg 1 ;
      tabgen . itg = 'TABLE' ;
      tabgen . itg . 'TYPE' = 'CHAINE' 'CH' ;
      tabgen . itg . 'SEG'  = ss ;
   'FINSI' ;
'FIN' ipp ;   
*
*  right gas channel
bpp = 'PLUS' bp (epg2 0. 0.) ;
ss  = 'DROIT' bp ng2 bpp ;
bp  = bpp ;
itg = '+' itg 1 ;
tabgen . itg = 'TABLE' ;
tabgen . itg . 'TYPE' = 'CHAINE' 'CH' ;
tabgen . itg . 'SEG'  = ss ;
*  right structure
bpp = 'PLUS'  bp (eps 0. 0.) ;
ss  = 'DROIT' bp ns bpp ;
bp  = bpp ;
itg = '+' itg 1 ;
tabgen . itg = 'TABLE' ;
tabgen . itg . 'TYPE' = 'CHAINE' 'ST' ;
tabgen . itg . 'SEG'  = ss ;
*
* Build the vertical generators
*
prot1 = 0.D0 0.D0 0.D0 ;
prot2 = 0.D0 0.D0 1.D0 ;
arot = 120.D0 ;
tabtmp = 'TABLE' ;
'REPETER' iitg ('DIME' tabgen) ;
   titg = tabgen . &iitg ;
   'SI' ('EGA' (titg . 'TYPE') 'ST') ;
      APPEND tabtmp 'GENSTR' (titg . 'SEG') ;
   'SINON' ;      
      APPEND tabtmp 'GENFU' (titg . 'SEG') ;
   'FINSI' ;
'FIN' iitg ;
genstr  = tabtmp . 'GENSTR' ;   
genfu = tabtmp . 'GENFU' ;
*
genstr  = 'TOURNER' genstr  arot prot1 prot2 ;
genfu   = 'TOURNER' genfu   arot prot1 prot2 ;
gench   = 'DROIT' ('POIN' genfu 'INITIAL') 1 ('POIN' genfu 'FINAL') ;
*gent = 'ET' (tabtmp . 'GENSTR') (tabtmp . 'GENFU') ;
*'LISTE' genstr ;
*'LISTE' genfu ;
*'LISTE' gent ;
*strgen = 'COULEUR' (genstr 'GENERATRICE' gent) 'BLAN' ;
*fugen  = 'COULEUR' (genfu  'GENERATRICE' gent) 'ROUG' ;
*'OPTION' 'DIME' 2 ;
*'TRACER' ('ET' genstr ('COULEUR' gent 'ROUG')) ;
*'TRACER' fugen ;
*'TRACER' (fugen 'ET' strgen) ;
*
* Construire une cellule élémentaire dans tabin1
*
tabin1 = 'TABLE' ;
tabin1 . 'INTERFACE' = 'TABLE' ;
tabin1 . 'CHANNEL'   = 'TABLE' ;
tabin1 . 'DH'        = 'TABLE' ;
ich = 0 ;
bp = 0. 0. 0. ;
'REPETER' iitg ('DIME' tabgen) ;
   titg = tabgen . &iitg ;
   'SI' ('EGA' (titg . 'TYPE') 'ST') ;
      li = titg . 'SEG' ;
      su = ('PLUS' ('ET' genstr genfu) bp) 'GENERATRICE' li ;
      APPEND tabin1 'STRUCTURE' su ;
      bp = 'POIN' li 'FINAL' ;
   'FINSI' ;
   'SI' ('EGA' (titg . 'TYPE') 'CH') ;
      li  = titg . 'SEG' ;
      ct  = 'CONTOUR' (('PLUS' genfu bp) 'GENERATRICE' li) ;
      dh  = '/' ('*' ('MESURE' ct 'SURF') 4.D0) ('MESURE' ct) ;
*      'MESSAGE' ('CHAINE' 'Dh = ' dh) ;
      li2 = 'DROIT' ('POIN' li 'INITIAL') 1 ('POIN' li 'FINAL') ;
      su  = ('PLUS' gench bp) 'GENERATRICE' li2 ;
      ich = '+' ich 1 ;
      APPEND (tabin1 . 'INTERFACE') ich ct ;
      APPEND (tabin1 . 'CHANNEL')   ich su ;
      tabin1 . 'DH' . ich = dh ;
      su  = ('PLUS' genstr bp) 'GENERATRICE' li ;
      APPEND tabin1 'STRUCTURE' su ;
      bp = 'POIN' li 'FINAL' ;
   'FINSI' ;
   'SI' ('EGA' (titg . 'TYPE') 'FU') ;
      li  = titg . 'SEG' ;
      su  = ('PLUS' genfu bp) 'GENERATRICE' li ;
      APPEND tabin1 'FUEL' su ;
      su  = ('PLUS' genstr bp) 'GENERATRICE' li ;
      APPEND tabin1 'STRUCTURE' su ;
      bp = 'POIN' li 'FINAL' ;
   'FINSI' ;
'FIN' iitg ;
*
* Construire un bloc de trois cellules 
*
tabin2 = 'TABLE' ;
tabin2 . 'INTERFACE' = 'TABLE' ;
tabin2 . 'CHANNEL'   = 'TABLE' ;
tabin2 . 'DH'        = 'TABLE' ;
iitf = 0 ;
p0 = 0. 0. 0. ;
pr = 0. 0. 1. ;
'REPETER' bclrot 3 ;
   arot = '*' 120.D0 ('-' &bclrot 1) ;
   append tabin2 'STRUCTURE'
      ('TOURNER' (tabin1 . 'STRUCTURE') arot p0 pr) ;
   append tabin2 'FUEL'
      ('TOURNER' (tabin1 . 'FUEL') arot p0 pr) ;
   'REPETER' bich ich ;
      iitf = '+' iitf 1 ;
      append (tabin2 . 'INTERFACE') iitf
         ('TOURNER' (tabin1 . 'INTERFACE' . &bich) arot p0 pr) ;
      append (tabin2 . 'CHANNEL') iitf
         ('TOURNER' (tabin1 . 'CHANNEL' . &bich) arot p0 pr) ;
      tabin2 . 'DH' . iitf = tabin1 . 'DH' . &bich ;         
   'FIN' bich ;
'FIN' bclrot ;
*
* Construire un bloc de 2 hexagones
* 
tabin3 = 'TABLE' ;
tabin3 . 'INTERFACE' = 'TABLE' ;
tabin3 . 'CHANNEL'   = 'TABLE' ;
tabin3 . 'DH'        = 'TABLE' ;
nitf = iitf ;
iitf = 0 ;
lrot = '/' ('*' l 1.5D0) ('COS' 30.) ;