Télécharger stokes_rima.dgibi

Retour à la liste

* fichier stokes_rima.dgibi
************************************************************************
************************************************************************
*
'OPTI' 'ECHO' 0 ;
*
'SAUTER' 2 'LIGNE' ;
'MESSAGE' ' Execution de stokes_rima.dgibi' ;
'SAUTER' 2 'LIGNE' ;
*
graph = faux ;
interact = faux ;
************************************************************************
* NOM         : STOKES_RIMA
* DESCRIPTION : Cas-test équation de Stokes incompressible
*               Méthode directe et pénalisation des contraintes
*               Cavité carrée entrainee
*               On utilise KRES et RESO avec les MATRIK et les
*               RIGIDITE et on vérifie que les résultats obtenus sont
*               identiques. 
*               Permet de tester également EXTR RIG1 MOT1 MOT2
*                                          EXTR RIG1 'DIAG'
*                                          CMCT
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 10/06/2011, version initiale
* HISTORIQUE : v1, 10/06/2011, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
 'OPTION' 'DIME' 2 'MODE' 'PLAN' 'ISOV' 'SURF' ;
************************************************************************
*
* 
*   PROCEDURES 
*
*
************************************************************************
*
*
*
'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 quafme
************************************************************************
* NOM         : QUAFME
* DESCRIPTION :
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 01/12/2004, version initiale
* HISTORIQUE : v1, 01/12/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' QUAFME ;
'ARGUMENT' tmail/'TABLE' ;
'SI' ('NON' ('EXISTE' tmail)) ;
   tmail = 'TABLE' ;
'FINSI' ;
'REPETER' bcl ;
   'ARGUMENT' mquad/'MAILLAGE' ;
   'SI' ('EXISTE' mquad) ;
      'REPETER' bbcl 1 ;
         'ARGUMENT' nom*'MOT' ;
         tmail . nom = 'TABLE' ;
         dimm = DEADUTIL 'DIMM' mquad ;
         'SI' ('EGA' dimm 0) ;
            mquaf = mquad ;
            tmail . nom = mquad ;
            'QUITTER' bbcl ;
         'FINSI' ;            
         typm = DEADUTIL 'TYPM' mquad ;
         'SI' ('EGA' typm 'QUAF') ;
            mquaf = mquad ;
            tmail . nom . 'QUAF' = mquaf ;
            'SI' ('EGA' dimm 2) ;
               tmail . nom . 'QUAI' = mquaf ;
               tmail . nom . 'QUAD' = mquaf ;
            'FINSI' ;
            'QUITTER' bbcl ;
         'FINSI' ;            
         'SI' ('EGA' typm 'QUAI') ;
            mquaf = 'CHANGER' mquad 'QUAF' ;
            tmail . nom . 'QUAF' = mquaf ;
            tmail . nom . 'QUAI' = mquad ;
            tmail . nom . 'LINE' = 'CHANGER' 'LINEAIRE' mquad ;
            'QUITTER' bbcl ;
         'FINSI' ;            
         'SI' ('EGA' typm 'LINE') ;
            mquai = 'CHANGER' 'QUADRATIQUE' mquad ;
            mquaf = 'CHANGER' mquai 'QUAF' ;
            tmail . nom . 'QUAF' = mquaf ;
            tmail . nom . 'QUAI' = mquai ;
            tmail . nom . 'LINE' = mquad ;
            'QUITTER' bbcl ;
         'FINSI' ;            
      'FIN' bbcl ;
*      ielim = '+' ielim 1 ;
*      telim . ielim = mquaf ;
   'SINON' ;
      'QUITTER' bcl ;
   'FINSI' ;
'FIN' bcl ;
*
*
telim = 'TABLE' 'ESCLAVE' ;
ielim = 0 ;
tidx = 'INDEX' tmail ;
'REPETER' iidx ('DIME' tidx) ;
   idx = tidx . &iidx ;
   val = tmail . idx ;
   tval = 'TYPE' val ;
   'SI' ('EGA' tval 'MAILLAGE') ;
      ielim = '+' ielim 1 ;
      telim . ielim = val ;
   'FINSI' ;      
   'SI' ('EGA' tval 'TABLE') ;
      'SI' ('EXISTE' val 'QUAF') ;
         ielim = '+' ielim 1 ;
         telim . ielim = val . 'QUAF' ;
      'FINSI' ;         
   'FINSI' ;      
'FIN' iidx ;
*
*
'ARGUMENT' tol*'FLOTTANT' ;
'ELIMINATION' ('ET' telim) tol ;
*
'RESPRO' tmail ;
'FINPROC' ;         
*
* End of procedure file QUAFME
*
*ENDPROCEDUR quafme
*BEGINPROCEDUR modenlin
************************************************************************
* NOM         : MODENLIN
* DESCRIPTION :
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, ??/??/2007, version initiale
* HISTORIQUE : v1, ??/??/2007, 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' MODENLIN ;
*
vdim = 'VALEUR' 'DIME' ;
vmod = 'VALEUR' 'MODE' ;
*
'SI' ('ET' ('NEG' vdim 2) ('NEG' vdim 3)) ;
   'ERREUR' ('CHAINE' 'vdim=' vdim) ;
'FINSI' ;   
*
*  Noms par défaut 
*
'SI' ('EGA' vdim 2) ;
   'SI' ('NEG' vmod 'AXIS') ;
      nomvec = 'MOTS' 'UX' 'UY' ;
      nomfor = 'MOTS' 'FX' 'FY' ;
   'SINON' ;      
      nomvec = 'MOTS' 'UR' 'UZ' ;
      nomfor = 'MOTS' 'FR' 'FZ' ;
   'FINSI' ;
'FINSI' ;   
'SI' ('EGA' vdim 3) ;
   nomvec = 'MOTS' 'UX' 'UY' 'UZ' ;
   nomfor = 'MOTS' 'FX' 'FY' 'FZ' ;
'FINSI' ;   
*
nomsca = 'MOTS' 'SCAL' ;
*nomflu = 'MOTS' 'SCAL' ;
*
'ARGUMENT' tmode/'TABLE' ;
'SI' ('NON' ('EXISTE' tmode)) ;
   tmode = 'TABLE' ;
*
* Initialisation des inconnues par défaut
*
   lmdisc = 'MOTS' 'CSTE' 'LINE' 'QUAF' ;
   'REPETER' iidisc ('DIME' lmdisc) ;
      mdisc = 'EXTRAIRE' lmdisc &iidisc ;
      tmode . mdisc = 'TABLE' ;
      tmode . mdisc . 'DISC' = mdisc ;
      tmode . mdisc . 'NOMPRI' = 'TABLE' ;
      tmode . mdisc . 'NOMPRI' . 1 = nomsca ;
      tmode . mdisc . 'NOMDUA' = 'TABLE' ;
      tmode . mdisc . 'NOMDUA' . 1 = nomsca ;
      mdiscv = 'CHAINE' mdisc 'V' ;
      tmode . mdiscv = 'TABLE' ;
      tmode . mdiscv . 'DISC' = mdisc ;
      tmode . mdiscv . 'NOMPRI' = 'TABLE' ;
      'REPETER' idim vdim ;
         TMODE . mdiscv. 'NOMPRI' . &idim =
            'MOTS' ('EXTRAIRE' nomvec &idim) ;
      'FIN' idim ;
      tmode . mdiscv . 'NOMDUA' = 'TABLE' ;
      'REPETER' idim vdim ;
         TMODE . mdiscv. 'NOMDUA' . &idim =
            'MOTS' ('EXTRAIRE' nomfor &idim) ;
      'FIN' idim ;
   'FIN' iidisc ;   
'FINSI' ;
*
* Lecture des mots clés et des inconnues
*
lmotcle = 'MOTS' 'GEOM' 'INCO' ;
ldiscdd = 'MOTS' 'LINM' 'CUBI' ;
ltypinc = 'MOTS' 'SCAL' 'VECT' ;
'REPETER' imotcle ;
   'ARGUMENT' motcle/'MOT' ;
   'SI' ('NON' ('EXISTE' motcle)) ; 'QUITTER' imotcle ; 'FINSI' ;
   'SI' ('NON' ('EXISTE' lmotcle motcle)) ;
      cherr = 'CHAINE' 'Keyword ' motcle ' unknown.' ; 'ERREUR' cherr ;
   'FINSI' ;
   'SI' ('EGA' motcle 'GEOM') ;
      'ARGUMENT' discg*'MOT' ;
      TMODE . 'GEOM' = 'TABLE' ;
      TMODE . 'GEOM' . 'DISC' = discg ;
   'FINSI' ;      
   'SI' ('EGA' motcle 'INCO') ;
      'ARGUMENT' nominc*'MOT' ;
      'SI' ('EGA' nominc 'GEOM') ;
         'ERREUR' 'GEOM nest pas un nom dinconnu acceptable' ;
      'FINSI' ;         
      TMODE . nominc = 'TABLE' ;
      'ARGUMENT' disinc*'MOT' ;
      TMODE . nominc . 'DISC' = disinc ;
      ldd = 'EXISTE' ldiscdd disinc ;
      'ARGUMENT' typinc*'MOT' ;
      'SI' ('NON' ('EXISTE' ltypinc typinc)) ;
         cherr = 'CHAINE' 'Type=' typinc ' unknown.' ;
         'ERREUR' cherr ;
      'FINSI' ;
      'SI' ('EGA' typinc 'SCAL') ;
         nbinc = 1 ;
      'SINON' ;         
         nbinc = vdim ;
      'FINSI' ;         
      TMODE . nominc . 'NOMPRI' = 'TABLE' ;
      'REPETER' ibinc nbinc ;
         'SI' ldd ;
            'ARGUMENT' nomcom*'LISTMOTS' ;
         'SINON' ;
            'ARGUMENT' nomcom*'MOT' ;
            nomcom = 'MOTS' nomcom ;
         'FINSI' ;
         TMODE . nominc . 'NOMPRI' . &ibinc = nomcom ;
      'FIN' ibinc ;
      TMODE . nominc . 'NOMDUA' = 'TABLE' ;
      'REPETER' ibinc nbinc ;
         'SI' ldd ;
            'ARGUMENT' nomcom*'LISTMOTS' ;
         'SINON' ;
            'ARGUMENT' nomcom*'MOT' ;
            nomcom = 'MOTS' nomcom ;
         'FINSI' ;
         TMODE . nominc . 'NOMDUA' . &ibinc = nomcom ;
      'FIN' ibinc ;
   'FINSI' ;      
'FIN' imotcle ;      
*
'RESPRO' tmode ;
*
* End of procedure file MODENLIN
*
'FINPROC' ;
*ENDPROCEDUR modenlin
*BEGINPROCEDUR grig2
************************************************************************
* NOM         : GRIG2
* DESCRIPTION :
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, ??/??/2007, version initiale
* HISTORIQUE : v1, ??/??/2007, 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' GRIG2 ;
'ARGUMENT'     _mt*'MAILLAGE' ;
'ARGUMENT' tdisc*'TABLE' ;
*
* Lectures
*
debug = FAUX ;
lmotcle = 'MOTS' 'NPRI' 'FPRI' 'CPRI' 'NDUA' 'FDUA' 'CDUA'
                 'NCOF' 'FCOF' 'CCOF' 'LAPN' 'GMBT' ;
* Il faut initialiser valt et valq, sinon on peut capturer ceux de
* la procédure appelante
valt  = 'valt' ; valq = 'valq' ;
llapn = 0 ;
'REPETER' imotcle ;
   'ARGUMENT' motcle/'MOT' ;
   'SI' ('NON' ('EXISTE' motcle)) ; 'QUITTER' imotcle ; 'FINSI' ;
   'SI' ('NON' ('EXISTE' lmotcle motcle)) ;
      cherr = 'CHAINE' 'Keyword ' motcle ' unknown.' ; 'ERREUR' cherr ;
   'FINSI' ;
   'SI' ('EGA' motcle 'NPRI') ; 'ARGUMENT' nomt*'MOT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'NDUA') ; 'ARGUMENT' nomq*'MOT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'NCOF') ; 'ARGUMENT' nomo*'MOT' ; 'FINSI' ;
   tst1 = 'EGA' motcle 'FPRI' ; tst2 = 'EGA' motcle 'FDUA' ;
   tst = tst1 'OU' tst2 ;
   'SI' tst ;
      'SI' tst1 ; tt = TDISC . nomt ; 'FINSI' ;
      'SI' tst2 ; tt = TDISC . nomq ; 'FINSI' ;
      isvec = ('>' ('DIME' (tt. 'NOMPRI')) 1) ;
      'SI' isvec ; 'ARGUMENT' valv*'LISTREEL' ; 'SINON' ;
         'ARGUMENT' valv*'FLOTTANT' ;
      'FINSI' ;
      'SI' tst1 ; valt = valv ; 'FINSI' ;
      'SI' tst2 ; valq = valv ; 'FINSI' ;
   'FINSI' ;
   'SI' ('EGA' motcle 'FCOF') ; 'ARGUMENT' valo*'FLOTTANT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'CPRI') ; 'ARGUMENT' valt*'CHPOINT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'CDUA') ; 'ARGUMENT' valq*'CHPOINT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'CCOF') ; 'ARGUMENT' valo*'CHPOINT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'LAPN') ; llapn = 1 ; 'FINSI' ;
   'SI' ('EGA' motcle 'GMBT') ; llapn = 2 ; 'FINSI' ;
'FIN' imotcle ;
*
* Tests
*
discg = TDISC . 'GEOM' . 'DISC' ;
'SI' ('EXISTE' tdisc 'methgau') ;
   methgau = tdisc . 'methgau' . 'rigi' ;
'SINON' ;
   methgau = 'GAU7' ;
'FINSI' ;
tnomt = TDISC . nomt ;
lvalt = 'NEG' ('TYPE' valt) 'MOT' ;
tnomq = TDISC . nomq ;
lvalq = 'NEG' ('TYPE' valq) 'MOT' ;
* Scalaire ou vecteur
ninct = 'DIME' (tnomt . 'NOMPRI') ;
nincq = 'DIME' (tnomq . 'NOMPRI') ;
'SI' ('NEG' ninct nincq) ;
   cherr = 'CHAINE'
   'les primales et duales nont pas le meme nombre de composantes' ;
   'ERREUR' cherr ;
'FINSI' ;
'SI' ('NEG' ninct ('VALEUR' 'DIME')) ;
   cherr = 'CHAINE'
   'les inconnues doivent etre vectorielles' ;
   'ERREUR' cherr ;
'FINSI' ;
*
ninc = ninct ;
*
lcof = 'EXISTE' TDISC nomo ;
'SI' lcof ; ncof = 1 ; tcof = TDISC . nomo ;
'SINON' ; ncof = 0 ;
'FINSI' ;
*
'SI' debug ;
'SI' lcof ; 'MESSAGE' 'Un coef a ete detecte' ;
'SINON' ; 'MESSAGE' 'pas de coef  detecte' ;
'FINSI' ;
'FINSI' ;
*
   vdim = 'VALEUR' 'DIME' ;
   vmod = 'VALEUR' 'MODE' ;
   idim = 0 ;
   'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'PLANDEFO')) ;
      idim = 2 ;
      iaxi = FAUX ;
   'FINSI' ;
   'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'AXIS')) ;
      idim = 2 ;
      iaxi = VRAI ;
   'FINSI' ;
   'SI' ('ET' ('EGA' vdim 3) ('EGA' vmod 'TRID')) ;
      idim = 3 ;
      iaxi = FAUX ;
   'FINSI' ;
   'SI' ('EGA' vdim 1) ;
      idim = 1 ;
      iaxi = FAUX ;
   'FINSI' ;
*   'MESSAGE' ('CHAINE' 'iaxi=' iaxi ); 
   'SI' ('EGA' idim 0) ;
      'ERREUR' ('CHAINE' 'vmod=' vmod ' et vdim=' vdim ' non prevu') ;
   'FINSI' ;
   'SI' iaxi ;
      rmt   = 'COORDONNEE' 1 _mt ;
      deupi = '*' PI 2.D0 ;
   'FINSI' ;
*
* Optimisation possible : construire la matrice par blocs
* qd valt et valq ne sont pas donnés
*
*
*Bug ?   numop  = ('**' ninc 2) '+' 1 ;
   numop  = '**' ninc 2 ;
   'SI' iaxi ; numop = '+' numop 1 ; 'FINSI' ;
   numder = idim ;
   numvar = ninc ;
   ncof = '+' ncof 1 ;
*delete   'SI' iaxi ; ncof = '+' ncof 1 ; 'FINSI' ;
   numdat = ncof ;
   numcof = ncof ;
*
   A = ININLIN numop numvar numdat numcof numder ;
*
   lvt = 'EGA' ('TYPE' valt) 'LISTREEL' ;
   'REPETER' iiinc ninc ;
      iinc = &iiinc ;
      A . 'VAR' . iinc . 'NOMDDL' = tnomt . 'NOMPRI' . iinc ;
      A . 'VAR' . iinc . 'DISC'   = tnomt . 'DISC'  ;
      'SI' lvalt ;
         'SI' lvt ;
            A . 'VAR' . iinc . 'VALEUR'   = 'EXTRAIRE' valt iinc ;
         'SINON' ;
            A . 'VAR' . iinc . 'VALEUR'   = valt ;
         'FINSI' ;            
      'FINSI' ;         
   'FIN' iiinc ;
*
   icof = 0 ;
   ll = 'LECT' ;
   icof = '+' icof 1 ;
   A . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
   A . 'DAT' . icof . 'DISC'   = 'CSTE' ;
   A . 'DAT' . icof . 'VALEUR' = 2.D0 ;
   A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
   A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
   ll  = 'LECT' ;
   ll2 = 'LECT' icof ;
   'SI' lcof ;
      icof = '+' icof 1 ;
      A . 'DAT' . icof . 'NOMDDL' = tcof . 'NOMPRI' . 1 ;
      A . 'DAT' . icof . 'DISC'   = tcof . 'DISC' ;
      A . 'DAT' . icof . 'VALEUR' = valo ;
      A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
      A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
      ll  = 'ET' ll  ('LECT' icof) ; 
      ll2 = 'ET' ll2 ('LECT' icof) ; 
   'FINSI' ;
*
   iop = 0 ;
   'REPETER' iidim idim ;
      'REPETER' jidim idim ;
         iop = '+' iop 1 ;
         'SI' ('EGA' &iidim &jidim) ;
            'SI' ('EGA' llapn 0) ;
               A . iop . &iidim . &jidim = ll2 ;
            'SINON' ;
               A . iop . &iidim . &jidim = ll ;
            'FINSI' ;
         'SINON' ;
            'SI' ('NEG' llapn 2) ;
               A . iop . &iidim . &jidim = ll ;
            'FINSI' ;
            'SI' ('NEG' llapn 1) ;
               A . iop . &jidim . &iidim = ll ;
            'FINSI' ;               
         'FINSI' ;            
      'FIN' jidim ;
   'FIN' iidim ;
   'SI' iaxi ;   
      iop = '+' iop 1 ;
      'SI' ('EGA' llapn 0) ;
         A . iop . 1 . 0 = ll2 ;
      'SINON' ;         
         A . iop . 1 . 0 = ll ;
      'FINSI' ;         
   'FINSI' ;
*
   'SI' iaxi ;
      numdat = 2 ;
      numcof = 2 ;
   'SINON' ;      
      numdat = 0 ;
      numcof = 0 ;
   'FINSI' ;      
*
   B = ININLIN numop numvar numdat numcof numder ;
*
   icof = 0 ;
*   
   'SI' iaxi ;
      icof = '+' icof 1 ;
      B . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
      B . 'DAT' . icof . 'DISC'   = discg ;
      B . 'DAT' . icof . 'VALEUR' = rmt ;
      B . 'COF' . icof . 'COMPOR' = 'IDEN' ;
      B . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
      icof = '+' icof 1 ;
      B . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
      B . 'DAT' . icof . 'DISC'   = 'CSTE' ;
      B . 'DAT' . icof . 'VALEUR' = deupi ;
      B . 'COF' . icof . 'COMPOR' = 'IDEN' ;
      B . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
      ll  = 'LECT' ('-' icof 1) icof ;
      mic = '*' ('-' icof 1) -1 ;
      llb = 'LECT' mic icof ;
   'FINSI' ;      
*
   lvq = 'EGA' ('TYPE' valq) 'LISTREEL' ;
   'REPETER' iiinc ninc ;
      iinc = &iiinc ;
      B . 'VAR' . iinc . 'NOMDDL' = tnomq . 'NOMDUA' . iinc ;
      B . 'VAR' . iinc . 'DISC'   = tnomq . 'DISC'  ;
      'SI' lvalq ;
         'SI' lvq ;
            B . 'VAR' . iinc . 'VALEUR'   = 'EXTRAIRE' valq iinc ;
         'SINON' ;
            B . 'VAR' . iinc . 'VALEUR'   = valq ;
         'FINSI' ;            
      'FINSI' ;         
   'FIN' iiinc ;
*
   iop = 0 ;
   'REPETER' iidim idim ;
      'REPETER' jidim idim ;
         iop = '+' iop 1 ;
         B . iop . &iidim . &jidim = ll ;
      'FIN' jidim ;
   'FIN' iidim ;
   'SI' iaxi ;   
      iop = '+' iop 1 ;
      B . iop . 1 . 0 = llb ;
*      B . iop . 1 . 0 = ll ;
   'FINSI' ;      
*
   mgrig = NLINP discg _mt A B methgau ;
* Integration par parties
*   mgrig = '*' mgrig -1.D0 ;
*
'RESPRO' mgrig ;
'FINPROC' ;
*
* End of procedure file GRIG2
*
*ENDPROCEDUR grig2
*BEGINPROCEDUR gdiv2
************************************************************************
* NOM         : GDIV2
* DESCRIPTION : Une matrice de masse
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v2, 14/03/2006, mise à jour NLIN évolué
* VERSION    : v1, 13/05/2004, version initiale
* HISTORIQUE : v1, 13/05/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' GDIV2 ;
'ARGUMENT'     _mt*'MAILLAGE' ;
'ARGUMENT'     _smt/'MAILLAGE' ;
'ARGUMENT' tdisc*'TABLE' ;
*
* Lectures
*
debug = FAUX ;
lmotcle = 'MOTS' 'NPRI' 'FPRI' 'CPRI' 'NDUA' 'FDUA' 'CDUA'
                 'NCOF' 'FCOF' 'CCOF' 'GBBT' 'GMBT' ;
* Il faut initialiser valt et valq, sinon on peut capturer ceux de
* la procédure appelante
valt = 'valt' ; valq = 'valq' ;
lbbt = 0 ;
*           
'REPETER' imotcle ;
   'ARGUMENT' motcle/'MOT' ;
   'SI' ('NON' ('EXISTE' motcle)) ; 'QUITTER' imotcle ; 'FINSI' ;
   'SI' ('NON' ('EXISTE' lmotcle motcle)) ;
      cherr = 'CHAINE' 'Keyword ' motcle ' unknown.' ; 'ERREUR' cherr ;
   'FINSI' ;
   'SI' ('EGA' motcle 'NPRI') ; 'ARGUMENT' nomt*'MOT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'NDUA') ; 'ARGUMENT' nomq*'MOT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'NCOF') ; 'ARGUMENT' nomo*'MOT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'FPRI') ; 'ARGUMENT' valt*'LISTREEL' ; 'FINSI' ;
   'SI' ('EGA' motcle 'FDUA') ; 'ARGUMENT' valq*'FLOTTANT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'FCOF') ; 'ARGUMENT' valo*'FLOTTANT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'CPRI') ; 'ARGUMENT' valt*'CHPOINT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'CDUA') ; 'ARGUMENT' valq*'CHPOINT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'CCOF') ; 'ARGUMENT' valo*'CHPOINT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'GBBT') ; lbbt = 1 ; 'FINSI' ;
   'SI' ('EGA' motcle 'GMBT') ; lbbt = 2 ; 'FINSI' ;
'FIN' imotcle ;
*
* Tests
*
discg = TDISC . 'GEOM' . 'DISC' ;
'SI' ('EXISTE' tdisc 'methgau') ;
   methgau = tdisc . 'methgau' . 'amor' ;
'SINON' ;
   methgau = 'GAU7' ;
'FINSI' ;
tnomt = TDISC . nomt ;
lvalt = 'NEG' ('TYPE' valt) 'MOT' ;
tnomq = TDISC . nomq ;
lvalq = 'NEG' ('TYPE' valq) 'MOT' ;
*
lcof = 'EXISTE' TDISC nomo ;
'SI' lcof ; ncof = 1 ; tcof = TDISC . nomo ;
'SINON' ; ncof = 0 ;
'FINSI' ;
*
'SI' debug ;
'SI' lcof ; 'MESSAGE' 'Un coef a ete detecte' ;
'SINON' ; 'MESSAGE' 'pas de coef  detecte' ;
'FINSI' ;
'FINSI' ;
*
   vdim = 'VALEUR' 'DIME' ;
   vmod = 'VALEUR' 'MODE' ;
   idim = 0 ;
   'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'PLANDEFO')) ;
      idim = 2 ;
      iaxi = FAUX ;
   'FINSI' ;
   'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'AXIS')) ;
      idim = 2 ;
      iaxi = VRAI ;
   'FINSI' ;
   'SI' ('ET' ('EGA' vdim 3) ('EGA' vmod 'TRID')) ;
      idim = 3 ;
      iaxi = FAUX ;
   'FINSI' ;
   'SI' ('EGA' vdim 1) ;
      idim = 1 ;
      iaxi = FAUX ;
   'FINSI' ;
   'SI' ('EGA' idim 0) ;
      'ERREUR' ('CHAINE' 'vmod=' vmod ' et vdim=' vdim ' non prevu') ;
   'FINSI' ;
   'SI' iaxi ;
      dp = ('*' PI 2.D0) ;
      rmt = 'COORDONNEE' 1 _mt ;
      ncof = ncof '+' 2 ;
   'FINSI' ;
* Scalaire ou vecteur
ninct = 'DIME' (tnomt . 'NOMPRI') ;
nincq = 'DIME' (tnomq . 'NOMPRI') ;
'SI' ('NEG' ninct idim) ;
   cherr = 'CHAINE'
   'la primale doit etre un vecteur' ;
   'ERREUR' cherr ;
'FINSI' ;   
'SI' ('NEG' nincq 1) ;
   cherr = 'CHAINE'
   'la duale doit etre un scalaire' ;
   'ERREUR' cherr ;
'FINSI' ;   
*
   numop = 1 ; numder = idim ; numvar = ninct ;
   numdat = ncof ; numcof = ncof ;
   A = ININLIN numop numvar numdat numcof numder ;
*   
   lvt = 'EGA' ('TYPE' valt) 'LISTREEL' ;
   'REPETER' iiinct ninct ;
      iinct = &iiinct ;
      A . 'VAR' . iinct . 'NOMDDL' = tnomt . 'NOMPRI' . iinct ;
      A . 'VAR' . iinct . 'DISC'   = tnomt . 'DISC'  ;
      'SI' lvalt ;
         'SI' lvt ;
            A . 'VAR' . iinct . 'VALEUR'   = 'EXTRAIRE' valt iinct ;
         'SINON' ;
            A . 'VAR' . iinct . 'VALEUR'   = valt ;
         'FINSI' ;            
      'FINSI' ;         
   'FIN' iiinct ;
*
   icof = 0 ;
   'SI' lcof ;
      icof = '+' icof 1 ;
      A . 'DAT' . icof . 'NOMDDL' = tcof . 'NOMPRI' . 1 ;
      A . 'DAT' . icof . 'DISC'   = tcof . 'DISC' ;
      A . 'DAT' . icof . 'VALEUR' = valo ;
      A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
      A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
      ll = 'LECT' 1 ;
   'SINON' ;
      ll = 'LECT' ;
   'FINSI' ;
*   
   'SI' iaxi ;
      icof = '+' icof 1 ;
      A . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
      A . 'DAT' . icof . 'DISC'   = 'CSTE' ;
      A . 'DAT' . icof . 'VALEUR' = dp ;
      A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
      A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
      icof = '+' icof 1 ;
      A . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
      A . 'DAT' . icof . 'DISC'   = discg ;
      A . 'DAT' . icof . 'VALEUR' = rmt ;
      A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
      A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
      lldpr = ll 'ET' ('LECT' ('-' icof 1) icof) ;
      lldp = ll 'ET' ('LECT' ('-' icof 1)) ;
   'FINSI' ;
*
   'SI' iaxi ;
      'REPETER' iidim idim ;
         A . 1 . &iidim . &iidim  = lldpr ;
      'FIN' iidim ;
      A . 1 . 1 . 0 = lldp ;
   'SINON' ;
      'REPETER' iidim idim ;
         A . 1 . &iidim . &iidim  = ll ;
      'FIN' iidim ;
   'FINSI' ;
*
   numvar = 1 ;
   numdat = 0 ;
   numcof = 0 ;
*
   B = ININLIN numop numvar numdat numcof numder ;
   B . 'VAR' . 1 . 'NOMDDL' = tnomq . 'NOMDUA' . 1 ;
   B . 'VAR' . 1 . 'DISC'   = tnomq . 'DISC'  ;
   'SI' lvalq ;
      B . 'VAR' . 1 . 'VALEUR'   = valq ;
   'FINSI' ;         
   B . 1 . 1 . 0  = 'LECT' ;
*
   'SI' ('OU' ('EGA' lbbt 0) ('EGA' lbbt 1)) ;
      'SI' ('EXISTE' _smt) ;
          mgdiv2 = 'NLIN' discg _mt _smt A B methgau ;
      'SINON' ;       
          mgdiv2 = NLINP discg _mt A B methgau ;
      'FINSI' ;
   'FINSI' ;      
   'SI' ('OU' ('EGA' lbbt 1) ('EGA' lbbt 2)) ;
      B . 'VAR' . 1 . 'NOMDDL' = tnomq . 'NOMPRI' . 1 ;
      'REPETER' iiinct ninct ;
         iinct = &iiinct ;
         A . 'VAR' . iinct . 'NOMDDL' = tnomt . 'NOMDUA' . iinct ;
      'FIN' iiinct ;         
      'SI' ('EXISTE' _smt) ;
          mgdiv3 = 'NLIN' discg _mt _smt B A methgau ;
      'SINON' ;       
          mgdiv3 = NLINP discg _mt B A methgau ;
      'FINSI' ;
   'FINSI' ;
   'SI' ('EGA' lbbt 0) ;
      mgdiv = mgdiv2 ;
   'FINSI' ;      
   'SI' ('EGA' lbbt 1) ;
      mgdiv = mgdiv2 'ET' mgdiv3 ;
   'FINSI' ;      
   'SI' ('EGA' lbbt 2) ;
      mgdiv = mgdiv3 ;
   'FINSI' ;      
'RESPRO' mgdiv ;
'FINPROC' ;
*
* End of procedure file GDIV2
*
*ENDPROCEDUR gdiv2
*BEGINPROCEDUR gnor2
************************************************************************
* NOM         : GNOR2
* DESCRIPTION : Calcule le champ de normales à une surface.
*               Peut servir à calculer une pression, un potentiel
*               lié à la gravité, un volume contenu dans une surface.
*               Attention à l'orientation de la surface !
*
*               Computes a field of normal to a surface.
*               Also useful to compute a pressure field,
*               a gravity potential field, a volume enclosed
*               by a surface.
*               WARNING : The orientation of the surface matters !
*
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 22/04/2011
* HISTORIQUE : v1, 22/04/2011, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
*
*
'DEBPROC' GNOR2 ;
'ARGUMENT'     _mt*'MAILLAGE' ;
'ARGUMENT' tdisc*'TABLE' ;
*
* Lectures
*
dim = 'VALEUR' 'DIME' ;
debug = FAUX ;
lmotcle = 'MOTS' 'NPRI' 'FPRI' 'CPRI' 'NDUA' 'FDUA' 'CDUA'
                 'NCOF' 'FCOF' 'CCOF' ;
* Il faut initialiser valt et valq, sinon on peut capturer ceux de
* la procédure appelante
valt = 'valt' ; valq = 'valq' ;
'REPETER' imotcle ;
   'ARGUMENT' motcle/'MOT' ;
   'SI' ('NON' ('EXISTE' motcle)) ; 'QUITTER' imotcle ; 'FINSI' ;
   'SI' ('NON' ('EXISTE' lmotcle motcle)) ;
      cherr = 'CHAINE' 'Keyword ' motcle ' unknown.' ; 'ERREUR' cherr ;
   'FINSI' ;
   'SI' ('EGA' motcle 'NPRI') ; 'ARGUMENT' nomt*'MOT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'NDUA') ; 'ARGUMENT' nomq*'MOT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'NCOF') ; 'ARGUMENT' nomo*'MOT' ; 'FINSI' ;
   tst1 = 'EGA' motcle 'FPRI' ; tst2 = 'EGA' motcle 'FDUA' ;
   tst3 = 'EGA' motcle 'FCOF' ;
   tst = tst1 'OU' tst2 'OU' tst3 ;
   'SI' tst ;
      'SI' tst1 ; tt = TDISC . nomt ; 'FINSI' ;
      'SI' tst2 ; tt = TDISC . nomq ; 'FINSI' ;
      'SI' tst3 ; tt = TDISC . nomo ; 'FINSI' ;
      isvec = ('>' ('DIME' (tt. 'NOMPRI')) 1) ;
      'SI' isvec ; 'ARGUMENT' valv*'LISTREEL' ; 'SINON' ;
         'ARGUMENT' valv*'FLOTTANT' ;
      'FINSI' ;
      'SI' tst1 ; valt = valv ; 'FINSI' ;
      'SI' tst2 ; valq = valv ; 'FINSI' ;
      'SI' tst3 ; valo = valv ; 'FINSI' ;
   'FINSI' ;
   'SI' ('EGA' motcle 'CPRI') ; 'ARGUMENT' valt*'CHPOINT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'CDUA') ; 'ARGUMENT' valq*'CHPOINT' ; 'FINSI' ;
   'SI' ('EGA' motcle 'CCOF') ; 'ARGUMENT' valo*'CHPOINT' ; 'FINSI' ;
'FIN' imotcle ;
*
* Tests
*
discg = TDISC . 'GEOM' . 'DISC' ;
'SI' ('EXISTE' tdisc 'methgau') ;
   methgau = tdisc . 'methgau' . 'mass' ;
'SINON' ;
   methgau = 'GAU7' ;
'FINSI' ;
tnomt = TDISC . nomt ;
lvalt = 'NEG' ('TYPE' valt) 'MOT' ;
tnomq = TDISC . nomq ;
lvalq = 'NEG' ('TYPE' valq) 'MOT' ;
* Scalaire ou vecteur
ninct = 'DIME' (tnomt . 'NOMPRI') ;
nincq = 'DIME' (tnomq . 'NOMDUA') ;
'SI' ('ET' ('NEG' ninct 1) ('NEG' ninct dim)) ;
   cherr = 'CHAINE'
   'la primale doit etre un scalaire ou un vecteur' ;
   'ERREUR' cherr ;
'FINSI' ;   
'SI' ('NEG' nincq dim) ;
   cherr = 'CHAINE'
   'la duale doit etre un vecteur' ;
   'ERREUR' cherr ;
'FINSI' ;   
*ninc = ninct ;
*
lcof = 'EXISTE' TDISC nomo ;
'SI' lcof ; tcof = TDISC . nomo ;
ncof = 'DIME' (tcof . 'NOMPRI') ;
'SINON' ; ncof = 0 ;
'FINSI' ;
*
'SI' debug ;
'SI' lcof ; 'MESSAGE' 'Un coef a ete detecte' ;
'SINON' ; 'MESSAGE' 'pas de coef  detecte' ;
'FINSI' ;
'FINSI' ;
*
   vdim = 'VALEUR' 'DIME' ;
   vmod = 'VALEUR' 'MODE' ;
   idim = 0 ;
   'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'PLANDEFO')) ;
      idim = 2 ;
      iaxi = FAUX ;
   'FINSI' ;
   'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'AXIS')) ;
      idim = 2 ;
      iaxi = VRAI ;
   'FINSI' ;
   'SI' ('ET' ('EGA' vdim 3) ('EGA' vmod 'TRID')) ;
      idim = 3 ;
      iaxi = FAUX ;
   'FINSI' ;
   'SI' ('EGA' vdim 1) ;
      idim = 1 ;
      iaxi = FAUX ;
   'FINSI' ;
*   'MESSAGE' ('CHAINE' 'iaxi=' iaxi ); 
   'SI' ('EGA' idim 0) ;
      'ERREUR' ('CHAINE' 'vmod=' vmod ' et vdim=' vdim ' non prevu') ;
   'FINSI' ;
   'SI' iaxi ;
      dprmt = '*' ('COORDONNEE' 1 _mt) ('*' PI 2.D0) ;
   'FINSI' ;
*
* Optimisation possible : construire la matrice par blocs
* qd valt et valq ne sont pas donnés
*
   numop = nincq ; numder = idim ; numvar = ninct ;
   numdat = ncof ; numcof = ncof ;
   A = ININLIN numop numvar numdat numcof numder ;
   'SI' lcof ;
      lvo = 'EGA' ('TYPE' valo) 'LISTREEL' ;
      'REPETER' iicof ncof ;
         icof = &iicof ;
         A . 'DAT' . icof . 'NOMDDL' = tcof . 'NOMPRI' . icof ;
         A . 'DAT' . icof . 'DISC'   = tcof . 'DISC' ;
         'SI' lvo ;
            A . 'DAT' . icof . 'VALEUR' = 'EXTRAIRE' valo icof ;
         'SINON' ;            
            A . 'DAT' . icof . 'VALEUR' = valo ;
         'FINSI' ;            
         A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
         A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
      'FIN' iicof ;         
   'SINON' ;
      ll = 'LECT' ;
   'FINSI' ;
   lvt = 'EGA' ('TYPE' valt) 'LISTREEL' ;
   'REPETER' iiincq nincq ;
      iincq = &iiincq ;
      iinct = 'MINIMUM' ('LECT' iincq ninct) ;
      A . 'VAR' . iinct . 'NOMDDL' = tnomt . 'NOMPRI' . iinct ;
      A . 'VAR' . iinct . 'DISC'   = tnomt . 'DISC'  ;
      'SI' lvalt ;
         'SI' lvt ;
            A . 'VAR' . iinct . 'VALEUR'   = 'EXTRAIRE' valt iinct ;
         'SINON' ;
            A . 'VAR' . iinct . 'VALEUR'   = valt ;
         'FINSI' ;            
      'FINSI' ;
      'SI' lcof ;
         icof = 'MINIMUM' ('LECT' iincq ncof) ;
         A . iincq . iinct . 0  = 'LECT' icof ;
      'SINON' ;         
         A . iincq . iinct . 0  = ll ;
      'FINSI' ;         
   'FIN' iiincq ;
*
   'SI' iaxi ;
      numdat = 1 ;
      numcof = dim '+' 1 ;
   'SINON' ;      
      numdat = 0 ;
      numcof = dim ;
   'FINSI' ;      
   numvar = nincq ;
   B = ININLIN numop numvar numdat numcof numder ;
   icof = 0 ;
   'REPETER' iiidim idim ;
      icof = '+' icof 1 ;
      B . 'COF' . icof . 'COMPOR' = 'CHAINE' 'VNOR' &iiidim ;
      B . 'COF' . icof . 'LDAT'   = 'LECT' ;
   'FIN' iiidim ;
*
   'SI' iaxi ;
      icof = '+' icof 1 ;
      B . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
      B . 'DAT' . 1 . 'DISC'   = discg ;
      B . 'DAT' . 1 . 'VALEUR' = dprmt ;
      B . 'COF' . icof . 'COMPOR' = 'IDEN' ;
      B . 'COF' . icof . 'LDAT'   = 'LECT' 1 ;
      ll = 'LECT' icof ;
   'SINON' ;
      ll = 'LECT' ;
   'FINSI' ;
   lvq = 'EGA' ('TYPE' valq) 'LISTREEL' ;
   'REPETER' iiincq nincq ;
      iincq = &iiincq ;
      B . 'VAR' . iincq . 'NOMDDL' = tnomq . 'NOMDUA' . iincq ;
      B . 'VAR' . iincq . 'DISC'   = tnomq . 'DISC'  ;
      'SI' lvalq ;
         'SI' lvq ;
            B . 'VAR' . iincq . 'VALEUR'   = 'EXTRAIRE' valq iincq ;
         'SINON' ;
            B . 'VAR' . iincq . 'VALEUR'   = valq ;
         'FINSI' ;            
      'FINSI' ;         
      B . iincq . iincq . 0  = ('LECT' iincq) 'ET' ll ;
   'FIN' iiincq ;
*
   mgnor2 = NLIN discg _mt A B methgau ;
*
'RESPRO' mgnor2 ;
'FINPROC' ;
*
* End of procedure file GNOR2
*
*ENDPROCEDUR gnor2
*
*
*************************************************************************
*
* 
*   MAIN : 1) MESH
*          2) COMPUTATIONAL LOOP
*          3) TESTs
*
************************************************************************
*
*
* Maillage et modèle
*
*nx = 1 ; ny =2 ;
nx = 20 ; ny= 20 ;
pA = 0. 0. ; pB = 1. 0. ; pC = 1. 1. ; pD = 0. 1. ;
ang = 0. ;
pA pB pC pD = 'TOURNER' pA pB pC pD ang pA ;
discu = 'QUAF' ; discp = 'LINE' ;
'OPTI' 'ELEM' 'QUA4' ;
tol = 1.D-6 ;
bas = 'DROIT' nx pA pB ; 
dro = 'DROIT' ny pB pC ; 
hau = 'DROIT' nx pC pD ; 
gau = 'DROIT' ny pD pA ;
mt = 'DALLER' bas dro hau gau ;
*
tmail = QUAFME mt 'mt' tol ;
tmail = QUAFME tmail bas 'bas' dro 'dro' hau 'hau' gau 'gau' tol ;
_mt = tmail . 'mt' . 'QUAF' ;
mt  = tmail . 'mt' . discu ;
mtp = tmail . 'mt' . discp ;
cmt = 'CONTOUR' mt ;
bas = tmail . 'bas' . discu ;
dro = tmail . 'dro' . discu ;
hau = tmail . 'hau' . discu ;
gau = tmail . 'gau' . discu ;
_bas = tmail . 'bas' . 'QUAF' ;
_dro = tmail . 'dro' . 'QUAF' ;
_hau = tmail . 'hau' . 'QUAF' ;
_gau = tmail . 'gau' . 'QUAF' ;
*
tmode = MODENLIN 'GEOM' 'LINE'
                 'INCO' 'UN' discu 'VECT' 'UX' 'UY' 'FX' 'FY'
                 'INCO' 'PN' discp 'SCAL' 'LXP'  'LXP'
                 ;
nompp = 'EXTRAIRE' (tmode . 'PN' .'NOMPRI' . 1) 1 ;                 
nomvp = 'ET' (tmode . 'UN' .'NOMPRI' . 1)
             (tmode . 'UN' .'NOMPRI' . 2) ; 
nomvd = 'ET' (tmode . 'UN' .'NOMDUA' . 1)
             (tmode . 'UN' .'NOMDUA' . 2) ; 
discg = tmode . 'GEOM' . 'DISC' ;
*
* Matrice de rigidité
*
*rig  = GRIG2 _mt tmode 'NPRI' 'UN' 'NDUA' 'UN' 'LAPN' ;
rig  = GRIG2 _mt tmode 'NPRI' 'UN' 'NDUA' 'UN' ;
*
* Matrice de la contrainte de divergence nulle et sa transposée
*
cnt1 = GDIV2 _mt tmode 'NPRI' 'UN' 'NDUA' 'PN' ;
cnt2 = GDIV2 _mt tmode 'NPRI' 'UN' 'NDUA' 'PN' 'GMBT' ;
mcnt = cnt1 'ET' cnt2 ;
*
* Conditions aux limites
*
* On relache la contrainte sur la vitesse normale en 1 point (pmil) pour
* fixer le mode de pression constante
pmil = * ('PLUS' pC pD) 0.5 ;
pmil = 'POIN' hau 'PROC' pmil ;
haur = 'DIFF' ('CHANGER' 'POI1' hau)
              ('MANUEL' 'POI1' pmil) ;
gaur = 'DIFF' ('CHANGER' 'POI1' gau)
              ('MANUEL' 'POI1' pA) ;
dror = 'DIFF' ('CHANGER' 'POI1' dro)
               ('MANUEL' 'POI1' pB) ;
xbas ybas = 'COORDONNEE' bas ;
sbas = '**' ('+' ('**' xbas 2.) ('**' ybas 2)) 0.5 ;
sbas = '**' ('SIN' ('*' sbas 180.)) 2. ;
ssbas = 'SIN' ('*' 90. sbas) ;
vAB = 'MOIN' pB pA ;
xab = 'COORDONNEE' 1 vAB ;
yab = 'COORDONNEE' 2 vAB ;
*
ssbas = '+' ('NOMC' 'UX' ('*' ssbas xab))
            ('NOMC' 'UY' ('*' ssbas yab)) ;
*
* Construction de la normale
*
vnorh  = GNOR2 _hau tmode 'NPRI' discg 'FPRI' 1. 'NDUA' 'UN' ;
vnorg  = GNOR2 _gau tmode 'NPRI' discg 'FPRI' 1. 'NDUA' 'UN' ;
vnord  = GNOR2 _dro tmode 'NPRI' discg 'FPRI' 1. 'NDUA' 'UN' ;
vnorh = 'NOMC' nomvd nomvp vnorh ;
vnorg = 'NOMC' nomvd nomvp vnorg ;
vnord = 'NOMC' nomvd nomvp vnord ;
*
mblh = 'BLOQUE' 'DEPL' 'DIRE' vnorh haur ;
mblg = 'BLOQUE' 'DEPL' 'DIRE' vnorg gaur ;
mbld = 'BLOQUE' 'DEPL' 'DIRE' vnord dror ;
mblb = 'BLOQUE' 'DEPL' bas ;
*
bltot = mblh 'ET' mblg 'ET' mbld 'ET' mblb ;
*bltot =  mblg 'ET' mbld 'ET' mblb ;
fbltot = 'DEPIMPOSE' mblb ssbas ;
*
* Matrice et second membre totaux
* Résolution
*
mtot =  rig 'ET' mcnt 'ET' bltot ;
ftot = fbltot ;
sol  = KRES  mtot ftot 'IMPR' 2 ;
tres = 'TABLE'  ;
*
* Lagrangien augmenté pour la matrice totale
*
'REPETER' imatrik 2 ;
   iimatrik=&imatrik ;
*   iimatrik = 1 ;
   lmatrik = 'EGA' iimatrik 1 ;
*  On récupère les sous-matrices à partir de la matrice
*  totale pour tester l'extraction
   'SI' lmatrik ;
      kmtot = 'KOPS' nincdupr ('KOPS' rima mtot) ;
      krig  = 'EXTRAIRE' kmtot nomvp nomvp ;
      kcnt2 = 'EXTRAIRE' kmtot ('MOTS' 'LXP')
                               nomvp ;
      kbltot = 'KOPS' nincdupr ('KOPS' rima bltot) ;
      kcntbl = 'EXTR' kbltot    ('MOTS' 'LX')
                                nomvp ;
      kfbltot = 'KOPS' nincdupr fbltot ;                             
   'SINON' ;
      kmtot = mtot ;
      krig = 'EXTRAIRE' kmtot nomvp nomvd ;
      kcnt2 = 'EXTRAIRE' kmtot ('MOTS' 'LXP')
                               nomvd ;
      kbltot = bltot ;
      kcntbl = 'EXTRAIRE' kbltot ('MOTS' 'LX')
                              nomvd ;
      kfbltot = fbltot ;                                 
   'FINSI' ;                                
*
   mlagd = 'KOPS' 'CMCT' kcnt2 kcnt2 ;
   mlagb  = 'KOPS' 'CMCT' kcntbl kcntbl ;
   flagb = '*' kbltot ('NOMC' ('MOTS' 'FLX')
                          ('MOTS' 'LX') kfbltot) ;
* Scaling
   alphad = 1.D4 ;
   alphab = 1.D8 ;
   drig = 'EXTR' krig 'DIAG' ;
   arig  = 'MAXIMUM' drig ;
   dlagd = 'EXTR' mlagd 'DIAG' ;
   alagd = 'MAXIMUM' dlagd ;
   dlagb = 'EXTRAIRE' mlagb 'DIAG' ;
   alagb = 'MAXIMUM' dlagb ;
   tres . ('CHAINE' 'drig' iimatrik) = drig ;
   tres . ('CHAINE' 'dlagd' iimatrik) = dlagd ;
   tres . ('CHAINE' 'dlagb' iimatrik) = dlagb ;
   alpd  = '*' ('/' arig alagd) alphad ;
   alpb  = '*' ('/' arig alagb) alphab ;
   kkmtot = krig 'ET' ('*' mlagd alpd)
      'ET' ('*' mlagb alpb) ;
   kkftot = ('*' flagb alpb) ;
   'SI' lmatrik ;
      sol2 = KRES  kkmtot kkftot 'IMPR' 2 'SCALING' 2 ;
   'SINON' ;      
      sol2 = 'RESOUD' kkmtot kkftot ;
   'FINSI' ;   
   tres . ('CHAINE' 'sol' iimatrik) = sol2 ;
'FIN' imatrik ;
*
* Post traitement
*
'SI' graph ;
   vit = 'EXCO' ('MOTS' 'UX' 'UY') sol ;
   vvit = 'VECT' vit 'DEPL' 'JAUN' ;
   vit2 = 'EXCO' ('MOTS' 'UX' 'UY') (tres . 'sol1') ;
   vvit2 = 'VECT' vit2 'DEPL' 'ROUG' ;
   vit3 = 'EXCO' ('MOTS' 'UX' 'UY') (tres . 'sol2') ;
   vvit3 = 'VECT' vit3 'DEPL' 'TURQ' ;
   vvitt = vvit 'ET' vvit2 'ET' vvit3 ;
   'TRACER' vvitt cmt 'TITR' 'Vitesses' ;
'FINSI' ;
*
* Tests
*
'SAUTER' 2 'LIGNE' ;
lpass = vrai ;
drig1 = tres . 'drig1' ;
drig2 = tres . 'drig2' ;
dlagb1 = tres . 'dlagb1' ;
dlagb2 = tres . 'dlagb2' ;
dlagd1 = tres . 'dlagd1' ;
dlagd2 = tres . 'dlagd2' ;
errv = 1.D-13 ;
err1 = '/' ('MAXIMUM' ('-' drig1 drig2) 'ABS')
           ('MAXIMUM' drig1 'ABS') ;
err2 = '/' ('MAXIMUM' ('-' dlagb1 dlagb2) 'ABS')
           ('MAXIMUM' dlagb1 'ABS') ;
err3 = '/' ('MAXIMUM' ('-' dlagd1 dlagd2) 'ABS')
           ('MAXIMUM' dlagd1 'ABS') ;
tst1 = '<' err1 errv ;
tst2 = '<' err2 errv ;
tst3 = '<' err3 errv ;
tst = tst1 'ET' tst2 'ET' tst3 ;
'MESSAGE' ('CHAINE' 'Test 1 sur la diagonale des matrices')   ;
'MESSAGE' '         err1 = ' err1  ;
'MESSAGE' '         err2 = ' err2  ;
'MESSAGE' '         err3 = ' err3  ;
'SI' tst ;
   'MESSAGE' 'Test 1 OK' ;
'SINON' ;   
   'MESSAGE' '!!! Test 1 not passed  ' ;
'FINSI' ;
lpass = lpass 'ET' tst ;
*
sol1 = tres . 'sol1' ;
sol2 = tres . 'sol2' ;
* La pénalisation provoque des erreurs d'arrondis
* suivant le scaling, le solveur,etc...
errv = 1.D-13 '*' alphad ;
err1 = '/' ('MAXIMUM' ('-' sol1 sol2) 'ABS')
           ('MAXIMUM' sol1 'ABS') ;
tst = '<' err1 errv ;
'MESSAGE' ('CHAINE' 'Test 2 MATRIK/RIGIDITE')   ;
'MESSAGE' '         err1 = ' err1  ;
'SI' tst ;
   'MESSAGE' 'Test 2 OK' ;
'SINON' ;   
   'MESSAGE' '!!! Test 2 not passed  ' ;
'FINSI' ;
lpass = lpass 'ET' tst ;
*
* La pénalisation donne des erreurs par rapport
* à la méthode directe
*
errv = '/' 1.D0 alphad ;
err1 = '/' ('MAXIMUM' ('-' sol sol2) 'ABS' 'AVEC' nomvp)
           ('MAXIMUM' sol 'ABS' 'AVEC' nomvp) ;
tst = '<' err1 errv ;
'MESSAGE' ('CHAINE' 'Test 3 Direct/Penalisation')   ;
'MESSAGE' '         err1 = ' err1  ;
'SI' tst ;
   'MESSAGE' 'Test 3 OK' ;
'SINON' ;   
   'MESSAGE' '!!! Test 3 not passed  ' ;
'FINSI' ;
lpass = lpass 'ET' tst ;
*
 
 
'SAUTER' 2 'LIGNE' ;
'SI' lpass ;
   'MESSAGE' 'Tout sest bien passe' ;
'SINON' ;
   'MESSAGE' 'Il y a eu des erreurs' ;
'FINSI' ;   
'SAUTER' 2 'LIGNE' ;
'SI' ('NON' lpass) ;
   'ERREUR' 5 ;
'FINSI' ;   
************************************************************************
*
* END OF TEST PART
* 
************************************************************************
'FINSI' ;
'SI' interact ;
   'OPTION' 'DONN' 5 'ECHO' 1 ;
'FINSI' ;   
*
 'FIN' ; 
 
 
 
 
 
 
 
 
 
 
 
 

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