Télécharger la plaquette Cast3M


fr en

Télécharger stokes_lagaug.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier stokes_lagaug.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. 'OPTI' 'ECHO' 0 ;
  6. *
  7. 'SAUTER' 2 'LIGNE' ;
  8. 'MESSAGE' ' Execution de stokes_lagaug.dgibi' ;
  9. 'SAUTER' 2 'LIGNE' ;
  10. *
  11. graph = faux ;
  12. interact = faux ;
  13. ************************************************************************
  14. * NOM         : STOKES_LAGAUG
  15. * DESCRIPTION : Cas-test équation de Stokes incompressible
  16. *               Méthode directe et de Lagrangien augmenté pour les
  17. *               contraintes de divergence.
  18. *               Cavité carrée entrainee
  19. *
  20. *
  21. * LANGAGE     : GIBIANE-CAST3M
  22. * AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  23. *               mél : gounand@semt2.smts.cea.fr
  24. **********************************************************************
  25. * VERSION    : v1, 10/06/2011, version initiale
  26. * HISTORIQUE : v1, 10/06/2011, création
  27. * HISTORIQUE :
  28. * HISTORIQUE :
  29. ************************************************************************
  30.  'OPTION' 'DIME' 2 'MODE' 'PLAN' 'ISOV' 'SURF' ;
  31. ************************************************************************
  32. *
  33. * 
  34. *   PROCEDURES 
  35. *
  36. *
  37. ************************************************************************
  38. *
  39. *
  40. *
  41. 'DEBPROC' ERRREL ;
  42. 'ARGUMENT' val*'FLOTTANT' ;
  43. 'ARGUMENT' valref*'FLOTTANT' ;
  44. *
  45. 'SI' ('<' ('ABS' valref) 1.D-10) ;
  46.    echref = 1.D0 ;
  47. 'SINON' ;
  48.    echref = valref ;   
  49. 'FINSI' ;
  50. *
  51. errabs = 'ABS' ('/' ('-' val valref) echref);
  52. *
  53. 'RESPRO' errabs ;
  54. *
  55. * End of procedure file ERRREL
  56. *
  57. 'FINPROC' ;
  58. *ENDPROCEDUR errrel
  59. *BEGINPROCEDUR quafme
  60. ************************************************************************
  61. * NOM         : QUAFME
  62. * DESCRIPTION :
  63. *
  64. *
  65. *
  66. * LANGAGE     : GIBIANE-CAST3M
  67. * AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  68. *               mél : gounand@semt2.smts.cea.fr
  69. **********************************************************************
  70. * VERSION    : v1, 01/12/2004, version initiale
  71. * HISTORIQUE : v1, 01/12/2004, création
  72. * HISTORIQUE :
  73. * HISTORIQUE :
  74. ************************************************************************
  75. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  76. * en cas de modification de ce sous-programme afin de faciliter
  77. * la maintenance !
  78. ************************************************************************
  79. *
  80. *
  81. 'DEBPROC' QUAFME ;
  82. 'ARGUMENT' tmail/'TABLE' ;
  83. 'SI' ('NON' ('EXISTE' tmail)) ;
  84.    tmail = 'TABLE' ;
  85. 'FINSI' ;
  86. 'REPETER' bcl ;
  87.    'ARGUMENT' mquad/'MAILLAGE' ;
  88.    'SI' ('EXISTE' mquad) ;
  89.       'REPETER' bbcl 1 ;
  90.          'ARGUMENT' nom*'MOT' ;
  91.          tmail . nom = 'TABLE' ;
  92.          dimm = DEADUTIL 'DIMM' mquad ;
  93.          'SI' ('EGA' dimm 0) ;
  94.             mquaf = mquad ;
  95.             tmail . nom = mquad ;
  96.             'QUITTER' bbcl ;
  97.          'FINSI' ;            
  98.          typm = DEADUTIL 'TYPM' mquad ;
  99.          'SI' ('EGA' typm 'QUAF') ;
  100.             mquaf = mquad ;
  101.             tmail . nom . 'QUAF' = mquaf ;
  102.             'SI' ('EGA' dimm 2) ;
  103.                tmail . nom . 'QUAI' = mquaf ;
  104.                tmail . nom . 'QUAD' = mquaf ;
  105.             'FINSI' ;
  106.             'QUITTER' bbcl ;
  107.          'FINSI' ;            
  108.          'SI' ('EGA' typm 'QUAI') ;
  109.             mquaf = 'CHANGER' mquad 'QUAF' ;
  110.             tmail . nom . 'QUAF' = mquaf ;
  111.             tmail . nom . 'QUAI' = mquad ;
  112.             tmail . nom . 'LINE' = 'CHANGER' 'LINEAIRE' mquad ;
  113.             'QUITTER' bbcl ;
  114.          'FINSI' ;            
  115.          'SI' ('EGA' typm 'LINE') ;
  116.             mquai = 'CHANGER' 'QUADRATIQUE' mquad ;
  117.             mquaf = 'CHANGER' mquai 'QUAF' ;
  118.             tmail . nom . 'QUAF' = mquaf ;
  119.             tmail . nom . 'QUAI' = mquai ;
  120.             tmail . nom . 'LINE' = mquad ;
  121.             'QUITTER' bbcl ;
  122.          'FINSI' ;            
  123.       'FIN' bbcl ;
  124. *      ielim = '+' ielim 1 ;
  125. *      telim . ielim = mquaf ;
  126.    'SINON' ;
  127.       'QUITTER' bcl ;
  128.    'FINSI' ;
  129. 'FIN' bcl ;
  130. *
  131. *
  132. telim = 'TABLE' 'ESCLAVE' ;
  133. ielim = 0 ;
  134. tidx = 'INDEX' tmail ;
  135. 'REPETER' iidx ('DIME' tidx) ;
  136.    idx = tidx . &iidx ;
  137.    val = tmail . idx ;
  138.    tval = 'TYPE' val ;
  139.    'SI' ('EGA' tval 'MAILLAGE') ;
  140.       ielim = '+' ielim 1 ;
  141.       telim . ielim = val ;
  142.    'FINSI' ;      
  143.    'SI' ('EGA' tval 'TABLE') ;
  144.       'SI' ('EXISTE' val 'QUAF') ;
  145.          ielim = '+' ielim 1 ;
  146.          telim . ielim = val . 'QUAF' ;
  147.       'FINSI' ;         
  148.    'FINSI' ;      
  149. 'FIN' iidx ;
  150. *
  151. *
  152. 'ARGUMENT' tol*'FLOTTANT' ;
  153. 'ELIMINATION' ('ET' telim) tol ;
  154. *
  155. 'RESPRO' tmail ;
  156. 'FINPROC' ;         
  157. *
  158. * End of procedure file QUAFME
  159. *
  160. *ENDPROCEDUR quafme
  161. *BEGINPROCEDUR modenlin
  162. ************************************************************************
  163. * NOM         : MODENLIN
  164. * DESCRIPTION :
  165. *
  166. *
  167. *
  168. * LANGAGE     : GIBIANE-CAST3M
  169. * AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  170. *               mél : gounand@semt2.smts.cea.fr
  171. **********************************************************************
  172. * VERSION    : v1, ??/??/2007, version initiale
  173. * HISTORIQUE : v1, ??/??/2007, création
  174. * HISTORIQUE :
  175. * HISTORIQUE :
  176. ************************************************************************
  177. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  178. * en cas de modification de ce sous-programme afin de faciliter
  179. * la maintenance !
  180. ************************************************************************
  181. *
  182. *
  183. 'DEBPROC' MODENLIN ;
  184. *
  185. vdim = 'VALEUR' 'DIME' ;
  186. vmod = 'VALEUR' 'MODE' ;
  187. *
  188. 'SI' ('ET' ('NEG' vdim 2) ('NEG' vdim 3)) ;
  189.    'ERREUR' ('CHAINE' 'vdim=' vdim) ;
  190. 'FINSI' ;   
  191. *
  192. *  Noms par défaut 
  193. *
  194. 'SI' ('EGA' vdim 2) ;
  195.    'SI' ('NEG' vmod 'AXIS') ;
  196.       nomvec = 'MOTS' 'UX' 'UY' ;
  197.       nomfor = 'MOTS' 'FX' 'FY' ;
  198.    'SINON' ;      
  199.       nomvec = 'MOTS' 'UR' 'UZ' ;
  200.       nomfor = 'MOTS' 'FR' 'FZ' ;
  201.    'FINSI' ;
  202. 'FINSI' ;   
  203. 'SI' ('EGA' vdim 3) ;
  204.    nomvec = 'MOTS' 'UX' 'UY' 'UZ' ;
  205.    nomfor = 'MOTS' 'FX' 'FY' 'FZ' ;
  206. 'FINSI' ;   
  207. *
  208. nomsca = 'MOTS' 'SCAL' ;
  209. *nomflu = 'MOTS' 'SCAL' ;
  210. *
  211. 'ARGUMENT' tmode/'TABLE' ;
  212. 'SI' ('NON' ('EXISTE' tmode)) ;
  213.    tmode = 'TABLE' ;
  214. *
  215. * Initialisation des inconnues par défaut
  216. *
  217.    lmdisc = 'MOTS' 'CSTE' 'LINE' 'QUAF' ;
  218.    'REPETER' iidisc ('DIME' lmdisc) ;
  219.       mdisc = 'EXTRAIRE' lmdisc &iidisc ;
  220.       tmode . mdisc = 'TABLE' ;
  221.       tmode . mdisc . 'DISC' = mdisc ;
  222.       tmode . mdisc . 'NOMPRI' = 'TABLE' ;
  223.       tmode . mdisc . 'NOMPRI' . 1 = nomsca ;
  224.       tmode . mdisc . 'NOMDUA' = 'TABLE' ;
  225.       tmode . mdisc . 'NOMDUA' . 1 = nomsca ;
  226.       mdiscv = 'CHAINE' mdisc 'V' ;
  227.       tmode . mdiscv = 'TABLE' ;
  228.       tmode . mdiscv . 'DISC' = mdisc ;
  229.       tmode . mdiscv . 'NOMPRI' = 'TABLE' ;
  230.       'REPETER' idim vdim ;
  231.          TMODE . mdiscv. 'NOMPRI' . &idim =
  232.             'MOTS' ('EXTRAIRE' nomvec &idim) ;
  233.       'FIN' idim ;
  234.       tmode . mdiscv . 'NOMDUA' = 'TABLE' ;
  235.       'REPETER' idim vdim ;
  236.          TMODE . mdiscv. 'NOMDUA' . &idim =
  237.             'MOTS' ('EXTRAIRE' nomfor &idim) ;
  238.       'FIN' idim ;
  239.    'FIN' iidisc ;   
  240. 'FINSI' ;
  241. *
  242. * Lecture des mots clés et des inconnues
  243. *
  244. lmotcle = 'MOTS' 'GEOM' 'INCO' ;
  245. ldiscdd = 'MOTS' 'LINM' 'CUBI' ;
  246. ltypinc = 'MOTS' 'SCAL' 'VECT' ;
  247. 'REPETER' imotcle ;
  248.    'ARGUMENT' motcle/'MOT' ;
  249.    'SI' ('NON' ('EXISTE' motcle)) ; 'QUITTER' imotcle ; 'FINSI' ;
  250.    'SI' ('NON' ('EXISTE' lmotcle motcle)) ;
  251.       cherr = 'CHAINE' 'Keyword ' motcle ' unknown.' ; 'ERREUR' cherr ;
  252.    'FINSI' ;
  253.    'SI' ('EGA' motcle 'GEOM') ;
  254.       'ARGUMENT' discg*'MOT' ;
  255.       TMODE . 'GEOM' = 'TABLE' ;
  256.       TMODE . 'GEOM' . 'DISC' = discg ;
  257.    'FINSI' ;      
  258.    'SI' ('EGA' motcle 'INCO') ;
  259.       'ARGUMENT' nominc*'MOT' ;
  260.       'SI' ('EGA' nominc 'GEOM') ;
  261.          'ERREUR' 'GEOM nest pas un nom dinconnu acceptable' ;
  262.       'FINSI' ;         
  263.       TMODE . nominc = 'TABLE' ;
  264.       'ARGUMENT' disinc*'MOT' ;
  265.       TMODE . nominc . 'DISC' = disinc ;
  266.       ldd = 'EXISTE' ldiscdd disinc ;
  267.       'ARGUMENT' typinc*'MOT' ;
  268.       'SI' ('NON' ('EXISTE' ltypinc typinc)) ;
  269.          cherr = 'CHAINE' 'Type=' typinc ' unknown.' ;
  270.          'ERREUR' cherr ;
  271.       'FINSI' ;
  272.       'SI' ('EGA' typinc 'SCAL') ;
  273.          nbinc = 1 ;
  274.       'SINON' ;         
  275.          nbinc = vdim ;
  276.       'FINSI' ;         
  277.       TMODE . nominc . 'NOMPRI' = 'TABLE' ;
  278.       'REPETER' ibinc nbinc ;
  279.          'SI' ldd ;
  280.             'ARGUMENT' nomcom*'LISTMOTS' ;
  281.          'SINON' ;
  282.             'ARGUMENT' nomcom*'MOT' ;
  283.             nomcom = 'MOTS' nomcom ;
  284.          'FINSI' ;
  285.          TMODE . nominc . 'NOMPRI' . &ibinc = nomcom ;
  286.       'FIN' ibinc ;
  287.       TMODE . nominc . 'NOMDUA' = 'TABLE' ;
  288.       'REPETER' ibinc nbinc ;
  289.          'SI' ldd ;
  290.             'ARGUMENT' nomcom*'LISTMOTS' ;
  291.          'SINON' ;
  292.             'ARGUMENT' nomcom*'MOT' ;
  293.             nomcom = 'MOTS' nomcom ;
  294.          'FINSI' ;
  295.          TMODE . nominc . 'NOMDUA' . &ibinc = nomcom ;
  296.       'FIN' ibinc ;
  297.    'FINSI' ;      
  298. 'FIN' imotcle ;      
  299. *
  300. 'RESPRO' tmode ;
  301. *
  302. * End of procedure file MODENLIN
  303. *
  304. 'FINPROC' ;
  305. *ENDPROCEDUR modenlin
  306. *BEGINPROCEDUR grig2
  307. ************************************************************************
  308. * NOM         : GRIG2
  309. * DESCRIPTION :
  310. *
  311. *
  312. *
  313. * LANGAGE     : GIBIANE-CAST3M
  314. * AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  315. *               mél : gounand@semt2.smts.cea.fr
  316. **********************************************************************
  317. * VERSION    : v1, ??/??/2007, version initiale
  318. * HISTORIQUE : v1, ??/??/2007, création
  319. * HISTORIQUE :
  320. * HISTORIQUE :
  321. ************************************************************************
  322. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  323. * en cas de modification de ce sous-programme afin de faciliter
  324. * la maintenance !
  325. ************************************************************************
  326. *
  327. *
  328. 'DEBPROC' GRIG2 ;
  329. 'ARGUMENT'     _mt*'MAILLAGE' ;
  330. 'ARGUMENT' tdisc*'TABLE' ;
  331. *
  332. * Lectures
  333. *
  334. debug = FAUX ;
  335. lmotcle = 'MOTS' 'NPRI' 'FPRI' 'CPRI' 'NDUA' 'FDUA' 'CDUA'
  336.                  'NCOF' 'FCOF' 'CCOF' 'LAPN' 'GMBT' ;
  337. * Il faut initialiser valt et valq, sinon on peut capturer ceux de
  338. * la procédure appelante
  339. valt  = 'valt' ; valq = 'valq' ;
  340. llapn = 0 ;
  341. 'REPETER' imotcle ;
  342.    'ARGUMENT' motcle/'MOT' ;
  343.    'SI' ('NON' ('EXISTE' motcle)) ; 'QUITTER' imotcle ; 'FINSI' ;
  344.    'SI' ('NON' ('EXISTE' lmotcle motcle)) ;
  345.       cherr = 'CHAINE' 'Keyword ' motcle ' unknown.' ; 'ERREUR' cherr ;
  346.    'FINSI' ;
  347.    'SI' ('EGA' motcle 'NPRI') ; 'ARGUMENT' nomt*'MOT' ; 'FINSI' ;
  348.    'SI' ('EGA' motcle 'NDUA') ; 'ARGUMENT' nomq*'MOT' ; 'FINSI' ;
  349.    'SI' ('EGA' motcle 'NCOF') ; 'ARGUMENT' nomo*'MOT' ; 'FINSI' ;
  350.    tst1 = 'EGA' motcle 'FPRI' ; tst2 = 'EGA' motcle 'FDUA' ;
  351.    tst = tst1 'OU' tst2 ;
  352.    'SI' tst ;
  353.       'SI' tst1 ; tt = TDISC . nomt ; 'FINSI' ;
  354.       'SI' tst2 ; tt = TDISC . nomq ; 'FINSI' ;
  355.       isvec = ('>' ('DIME' (tt. 'NOMPRI')) 1) ;
  356.       'SI' isvec ; 'ARGUMENT' valv*'LISTREEL' ; 'SINON' ;
  357.          'ARGUMENT' valv*'FLOTTANT' ;
  358.       'FINSI' ;
  359.       'SI' tst1 ; valt = valv ; 'FINSI' ;
  360.       'SI' tst2 ; valq = valv ; 'FINSI' ;
  361.    'FINSI' ;
  362.    'SI' ('EGA' motcle 'FCOF') ; 'ARGUMENT' valo*'FLOTTANT' ; 'FINSI' ;
  363.    'SI' ('EGA' motcle 'CPRI') ; 'ARGUMENT' valt*'CHPOINT' ; 'FINSI' ;
  364.    'SI' ('EGA' motcle 'CDUA') ; 'ARGUMENT' valq*'CHPOINT' ; 'FINSI' ;
  365.    'SI' ('EGA' motcle 'CCOF') ; 'ARGUMENT' valo*'CHPOINT' ; 'FINSI' ;
  366.    'SI' ('EGA' motcle 'LAPN') ; llapn = 1 ; 'FINSI' ;
  367.    'SI' ('EGA' motcle 'GMBT') ; llapn = 2 ; 'FINSI' ;
  368. 'FIN' imotcle ;
  369. *
  370. * Tests
  371. *
  372. discg = TDISC . 'GEOM' . 'DISC' ;
  373. 'SI' ('EXISTE' tdisc 'methgau') ;
  374.    methgau = tdisc . 'methgau' . 'rigi' ;
  375. 'SINON' ;
  376.    methgau = 'GAU7' ;
  377. 'FINSI' ;
  378. tnomt = TDISC . nomt ;
  379. lvalt = 'NEG' ('TYPE' valt) 'MOT' ;
  380. tnomq = TDISC . nomq ;
  381. lvalq = 'NEG' ('TYPE' valq) 'MOT' ;
  382. * Scalaire ou vecteur
  383. ninct = 'DIME' (tnomt . 'NOMPRI') ;
  384. nincq = 'DIME' (tnomq . 'NOMPRI') ;
  385. 'SI' ('NEG' ninct nincq) ;
  386.    cherr = 'CHAINE'
  387.    'les primales et duales nont pas le meme nombre de composantes' ;
  388.    'ERREUR' cherr ;
  389. 'FINSI' ;
  390. 'SI' ('NEG' ninct ('VALEUR' 'DIME')) ;
  391.    cherr = 'CHAINE'
  392.    'les inconnues doivent etre vectorielles' ;
  393.    'ERREUR' cherr ;
  394. 'FINSI' ;
  395. *
  396. ninc = ninct ;
  397. *
  398. lcof = 'EXISTE' TDISC nomo ;
  399. 'SI' lcof ; ncof = 1 ; tcof = TDISC . nomo ;
  400. 'SINON' ; ncof = 0 ;
  401. 'FINSI' ;
  402. *
  403. 'SI' debug ;
  404. 'SI' lcof ; 'MESSAGE' 'Un coef a ete detecte' ;
  405. 'SINON' ; 'MESSAGE' 'pas de coef  detecte' ;
  406. 'FINSI' ;
  407. 'FINSI' ;
  408. *
  409.    vdim = 'VALEUR' 'DIME' ;
  410.    vmod = 'VALEUR' 'MODE' ;
  411.    idim = 0 ;
  412.    'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'PLANDEFO')) ;
  413.       idim = 2 ;
  414.       iaxi = FAUX ;
  415.    'FINSI' ;
  416.    'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'AXIS')) ;
  417.       idim = 2 ;
  418.       iaxi = VRAI ;
  419.    'FINSI' ;
  420.    'SI' ('ET' ('EGA' vdim 3) ('EGA' vmod 'TRID')) ;
  421.       idim = 3 ;
  422.       iaxi = FAUX ;
  423.    'FINSI' ;
  424.    'SI' ('EGA' vdim 1) ;
  425.       idim = 1 ;
  426.       iaxi = FAUX ;
  427.    'FINSI' ;
  428. *   'MESSAGE' ('CHAINE' 'iaxi=' iaxi ); 
  429.    'SI' ('EGA' idim 0) ;
  430.       'ERREUR' ('CHAINE' 'vmod=' vmod ' et vdim=' vdim ' non prevu') ;
  431.    'FINSI' ;
  432.    'SI' iaxi ;
  433.       rmt   = 'COORDONNEE' 1 _mt ;
  434.       deupi = '*' PI 2.D0 ;
  435.    'FINSI' ;
  436. *
  437. * Optimisation possible : construire la matrice par blocs
  438. * qd valt et valq ne sont pas donnés
  439. *
  440. *
  441. *Bug ?   numop  = ('**' ninc 2) '+' 1 ;
  442.    numop  = '**' ninc 2 ;
  443.    'SI' iaxi ; numop = '+' numop 1 ; 'FINSI' ;
  444.    numder = idim ;
  445.    numvar = ninc ;
  446.    ncof = '+' ncof 1 ;
  447. *delete   'SI' iaxi ; ncof = '+' ncof 1 ; 'FINSI' ;
  448.    numdat = ncof ;
  449.    numcof = ncof ;
  450. *
  451.    A = ININLIN numop numvar numdat numcof numder ;
  452. *
  453.    lvt = 'EGA' ('TYPE' valt) 'LISTREEL' ;
  454.    'REPETER' iiinc ninc ;
  455.       iinc = &iiinc ;
  456.       A . 'VAR' . iinc . 'NOMDDL' = tnomt . 'NOMPRI' . iinc ;
  457.       A . 'VAR' . iinc . 'DISC'   = tnomt . 'DISC'  ;
  458.       'SI' lvalt ;
  459.          'SI' lvt ;
  460.             A . 'VAR' . iinc . 'VALEUR'   = 'EXTRAIRE' valt iinc ;
  461.          'SINON' ;
  462.             A . 'VAR' . iinc . 'VALEUR'   = valt ;
  463.          'FINSI' ;            
  464.       'FINSI' ;         
  465.    'FIN' iiinc ;
  466. *
  467.    icof = 0 ;
  468.    ll = 'LECT' ;
  469.    icof = '+' icof 1 ;
  470.    A . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
  471.    A . 'DAT' . icof . 'DISC'   = 'CSTE' ;
  472.    A . 'DAT' . icof . 'VALEUR' = 2.D0 ;
  473.    A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
  474.    A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
  475.    ll  = 'LECT' ;
  476.    ll2 = 'LECT' icof ;
  477.    'SI' lcof ;
  478.       icof = '+' icof 1 ;
  479.       A . 'DAT' . icof . 'NOMDDL' = tcof . 'NOMPRI' . 1 ;
  480.       A . 'DAT' . icof . 'DISC'   = tcof . 'DISC' ;
  481.       A . 'DAT' . icof . 'VALEUR' = valo ;
  482.       A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
  483.       A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
  484.       ll  = 'ET' ll  ('LECT' icof) ; 
  485.       ll2 = 'ET' ll2 ('LECT' icof) ; 
  486.    'FINSI' ;
  487. *
  488.    iop = 0 ;
  489.    'REPETER' iidim idim ;
  490.       'REPETER' jidim idim ;
  491.          iop = '+' iop 1 ;
  492.          'SI' ('EGA' &iidim &jidim) ;
  493.             'SI' ('EGA' llapn 0) ;
  494.                A . iop . &iidim . &jidim = ll2 ;
  495.             'SINON' ;
  496.                A . iop . &iidim . &jidim = ll ;
  497.             'FINSI' ;
  498.          'SINON' ;
  499.             'SI' ('NEG' llapn 2) ;
  500.                A . iop . &iidim . &jidim = ll ;
  501.             'FINSI' ;
  502.             'SI' ('NEG' llapn 1) ;
  503.                A . iop . &jidim . &iidim = ll ;
  504.             'FINSI' ;               
  505.          'FINSI' ;            
  506.       'FIN' jidim ;
  507.    'FIN' iidim ;
  508.    'SI' iaxi ;   
  509.       iop = '+' iop 1 ;
  510.       'SI' ('EGA' llapn 0) ;
  511.          A . iop . 1 . 0 = ll2 ;
  512.       'SINON' ;         
  513.          A . iop . 1 . 0 = ll ;
  514.       'FINSI' ;         
  515.    'FINSI' ;
  516. *
  517.    'SI' iaxi ;
  518.       numdat = 2 ;
  519.       numcof = 2 ;
  520.    'SINON' ;      
  521.       numdat = 0 ;
  522.       numcof = 0 ;
  523.    'FINSI' ;      
  524. *
  525.    B = ININLIN numop numvar numdat numcof numder ;
  526. *
  527.    icof = 0 ;
  528. *   
  529.    'SI' iaxi ;
  530.       icof = '+' icof 1 ;
  531.       B . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
  532.       B . 'DAT' . icof . 'DISC'   = discg ;
  533.       B . 'DAT' . icof . 'VALEUR' = rmt ;
  534.       B . 'COF' . icof . 'COMPOR' = 'IDEN' ;
  535.       B . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
  536.       icof = '+' icof 1 ;
  537.       B . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
  538.       B . 'DAT' . icof . 'DISC'   = 'CSTE' ;
  539.       B . 'DAT' . icof . 'VALEUR' = deupi ;
  540.       B . 'COF' . icof . 'COMPOR' = 'IDEN' ;
  541.       B . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
  542.       ll  = 'LECT' ('-' icof 1) icof ;
  543.       mic = '*' ('-' icof 1) -1 ;
  544.       llb = 'LECT' mic icof ;
  545.    'FINSI' ;      
  546. *
  547.    lvq = 'EGA' ('TYPE' valq) 'LISTREEL' ;
  548.    'REPETER' iiinc ninc ;
  549.       iinc = &iiinc ;
  550.       B . 'VAR' . iinc . 'NOMDDL' = tnomq . 'NOMDUA' . iinc ;
  551.       B . 'VAR' . iinc . 'DISC'   = tnomq . 'DISC'  ;
  552.       'SI' lvalq ;
  553.          'SI' lvq ;
  554.             B . 'VAR' . iinc . 'VALEUR'   = 'EXTRAIRE' valq iinc ;
  555.          'SINON' ;
  556.             B . 'VAR' . iinc . 'VALEUR'   = valq ;
  557.          'FINSI' ;            
  558.       'FINSI' ;         
  559.    'FIN' iiinc ;
  560. *
  561.    iop = 0 ;
  562.    'REPETER' iidim idim ;
  563.       'REPETER' jidim idim ;
  564.          iop = '+' iop 1 ;
  565.          B . iop . &iidim . &jidim = ll ;
  566.       'FIN' jidim ;
  567.    'FIN' iidim ;
  568.    'SI' iaxi ;   
  569.       iop = '+' iop 1 ;
  570.       B . iop . 1 . 0 = llb ;
  571. *      B . iop . 1 . 0 = ll ;
  572.    'FINSI' ;      
  573. *
  574.    mgrig = NLINP discg _mt A B methgau ;
  575. * Integration par parties
  576. *   mgrig = '*' mgrig -1.D0 ;
  577. *
  578. 'RESPRO' mgrig ;
  579. 'FINPROC' ;
  580. *
  581. * End of procedure file GRIG2
  582. *
  583. *ENDPROCEDUR grig2
  584. *BEGINPROCEDUR gdiv2
  585. ************************************************************************
  586. * NOM         : GDIV2
  587. * DESCRIPTION : Une matrice de masse
  588. *
  589. *
  590. *
  591. * LANGAGE     : GIBIANE-CAST3M
  592. * AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  593. *               mél : gounand@semt2.smts.cea.fr
  594. **********************************************************************
  595. * VERSION    : v2, 14/03/2006, mise à jour NLIN évolué
  596. * VERSION    : v1, 13/05/2004, version initiale
  597. * HISTORIQUE : v1, 13/05/2004, création
  598. * HISTORIQUE :
  599. * HISTORIQUE :
  600. ************************************************************************
  601. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  602. * en cas de modification de ce sous-programme afin de faciliter
  603. * la maintenance !
  604. ************************************************************************
  605. *
  606. *
  607. 'DEBPROC' GDIV2 ;
  608. 'ARGUMENT'     _mt*'MAILLAGE' ;
  609. 'ARGUMENT'     _smt/'MAILLAGE' ;
  610. 'ARGUMENT' tdisc*'TABLE' ;
  611. *
  612. * Lectures
  613. *
  614. debug = FAUX ;
  615. lmotcle = 'MOTS' 'NPRI' 'FPRI' 'CPRI' 'NDUA' 'FDUA' 'CDUA'
  616.                  'NCOF' 'FCOF' 'CCOF' 'GBBT' 'GMBT' ;
  617. * Il faut initialiser valt et valq, sinon on peut capturer ceux de
  618. * la procédure appelante
  619. valt = 'valt' ; valq = 'valq' ;
  620. lbbt = 0 ;
  621. *           
  622. 'REPETER' imotcle ;
  623.    'ARGUMENT' motcle/'MOT' ;
  624.    'SI' ('NON' ('EXISTE' motcle)) ; 'QUITTER' imotcle ; 'FINSI' ;
  625.    'SI' ('NON' ('EXISTE' lmotcle motcle)) ;
  626.       cherr = 'CHAINE' 'Keyword ' motcle ' unknown.' ; 'ERREUR' cherr ;
  627.    'FINSI' ;
  628.    'SI' ('EGA' motcle 'NPRI') ; 'ARGUMENT' nomt*'MOT' ; 'FINSI' ;
  629.    'SI' ('EGA' motcle 'NDUA') ; 'ARGUMENT' nomq*'MOT' ; 'FINSI' ;
  630.    'SI' ('EGA' motcle 'NCOF') ; 'ARGUMENT' nomo*'MOT' ; 'FINSI' ;
  631.    'SI' ('EGA' motcle 'FPRI') ; 'ARGUMENT' valt*'LISTREEL' ; 'FINSI' ;
  632.    'SI' ('EGA' motcle 'FDUA') ; 'ARGUMENT' valq*'FLOTTANT' ; 'FINSI' ;
  633.    'SI' ('EGA' motcle 'FCOF') ; 'ARGUMENT' valo*'FLOTTANT' ; 'FINSI' ;
  634.    'SI' ('EGA' motcle 'CPRI') ; 'ARGUMENT' valt*'CHPOINT' ; 'FINSI' ;
  635.    'SI' ('EGA' motcle 'CDUA') ; 'ARGUMENT' valq*'CHPOINT' ; 'FINSI' ;
  636.    'SI' ('EGA' motcle 'CCOF') ; 'ARGUMENT' valo*'CHPOINT' ; 'FINSI' ;
  637.    'SI' ('EGA' motcle 'GBBT') ; lbbt = 1 ; 'FINSI' ;
  638.    'SI' ('EGA' motcle 'GMBT') ; lbbt = 2 ; 'FINSI' ;
  639. 'FIN' imotcle ;
  640. *
  641. * Tests
  642. *
  643. discg = TDISC . 'GEOM' . 'DISC' ;
  644. 'SI' ('EXISTE' tdisc 'methgau') ;
  645.    methgau = tdisc . 'methgau' . 'amor' ;
  646. 'SINON' ;
  647.    methgau = 'GAU7' ;
  648. 'FINSI' ;
  649. tnomt = TDISC . nomt ;
  650. lvalt = 'NEG' ('TYPE' valt) 'MOT' ;
  651. tnomq = TDISC . nomq ;
  652. lvalq = 'NEG' ('TYPE' valq) 'MOT' ;
  653. *
  654. lcof = 'EXISTE' TDISC nomo ;
  655. 'SI' lcof ; ncof = 1 ; tcof = TDISC . nomo ;
  656. 'SINON' ; ncof = 0 ;
  657. 'FINSI' ;
  658. *
  659. 'SI' debug ;
  660. 'SI' lcof ; 'MESSAGE' 'Un coef a ete detecte' ;
  661. 'SINON' ; 'MESSAGE' 'pas de coef  detecte' ;
  662. 'FINSI' ;
  663. 'FINSI' ;
  664. *
  665.    vdim = 'VALEUR' 'DIME' ;
  666.    vmod = 'VALEUR' 'MODE' ;
  667.    idim = 0 ;
  668.    'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'PLANDEFO')) ;
  669.       idim = 2 ;
  670.       iaxi = FAUX ;
  671.    'FINSI' ;
  672.    'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'AXIS')) ;
  673.       idim = 2 ;
  674.       iaxi = VRAI ;
  675.    'FINSI' ;
  676.    'SI' ('ET' ('EGA' vdim 3) ('EGA' vmod 'TRID')) ;
  677.       idim = 3 ;
  678.       iaxi = FAUX ;
  679.    'FINSI' ;
  680.    'SI' ('EGA' vdim 1) ;
  681.       idim = 1 ;
  682.       iaxi = FAUX ;
  683.    'FINSI' ;
  684.    'SI' ('EGA' idim 0) ;
  685.       'ERREUR' ('CHAINE' 'vmod=' vmod ' et vdim=' vdim ' non prevu') ;
  686.    'FINSI' ;
  687.    'SI' iaxi ;
  688.       dp = ('*' PI 2.D0) ;
  689.       rmt = 'COORDONNEE' 1 _mt ;
  690.       ncof = ncof '+' 2 ;
  691.    'FINSI' ;
  692. * Scalaire ou vecteur
  693. ninct = 'DIME' (tnomt . 'NOMPRI') ;
  694. nincq = 'DIME' (tnomq . 'NOMPRI') ;
  695. 'SI' ('NEG' ninct idim) ;
  696.    cherr = 'CHAINE'
  697.    'la primale doit etre un vecteur' ;
  698.    'ERREUR' cherr ;
  699. 'FINSI' ;   
  700. 'SI' ('NEG' nincq 1) ;
  701.    cherr = 'CHAINE'
  702.    'la duale doit etre un scalaire' ;
  703.    'ERREUR' cherr ;
  704. 'FINSI' ;   
  705. *
  706.    numop = 1 ; numder = idim ; numvar = ninct ;
  707.    numdat = ncof ; numcof = ncof ;
  708.    A = ININLIN numop numvar numdat numcof numder ;
  709. *   
  710.    lvt = 'EGA' ('TYPE' valt) 'LISTREEL' ;
  711.    'REPETER' iiinct ninct ;
  712.       iinct = &iiinct ;
  713.       A . 'VAR' . iinct . 'NOMDDL' = tnomt . 'NOMPRI' . iinct ;
  714.       A . 'VAR' . iinct . 'DISC'   = tnomt . 'DISC'  ;
  715.       'SI' lvalt ;
  716.          'SI' lvt ;
  717.             A . 'VAR' . iinct . 'VALEUR'   = 'EXTRAIRE' valt iinct ;
  718.          'SINON' ;
  719.             A . 'VAR' . iinct . 'VALEUR'   = valt ;
  720.          'FINSI' ;            
  721.       'FINSI' ;         
  722.    'FIN' iiinct ;
  723. *
  724.    icof = 0 ;
  725.    'SI' lcof ;
  726.       icof = '+' icof 1 ;
  727.       A . 'DAT' . icof . 'NOMDDL' = tcof . 'NOMPRI' . 1 ;
  728.       A . 'DAT' . icof . 'DISC'   = tcof . 'DISC' ;
  729.       A . 'DAT' . icof . 'VALEUR' = valo ;
  730.       A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
  731.       A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
  732.       ll = 'LECT' 1 ;
  733.    'SINON' ;
  734.       ll = 'LECT' ;
  735.    'FINSI' ;
  736. *   
  737.    'SI' iaxi ;
  738.       icof = '+' icof 1 ;
  739.       A . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
  740.       A . 'DAT' . icof . 'DISC'   = 'CSTE' ;
  741.       A . 'DAT' . icof . 'VALEUR' = dp ;
  742.       A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
  743.       A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
  744.       icof = '+' icof 1 ;
  745.       A . 'DAT' . icof . 'NOMDDL' = 'MOTS' 'SCAL' ;
  746.       A . 'DAT' . icof . 'DISC'   = discg ;
  747.       A . 'DAT' . icof . 'VALEUR' = rmt ;
  748.       A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
  749.       A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
  750.       lldpr = ll 'ET' ('LECT' ('-' icof 1) icof) ;
  751.       lldp = ll 'ET' ('LECT' ('-' icof 1)) ;
  752.    'FINSI' ;
  753. *
  754.    'SI' iaxi ;
  755.       'REPETER' iidim idim ;
  756.          A . 1 . &iidim . &iidim  = lldpr ;
  757.       'FIN' iidim ;
  758.       A . 1 . 1 . 0 = lldp ;
  759.    'SINON' ;
  760.       'REPETER' iidim idim ;
  761.          A . 1 . &iidim . &iidim  = ll ;
  762.       'FIN' iidim ;
  763.    'FINSI' ;
  764. *
  765.    numvar = 1 ;
  766.    numdat = 0 ;
  767.    numcof = 0 ;
  768. *
  769.    B = ININLIN numop numvar numdat numcof numder ;
  770.    B . 'VAR' . 1 . 'NOMDDL' = tnomq . 'NOMDUA' . 1 ;
  771.    B . 'VAR' . 1 . 'DISC'   = tnomq . 'DISC'  ;
  772.    'SI' lvalq ;
  773.       B . 'VAR' . 1 . 'VALEUR'   = valq ;
  774.    'FINSI' ;         
  775.    B . 1 . 1 . 0  = 'LECT' ;
  776. *
  777.    'SI' ('OU' ('EGA' lbbt 0) ('EGA' lbbt 1)) ;
  778.       'SI' ('EXISTE' _smt) ;
  779.           mgdiv2 = 'NLIN' discg _mt _smt A B methgau ;
  780.       'SINON' ;       
  781.           mgdiv2 = NLINP discg _mt A B methgau ;
  782.       'FINSI' ;
  783.    'FINSI' ;      
  784.    'SI' ('OU' ('EGA' lbbt 1) ('EGA' lbbt 2)) ;
  785.       B . 'VAR' . 1 . 'NOMDDL' = tnomq . 'NOMPRI' . 1 ;
  786.       'REPETER' iiinct ninct ;
  787.          iinct = &iiinct ;
  788.          A . 'VAR' . iinct . 'NOMDDL' = tnomt . 'NOMDUA' . iinct ;
  789.       'FIN' iiinct ;         
  790.       'SI' ('EXISTE' _smt) ;
  791.           mgdiv3 = 'NLIN' discg _mt _smt B A methgau ;
  792.       'SINON' ;       
  793.           mgdiv3 = NLINP discg _mt B A methgau ;
  794.       'FINSI' ;
  795.    'FINSI' ;
  796.    'SI' ('EGA' lbbt 0) ;
  797.       mgdiv = mgdiv2 ;
  798.    'FINSI' ;      
  799.    'SI' ('EGA' lbbt 1) ;
  800.       mgdiv = mgdiv2 'ET' mgdiv3 ;
  801.    'FINSI' ;      
  802.    'SI' ('EGA' lbbt 2) ;
  803.       mgdiv = mgdiv3 ;
  804.    'FINSI' ;      
  805. 'RESPRO' mgdiv ;
  806. 'FINPROC' ;
  807. *
  808. * End of procedure file GDIV2
  809. *
  810. *ENDPROCEDUR gdiv2
  811. *BEGINPROCEDUR gnor2
  812. ************************************************************************
  813. * NOM         : GNOR2
  814. * DESCRIPTION : Calcule le champ de normales à une surface.
  815. *               Peut servir à calculer une pression, un potentiel
  816. *               lié à la gravité, un volume contenu dans une surface.
  817. *               Attention à l'orientation de la surface !
  818. *
  819. *               Computes a field of normal to a surface.
  820. *               Also useful to compute a pressure field,
  821. *               a gravity potential field, a volume enclosed
  822. *               by a surface.
  823. *               WARNING : The orientation of the surface matters !
  824. *
  825. *
  826. *
  827. *
  828. * LANGAGE     : GIBIANE-CAST3M
  829. * AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  830. *               mél : gounand@semt2.smts.cea.fr
  831. **********************************************************************
  832. * VERSION    : v1, 22/04/2011
  833. * HISTORIQUE : v1, 22/04/2011, création
  834. * HISTORIQUE :
  835. * HISTORIQUE :
  836. ************************************************************************
  837. *
  838. *
  839. 'DEBPROC' GNOR2 ;
  840. 'ARGUMENT'     _mt*'MAILLAGE' ;
  841. 'ARGUMENT' tdisc*'TABLE' ;
  842. *
  843. * Lectures
  844. *
  845. dim = 'VALEUR' 'DIME' ;
  846. debug = FAUX ;
  847. lmotcle = 'MOTS' 'NPRI' 'FPRI' 'CPRI' 'NDUA' 'FDUA' 'CDUA'
  848.                  'NCOF' 'FCOF' 'CCOF' ;
  849. * Il faut initialiser valt et valq, sinon on peut capturer ceux de
  850. * la procédure appelante
  851. valt = 'valt' ; valq = 'valq' ;
  852. 'REPETER' imotcle ;
  853.    'ARGUMENT' motcle/'MOT' ;
  854.    'SI' ('NON' ('EXISTE' motcle)) ; 'QUITTER' imotcle ; 'FINSI' ;
  855.    'SI' ('NON' ('EXISTE' lmotcle motcle)) ;
  856.       cherr = 'CHAINE' 'Keyword ' motcle ' unknown.' ; 'ERREUR' cherr ;
  857.    'FINSI' ;
  858.    'SI' ('EGA' motcle 'NPRI') ; 'ARGUMENT' nomt*'MOT' ; 'FINSI' ;
  859.    'SI' ('EGA' motcle 'NDUA') ; 'ARGUMENT' nomq*'MOT' ; 'FINSI' ;
  860.    'SI' ('EGA' motcle 'NCOF') ; 'ARGUMENT' nomo*'MOT' ; 'FINSI' ;
  861.    tst1 = 'EGA' motcle 'FPRI' ; tst2 = 'EGA' motcle 'FDUA' ;
  862.    tst3 = 'EGA' motcle 'FCOF' ;
  863.    tst = tst1 'OU' tst2 'OU' tst3 ;
  864.    'SI' tst ;
  865.       'SI' tst1 ; tt = TDISC . nomt ; 'FINSI' ;
  866.       'SI' tst2 ; tt = TDISC . nomq ; 'FINSI' ;
  867.       'SI' tst3 ; tt = TDISC . nomo ; 'FINSI' ;
  868.       isvec = ('>' ('DIME' (tt. 'NOMPRI')) 1) ;
  869.       'SI' isvec ; 'ARGUMENT' valv*'LISTREEL' ; 'SINON' ;
  870.          'ARGUMENT' valv*'FLOTTANT' ;
  871.       'FINSI' ;
  872.       'SI' tst1 ; valt = valv ; 'FINSI' ;
  873.       'SI' tst2 ; valq = valv ; 'FINSI' ;
  874.       'SI' tst3 ; valo = valv ; 'FINSI' ;
  875.    'FINSI' ;
  876.    'SI' ('EGA' motcle 'CPRI') ; 'ARGUMENT' valt*'CHPOINT' ; 'FINSI' ;
  877.    'SI' ('EGA' motcle 'CDUA') ; 'ARGUMENT' valq*'CHPOINT' ; 'FINSI' ;
  878.    'SI' ('EGA' motcle 'CCOF') ; 'ARGUMENT' valo*'CHPOINT' ; 'FINSI' ;
  879. 'FIN' imotcle ;
  880. *
  881. * Tests
  882. *
  883. discg = TDISC . 'GEOM' . 'DISC' ;
  884. 'SI' ('EXISTE' tdisc 'methgau') ;
  885.    methgau = tdisc . 'methgau' . 'mass' ;
  886. 'SINON' ;
  887.    methgau = 'GAU7' ;
  888. 'FINSI' ;
  889. tnomt = TDISC . nomt ;
  890. lvalt = 'NEG' ('TYPE' valt) 'MOT' ;
  891. tnomq = TDISC . nomq ;
  892. lvalq = 'NEG' ('TYPE' valq) 'MOT' ;
  893. * Scalaire ou vecteur
  894. ninct = 'DIME' (tnomt . 'NOMPRI') ;
  895. nincq = 'DIME' (tnomq . 'NOMDUA') ;
  896. 'SI' ('ET' ('NEG' ninct 1) ('NEG' ninct dim)) ;
  897.    cherr = 'CHAINE'
  898.    'la primale doit etre un scalaire ou un vecteur' ;
  899.    'ERREUR' cherr ;
  900. 'FINSI' ;   
  901. 'SI' ('NEG' nincq dim) ;
  902.    cherr = 'CHAINE'
  903.    'la duale doit etre un vecteur' ;
  904.    'ERREUR' cherr ;
  905. 'FINSI' ;   
  906. *ninc = ninct ;
  907. *
  908. lcof = 'EXISTE' TDISC nomo ;
  909. 'SI' lcof ; tcof = TDISC . nomo ;
  910. ncof = 'DIME' (tcof . 'NOMPRI') ;
  911. 'SINON' ; ncof = 0 ;
  912. 'FINSI' ;
  913. *
  914. 'SI' debug ;
  915. 'SI' lcof ; 'MESSAGE' 'Un coef a ete detecte' ;
  916. 'SINON' ; 'MESSAGE' 'pas de coef  detecte' ;
  917. 'FINSI' ;
  918. 'FINSI' ;
  919. *
  920.    vdim = 'VALEUR' 'DIME' ;
  921.    vmod = 'VALEUR' 'MODE' ;
  922.    idim = 0 ;
  923.    'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'PLANDEFO')) ;
  924.       idim = 2 ;
  925.       iaxi = FAUX ;
  926.    'FINSI' ;
  927.    'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'AXIS')) ;
  928.       idim = 2 ;
  929.       iaxi = VRAI ;
  930.    'FINSI' ;
  931.    'SI' ('ET' ('EGA' vdim 3) ('EGA' vmod 'TRID')) ;
  932.       idim = 3 ;
  933.       iaxi = FAUX ;
  934.    'FINSI' ;
  935.    'SI' ('EGA' vdim 1) ;
  936.       idim = 1 ;
  937.       iaxi = FAUX ;
  938.    'FINSI' ;
  939. *   'MESSAGE' ('CHAINE' 'iaxi=' iaxi ); 
  940.    'SI' ('EGA' idim 0) ;
  941.       'ERREUR' ('CHAINE' 'vmod=' vmod ' et vdim=' vdim ' non prevu') ;
  942.    'FINSI' ;
  943.    'SI' iaxi ;
  944.       dprmt = '*' ('COORDONNEE' 1 _mt) ('*' PI 2.D0) ;
  945.    'FINSI' ;
  946. *
  947. * Optimisation possible : construire la matrice par blocs
  948. * qd valt et valq ne sont pas donnés
  949. *
  950.    numop = nincq ; numder = idim ; numvar = ninct ;
  951.    numdat = ncof ; numcof = ncof ;
  952.    A = ININLIN numop numvar numdat numcof numder ;
  953.    'SI' lcof ;
  954.       lvo = 'EGA' ('TYPE' valo) 'LISTREEL' ;
  955.       'REPETER' iicof ncof ;
  956.          icof = &iicof ;
  957.          A . 'DAT' . icof . 'NOMDDL' = tcof . 'NOMPRI' . icof ;
  958.          A . 'DAT' . icof . 'DISC'   = tcof . 'DISC' ;
  959.          'SI' lvo ;
  960.             A . 'DAT' . icof . 'VALEUR' = 'EXTRAIRE' valo icof ;
  961.          'SINON' ;            
  962.             A . 'DAT' . icof . 'VALEUR' = valo ;
  963.          'FINSI' ;            
  964.          A . 'COF' . icof . 'COMPOR' = 'IDEN' ;
  965.          A . 'COF' . icof . 'LDAT'   = 'LECT' icof ;
  966.       'FIN' iicof ;         
  967.    'SINON' ;
  968.       ll = 'LECT' ;
  969.    'FINSI' ;
  970.    lvt = 'EGA' ('TYPE' valt) 'LISTREEL' ;
  971.    'REPETER' iiincq nincq ;
  972.       iincq = &iiincq ;
  973.       iinct = 'MINIMUM' ('LECT' iincq ninct) ;
  974.       A . 'VAR' . iinct . 'NOMDDL' = tnomt . 'NOMPRI' . iinct ;
  975.       A . 'VAR' . iinct . 'DISC'   = tnomt . 'DISC'  ;
  976.       'SI' lvalt ;
  977.          'SI' lvt ;
  978.             A . 'VAR' . iinct . 'VALEUR'   = 'EXTRAIRE' valt iinct ;
  979.          'SINON' ;
  980.             A . 'VAR' . iinct . 'VALEUR'   = valt ;
  981.          'FINSI' ;            
  982.       'FINSI' ;
  983.       'SI' lcof ;
  984.          icof = 'MINIMUM' ('LECT' iincq ncof) ;
  985.          A . iincq . iinct . 0  = 'LECT' icof ;
  986.       'SINON' ;         
  987.          A . iincq . iinct . 0  = ll ;
  988.       'FINSI' ;         
  989.    'FIN' iiincq ;
  990. *
  991.    'SI' iaxi ;
  992.       numdat = 1 ;
  993.       numcof = dim '+' 1 ;
  994.    'SINON' ;      
  995.       numdat = 0 ;
  996.       numcof = dim ;
  997.    'FINSI' ;      
  998.    numvar = nincq ;
  999.    B = ININLIN numop numvar numdat numcof numder ;
  1000.    icof = 0 ;
  1001.    'REPETER' iiidim idim ;
  1002.       icof = '+' icof 1 ;
  1003.       B . 'COF' . icof . 'COMPOR' = 'CHAINE' 'VNOR' &iiidim ;
  1004.       B . 'COF' . icof . 'LDAT'   = 'LECT' ;
  1005.    'FIN' iiidim ;
  1006. *
  1007.    'SI' iaxi ;
  1008.       icof = '+' icof 1 ;
  1009.       B . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  1010.       B . 'DAT' . 1 . 'DISC'   = discg ;
  1011.       B . 'DAT' . 1 . 'VALEUR' = dprmt ;
  1012.       B . 'COF' . icof . 'COMPOR' = 'IDEN' ;
  1013.       B . 'COF' . icof . 'LDAT'   = 'LECT' 1 ;
  1014.       ll = 'LECT' icof ;
  1015.    'SINON' ;
  1016.       ll = 'LECT' ;
  1017.    'FINSI' ;
  1018.    lvq = 'EGA' ('TYPE' valq) 'LISTREEL' ;
  1019.    'REPETER' iiincq nincq ;
  1020.       iincq = &iiincq ;
  1021.       B . 'VAR' . iincq . 'NOMDDL' = tnomq . 'NOMDUA' . iincq ;
  1022.       B . 'VAR' . iincq . 'DISC'   = tnomq . 'DISC'  ;
  1023.       'SI' lvalq ;
  1024.          'SI' lvq ;
  1025.             B . 'VAR' . iincq . 'VALEUR'   = 'EXTRAIRE' valq iincq ;
  1026.          'SINON' ;
  1027.             B . 'VAR' . iincq . 'VALEUR'   = valq ;
  1028.          'FINSI' ;            
  1029.       'FINSI' ;         
  1030.       B . iincq . iincq . 0  = ('LECT' iincq) 'ET' ll ;
  1031.    'FIN' iiincq ;
  1032. *
  1033.    mgnor2 = NLIN discg _mt A B methgau ;
  1034. *
  1035. 'RESPRO' mgnor2 ;
  1036. 'FINPROC' ;
  1037. *
  1038. * End of procedure file GNOR2
  1039. *
  1040. *ENDPROCEDUR gnor2
  1041. *BEGINPROCEDUR gmass2
  1042. ************************************************************************
  1043. * NOM         : GMASS2
  1044. * DESCRIPTION : Une matrice de masse
  1045. *
  1046. *
  1047. *
  1048. * LANGAGE     : GIBIANE-CAST3M
  1049. * AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  1050. *               mél : gounand@semt2.smts.cea.fr
  1051. **********************************************************************
  1052. * VERSION    : v2, 14/03/2006, mise à jour NLIN évolué
  1053. * VERSION    : v1, 13/05/2004, version initiale
  1054. * HISTORIQUE : v1, 13/05/2004, création
  1055. * HISTORIQUE :
  1056. * HISTORIQUE :
  1057. ************************************************************************
  1058. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  1059. * en cas de modification de ce sous-programme afin de faciliter
  1060. * la maintenance !
  1061. ************************************************************************
  1062. *
  1063. *
  1064. 'DEBPROC' GMASS2 ;
  1065. 'ARGUMENT'     _mt*'MAILLAGE' ;
  1066. 'ARGUMENT'     _smt/'MAILLAGE' ;
  1067. 'ARGUMENT' tdisc*'TABLE' ;
  1068. *
  1069. * Lectures
  1070. *
  1071. debug = FAUX ;
  1072. lmotcle = 'MOTS' 'NPRI' 'FPRI' 'CPRI' 'NDUA' 'FDUA' 'CDUA'
  1073.                  'NCOF' 'FCOF' 'CCOF' ;
  1074. * Il faut initialiser valt et valq, sinon on peut capturer ceux de
  1075. * la procédure appelante
  1076. valt = 'valt' ; valq = 'valq' ;
  1077. 'REPETER' imotcle ;
  1078.    'ARGUMENT' motcle/'MOT' ;
  1079.    'SI' ('NON' ('EXISTE' motcle)) ; 'QUITTER' imotcle ; 'FINSI' ;
  1080.    'SI' ('NON' ('EXISTE' lmotcle motcle)) ;
  1081.       cherr = 'CHAINE' 'Keyword ' motcle ' unknown.' ; 'ERREUR' cherr ;
  1082.    'FINSI' ;
  1083.    'SI' ('EGA' motcle 'NPRI') ; 'ARGUMENT' nomt*'MOT' ; 'FINSI' ;
  1084.    'SI' ('EGA' motcle 'NDUA') ; 'ARGUMENT' nomq*'MOT' ; 'FINSI' ;
  1085.    'SI' ('EGA' motcle 'NCOF') ; 'ARGUMENT' nomo*'MOT' ; 'FINSI' ;
  1086.    tst1 = 'EGA' motcle 'FPRI' ; tst2 = 'EGA' motcle 'FDUA' ;
  1087.    tst = tst1 'OU' tst2 ;
  1088.    'SI' tst ;
  1089.       'SI' tst1 ; tt = TDISC . nomt ; 'FINSI' ;
  1090.       'SI' tst2 ; tt = TDISC . nomq ; 'FINSI' ;
  1091.       isvec = ('>' ('DIME' (tt. 'NOMPRI')) 1) ;
  1092.       'SI' isvec ; 'ARGUMENT' valv*'LISTREEL' ; 'SINON' ;
  1093.          'ARGUMENT' valv*'FLOTTANT' ;
  1094.       'FINSI' ;
  1095.       'SI' tst1 ; valt = valv ; 'FINSI' ;
  1096.       'SI' tst2 ; valq = valv ; 'FINSI' ;
  1097.    'FINSI' ;
  1098.    'SI' ('EGA' motcle 'FCOF') ; 'ARGUMENT' valo*'FLOTTANT' ; 'FINSI' ;
  1099.    'SI' ('EGA' motcle 'CPRI') ; 'ARGUMENT' valt*'CHPOINT' ; 'FINSI' ;
  1100.    'SI' ('EGA' motcle 'CDUA') ; 'ARGUMENT' valq*'CHPOINT' ; 'FINSI' ;
  1101.    'SI' ('EGA' motcle 'CCOF') ; 'ARGUMENT' valo*'CHPOINT' ; 'FINSI' ;
  1102. 'FIN' imotcle ;
  1103. *
  1104. * Tests
  1105. *
  1106. discg = TDISC . 'GEOM' . 'DISC' ;
  1107. 'SI' ('EXISTE' tdisc 'methgau') ;
  1108.    methgau = tdisc . 'methgau' . 'mass' ;
  1109. 'SINON' ;
  1110.    methgau = 'GAU7' ;
  1111. 'FINSI' ;
  1112. tnomt = TDISC . nomt ;
  1113. lvalt = 'NEG' ('TYPE' valt) 'MOT' ;
  1114. tnomq = TDISC . nomq ;
  1115. lvalq = 'NEG' ('TYPE' valq) 'MOT' ;
  1116. * Scalaire ou vecteur
  1117. ninct = 'DIME' (tnomt . 'NOMPRI') ;
  1118. nincq = 'DIME' (tnomq . 'NOMPRI') ;
  1119. 'SI' ('NEG' ninct nincq) ;
  1120.    cherr = 'CHAINE'
  1121.    'les primales et duales nont pas le meme nombre de composantes' ;
  1122.    'ERREUR' cherr ;
  1123. 'FINSI' ;   
  1124. ninc = ninct ;
  1125. *
  1126. lcof = 'EXISTE' TDISC nomo ;
  1127. 'SI' lcof ; ncof = 1 ; tcof = TDISC . nomo ;
  1128. 'SINON' ; ncof = 0 ;
  1129. 'FINSI' ;
  1130. *
  1131. 'SI' debug ;
  1132. 'SI' lcof ; 'MESSAGE' 'Un coef a ete detecte' ;
  1133. 'SINON' ; 'MESSAGE' 'pas de coef  detecte' ;
  1134. 'FINSI' ;
  1135. 'FINSI' ;
  1136. *
  1137.    vdim = 'VALEUR' 'DIME' ;
  1138.    vmod = 'VALEUR' 'MODE' ;
  1139.    idim = 0 ;
  1140.    'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'PLANDEFO')) ;
  1141.       idim = 2 ;
  1142.       iaxi = FAUX ;
  1143.    'FINSI' ;
  1144.    'SI' ('ET' ('EGA' vdim 2) ('EGA' vmod 'AXIS')) ;
  1145.       idim = 2 ;
  1146.       iaxi = VRAI ;
  1147.    'FINSI' ;
  1148.    'SI' ('ET' ('EGA' vdim 3) ('EGA' vmod 'TRID')) ;
  1149.       idim = 3 ;
  1150.       iaxi = FAUX ;
  1151.    'FINSI' ;
  1152.    'SI' ('EGA' vdim 1) ;
  1153.       idim = 1 ;
  1154.       iaxi = FAUX ;
  1155.    'FINSI' ;
  1156. *   'MESSAGE' ('CHAINE' 'iaxi=' iaxi ); 
  1157.    'SI' ('EGA' idim 0) ;
  1158.       'ERREUR' ('CHAINE' 'vmod=' vmod ' et vdim=' vdim ' non prevu') ;
  1159.    'FINSI' ;
  1160.    'SI' iaxi ;
  1161.       dprmt = '*' ('COORDONNEE' 1 _mt) ('*' PI 2.D0) ;
  1162.    'FINSI' ;
  1163. *
  1164. * Optimisation possible : construire la matrice par blocs
  1165. * qd valt et valq ne sont pas donnés
  1166. *
  1167.    numop = ninc ; numder = idim ; numvar = ninc ;
  1168.    numdat = ncof ; numcof = ncof ;
  1169.    A = ININLIN numop numvar numdat numcof numder ;
  1170.    'SI' lcof ;
  1171.       A . 'DAT' . 1 . 'NOMDDL' = tcof . 'NOMPRI' . 1 ;
  1172.       A . 'DAT' . 1 . 'DISC'   = tcof . 'DISC' ;
  1173.       A . 'DAT' . 1 . 'VALEUR' = valo ;
  1174.       A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  1175.       A . 'COF' . 1 . 'LDAT'   = 'LECT' 1 ;
  1176.       ll = 'LECT' 1 ;
  1177.    'SINON' ;
  1178.       ll = 'LECT' ;
  1179.    'FINSI' ;
  1180.    lvt = 'EGA' ('TYPE' valt) 'LISTREEL' ;
  1181.    'REPETER' iiinc ninc ;
  1182.       iinc = &iiinc ;
  1183.       A . 'VAR' . iinc . 'NOMDDL' = tnomt . 'NOMPRI' . iinc ;
  1184.       A . 'VAR' . iinc . 'DISC'   = tnomt . 'DISC'  ;
  1185.       'SI' lvalt ;
  1186.          'SI' lvt ;
  1187.             A . 'VAR' . iinc . 'VALEUR'   = 'EXTRAIRE' valt iinc ;
  1188.          'SINON' ;
  1189.             A . 'VAR' . iinc . 'VALEUR'   = valt ;
  1190.          'FINSI' ;            
  1191.       'FINSI' ;         
  1192.       A . iinc . iinc . 0  = ll ;
  1193.    'FIN' iiinc ;
  1194. *
  1195.    'SI' iaxi ;
  1196.       numdat = 1 ;
  1197.       numcof = 1 ;
  1198.    'SINON' ;      
  1199.       numdat = 0 ;
  1200.       numcof = 0 ;
  1201.    'FINSI' ;      
  1202.    B = ININLIN numop numvar numdat numcof numder ;
  1203.    'SI' iaxi ;
  1204.       B . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  1205.       B . 'DAT' . 1 . 'DISC'   = discg ;
  1206.       B . 'DAT' . 1 . 'VALEUR' = dprmt ;
  1207.       B . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  1208.       B . 'COF' . 1 . 'LDAT'   = 'LECT' 1 ;
  1209.       ll = 'LECT' 1 ;
  1210.    'SINON' ;
  1211.       ll = 'LECT' ;
  1212.    'FINSI' ;
  1213.    lvq = 'EGA' ('TYPE' valq) 'LISTREEL' ;
  1214.    'REPETER' iiinc ninc ;
  1215.       iinc = &iiinc ;
  1216.       B . 'VAR' . iinc . 'NOMDDL' = tnomq . 'NOMDUA' . iinc ;
  1217.       B . 'VAR' . iinc . 'DISC'   = tnomq . 'DISC'  ;
  1218.       'SI' lvalq ;
  1219.          'SI' lvq ;
  1220.             B . 'VAR' . iinc . 'VALEUR'   = 'EXTRAIRE' valq iinc ;
  1221.          'SINON' ;
  1222.             B . 'VAR' . iinc . 'VALEUR'   = valq ;
  1223.          'FINSI' ;            
  1224.       'FINSI' ;         
  1225.       B . iinc . iinc . 0  = ll ;
  1226.    'FIN' iiinc ;
  1227. *
  1228.    'SI' ('EXISTE' _smt) ;
  1229.        mgmass2 = 'NLIN' discg _mt _smt A B methgau ;
  1230.    'SINON' ;       
  1231.        mgmass2 = NLINP discg _mt A B methgau ;
  1232.    'FINSI' ;       
  1233. *
  1234. 'RESPRO' mgmass2 ;
  1235. 'FINPROC' ;
  1236. *
  1237. * End of procedure file GMASS2
  1238. *
  1239. *ENDPROCEDUR gmass2
  1240. *
  1241. *
  1242. *************************************************************************
  1243. *
  1244. * 
  1245. *   MAIN : 1) MESH
  1246. *          2) COMPUTATIONAL LOOP
  1247. *          3) TESTs
  1248. *
  1249. ************************************************************************
  1250. *
  1251. *
  1252. * Maillage et modèle
  1253. *
  1254. *nx = 1 ; ny =1 ;
  1255. *nx = 1 ; ny =2 ;
  1256. nx = 10 ; ny= 10 ; 
  1257. pA = 0. 0. ; pB = 1. 0. ; pC = 1. 1. ; pD = 0. 1. ;
  1258. ang = 32.1 ;
  1259. pA pB pC pD = 'TOURNER' pA pB pC pD ang pA ;
  1260. discu = 'QUAF' ; discp = 'LINE' ;
  1261. 'OPTI' 'ELEM' 'QUA4' ;
  1262. tol = 1.D-6 ;
  1263. bas = 'DROIT' nx pA pB ; 
  1264. dro = 'DROIT' ny pB pC ; 
  1265. hau = 'DROIT' nx pC pD ; 
  1266. gau = 'DROIT' ny pD pA ;
  1267. mt = 'DALLER' bas dro hau gau ;
  1268. *
  1269. tmail = QUAFME mt 'mt' tol ;
  1270. tmail = QUAFME tmail bas 'bas' dro 'dro' hau 'hau' gau 'gau' tol ;
  1271. _mt = tmail . 'mt' . 'QUAF' ;
  1272. mt  = tmail . 'mt' . discu ;
  1273. mtp = tmail . 'mt' . discp ;
  1274. cmt = 'CONTOUR' mt ;
  1275. bas = tmail . 'bas' . discu ;
  1276. dro = tmail . 'dro' . discu ;
  1277. hau = tmail . 'hau' . discu ;
  1278. gau = tmail . 'gau' . discu ;
  1279. _bas = tmail . 'bas' . 'QUAF' ;
  1280. _dro = tmail . 'dro' . 'QUAF' ;
  1281. _hau = tmail . 'hau' . 'QUAF' ;
  1282. _gau = tmail . 'gau' . 'QUAF' ;
  1283. *
  1284. tmode = MODENLIN 'GEOM' 'LINE'
  1285.                  'INCO' 'UN' discu 'VECT' 'UX' 'UY' 'FX' 'FY'
  1286.                  'INCO' 'PN' discp 'SCAL' 'LXP'  'LXP'
  1287.                  ;
  1288. nompp = 'EXTRAIRE' (tmode . 'PN' .'NOMPRI' . 1) ('LECT' 1) ;                 
  1289. nompd = 'EXTRAIRE' (tmode . 'PN' .'NOMDUA' . 1) ('LECT' 1) ;                 
  1290. nomvp = 'ET' (tmode . 'UN' .'NOMPRI' . 1)
  1291.              (tmode . 'UN' .'NOMPRI' . 2) ; 
  1292. nomvd = 'ET' (tmode . 'UN' .'NOMDUA' . 1)
  1293.              (tmode . 'UN' .'NOMDUA' . 2) ; 
  1294. discg = tmode . 'GEOM' . 'DISC' ;
  1295. *
  1296. * Matrice de rigidité
  1297. *
  1298. *rig  = GRIG2 _mt tmode 'NPRI' 'UN' 'NDUA' 'UN' 'LAPN' ;
  1299. rig  = GRIG2 _mt tmode 'NPRI' 'UN' 'NDUA' 'UN' ;
  1300. *
  1301. * Matrice de la contrainte de divergence nulle et sa transposée
  1302. *
  1303. cnt1 = GDIV2 _mt tmode 'NPRI' 'UN' 'NDUA' 'PN' ;
  1304. cnt2 = GDIV2 _mt tmode 'NPRI' 'UN' 'NDUA' 'PN' 'GMBT' ;
  1305. mcnt = cnt1 'ET' cnt2 ;
  1306. *
  1307. * Conditions aux limites
  1308. *
  1309. * On relache la contrainte sur la vitesse normale en 1 point (pmil) pour
  1310. * fixer le mode de pression constante
  1311. pmil = * ('PLUS' pC pD) 0.5 ;
  1312. pmil = 'POIN' hau 'PROC' pmil ;
  1313. haur = 'DIFF' ('CHANGER' 'POI1' hau)
  1314.               (pmil 'ET' pC 'ET' pD) ;
  1315. gaur = 'DIFF' ('CHANGER' 'POI1' gau)
  1316.               ('ET' pA pD) ;
  1317. dror = 'DIFF' ('CHANGER' 'POI1' dro)
  1318.               ('ET' pB pC) ;
  1319. xbas ybas = 'COORDONNEE' bas ;
  1320. sbas = '**' ('+' ('**' xbas 2.) ('**' ybas 2)) 0.5 ;
  1321. sbas = '**' ('SIN' ('*' sbas 180.)) 2. ;
  1322. ssbas = 'SIN' ('*' 90. sbas) ;
  1323. vAB = 'MOIN' pB pA ;
  1324. xab = 'COORDONNEE' 1 vAB ;
  1325. yab = 'COORDONNEE' 2 vAB ;
  1326. *
  1327. ssbas = '+' ('NOMC' 'UX' ('*' ssbas xab))
  1328.             ('NOMC' 'UY' ('*' ssbas yab)) ;
  1329. *
  1330. * Construction de la normale
  1331. *
  1332. vnorh  = GNOR2 _hau tmode 'NPRI' discg 'FPRI' 1. 'NDUA' 'UN' ;
  1333. vnorg  = GNOR2 _gau tmode 'NPRI' discg 'FPRI' 1. 'NDUA' 'UN' ;
  1334. vnord  = GNOR2 _dro tmode 'NPRI' discg 'FPRI' 1. 'NDUA' 'UN' ;
  1335. vnorh = 'NOMC' nomvd nomvp vnorh ;
  1336. vnorg = 'NOMC' nomvd nomvp vnorg ;
  1337. vnord = 'NOMC' nomvd nomvp vnord ;
  1338. *
  1339. mblh = 'BLOQUE' 'DEPL' 'DIRE' vnorh haur ;
  1340. mblg = 'BLOQUE' 'DEPL' 'DIRE' vnorg gaur ;
  1341. mbld = 'BLOQUE' 'DEPL' 'DIRE' vnord dror ;
  1342. mblb = 'BLOQUE' 'DEPL' bas ;
  1343. mblp = 'BLOQUE' 'DEPL' (pC 'ET' pD) ;
  1344. bltot = mblh 'ET' mblg 'ET' mbld 'ET' mblb 'ET' mblp ;
  1345. fbltot = 'DEPIMPOSE' mblb ssbas ;
  1346. *
  1347. * Matrice et second membre totaux
  1348. * Résolution
  1349. *
  1350. mtot = rig 'ET' mcnt 'ET' bltot ;
  1351. *mtot = rig 'ET' bltot ;
  1352. ftot = fbltot ;
  1353. *'LISTE' ('EXTRAIRE' mtot 'COMP') ;
  1354. *'LISTE' ('EXTRAIRE' mtot 'COMP' 'DUAL') ;
  1355. *'LISTE' ('EXTRAIRE' ftot 'COMP') ;
  1356. * opti impi 1 ;
  1357. lpass = VRAI ;
  1358. * On recopie la matrice car enchaîner les deux options peut poser
  1359. * problème 
  1360. mtot1 = mtot '*' 1. ; mtot2 = mtot '*' 1. ;
  1361. sol1 = 'KRES' mtot1 ftot 'LDEPE' FAUX 'IMPR' 2 ;
  1362. sol2 = 'KRES' mtot2 ftot 'LDEPE' VRAI 'IMPR' 2 ;
  1363. vit = 'EXCO' ('MOTS' 'UX' 'UY') sol1 ;
  1364. 'SI' graph ;
  1365.    vvit = 'VECT' vit 'DEPL' 'JAUN' ;
  1366.    'TRACER' vvit cmt 'TITR' 'Vitesses' ;
  1367. 'FINSI' ;   
  1368. *
  1369. * Comme l'élimination des dépendances est un préconditionnement,
  1370. * on peut vérifier qu'il marche
  1371. *
  1372. sol3 = 'KRES' mtot1 ftot 'LDEPE' FAUX 'IMPR' 2 ;
  1373. sol4 = 'KRES' mtot2 ftot 'LDEPE' VRAI 'IMPR' 2 ;
  1374. errv = 1.D-14 ;
  1375. err1 = 'MAXIMUM' ('-' sol1 sol3 'ABS') ;
  1376. err2 = 'MAXIMUM' ('-' sol2 sol4 'ABS') ;
  1377. tst1 = '<' err1 errv ; tst2 = '<' err2 errv ;
  1378. tst = tst1 'ET' tst2 ;
  1379. 'MESSAGE' ('CHAINE' 'Test 1 sur la reproductibilité de lelimination')  ;
  1380. 'MESSAGE' '         err1 = ' err1 ' tol=' errv ;
  1381. 'MESSAGE' '         err2 = ' err2 ' tol=' errv  ;
  1382. 'SI' tst ;
  1383.    'MESSAGE' 'Test 1 OK' ;
  1384. 'SINON' ;   
  1385.    'MESSAGE' '!!! Test 1 not passed  ' ;
  1386. 'FINSI' ;
  1387. lpass = lpass 'ET' tst ;
  1388. *
  1389. * Vérification de CONDENSE et EVAPORE
  1390. *
  1391. mtotb = '*' mtot 1.D0 ;
  1392. * On enlève ce test pour ne pas modifier mtot
  1393. * mtotc ftotc ftot1 = 'KOPS' 'CONDENSE' mtot ftot ;
  1394. * solc1 = 'KRES' mtotc ftotc 'LDEPE' FAUX 'IMPR' 2  ;
  1395. * * Vérification de la résolution du système réduit (cf. fiche ano 8186)
  1396. * sol5 = 'KOPS' 'EVAPORE' solc1 mtot ftot ftot1 ;
  1397. mtotc ftotc ftot1 = 'KOPS' 'CONDENSE' mtotb ftot ;
  1398. solc2 = 'KRES' mtotc ftotc 'LDEPE' FAUX 'IMPR' 2  ;
  1399. err0 = '/' ('MAXIMUM' ('-' ftotc ('*' mtotc solc2)) 'ABS')
  1400.            ('MAXIMUM' ftotc 'ABS') ;
  1401. sol6 = 'KOPS' 'EVAPORE' solc2 mtotb ftot ftot1 ;
  1402. errv = 1.D-14 ;
  1403. * err1 = 'MAXIMUM' ('-' sol2 sol5 'ABS') ;
  1404. err2 = 'MAXIMUM' ('-' sol2 sol6 'ABS') ;
  1405. tst0 = '<' err0 errv ;
  1406. * tst1 = '<' err1 errv ;
  1407. tst2 = '<' err2 errv ;
  1408. * tst = tst0 'ET' tst1 'ET' tst2 ;
  1409. tst = tst0 'ET' tst2 ;
  1410. 'MESSAGE' ('CHAINE' 'Test 2 sur CONDENSE et EVAPORE')  ;
  1411. 'MESSAGE' '         err0 = ' err0 ' tol=' errv ;
  1412. * 'MESSAGE' '         err1 = ' err1 ' tol=' errv ;
  1413. 'MESSAGE' '         err2 = ' err2 ' tol=' errv ;
  1414. 'SI' tst ;
  1415.    'MESSAGE' 'Test 2 OK' ;
  1416. 'SINON' ;   
  1417.    'MESSAGE' '!!! Test 2 not passed  ' ;
  1418. 'FINSI' ;
  1419. *'OPTI' 'DONN'  5 ;
  1420. lpass = lpass 'ET' tst ;
  1421. *
  1422. * On vérifie que l'élimination des conditions aux limites sur
  1423. * la vitesse s'est bien faite.
  1424. * Les résultats avec et sans élimination des conditions aux limites ne
  1425. * sont pas si proches car la matrice Vitesse-Pression est mal
  1426. * conditionnée (l'erreur est plus élevée en pression)
  1427. *
  1428. errv = 1.D-13 ; errp = 1.D-11 ; errlx = 4.D-13 ;
  1429. err1 = 'MAXIMUM' ('-' sol1 sol2) 'ABS' 'AVEC' nomvp ;
  1430. err2 = 'MAXIMUM' ('-' sol1 sol2) 'ABS' 'AVEC' nompp ;
  1431. err3 = 'MAXIMUM' ('-' sol1 sol2) 'ABS' 'SANS' ('ET' nomvp nompp) ;
  1432. tst1 = '<' err1 errv ;
  1433. tst2 = '<' err2 errp ;
  1434. tst3 = '<' err3 errlx ;
  1435. tst = tst1 'ET' tst2 'ET' tst3 ;
  1436. 'MESSAGE' ('CHAINE' 'Test 3 sur la solution avec OU sans elimination'
  1437.    ' des conditions aux limites')  ;
  1438. 'MESSAGE' '         err UX  = ' err1 ' tol=' errv  ;
  1439. 'MESSAGE' '         err P   = ' err2 ' tol=' errp  ;
  1440. 'MESSAGE' '         err LX  = ' err3 ' tol=' errlx ;
  1441. 'SI' tst ;
  1442.    'MESSAGE' 'Test 3 OK' ;
  1443. 'SINON' ;   
  1444.    'MESSAGE' '!!! Test 3 not passed  ' ;
  1445. 'FINSI' ;
  1446. lpass = lpass 'ET' tst ;
  1447. *
  1448. *
  1449. * On essaie le raffinement itératif pour rapprocher
  1450. * les solutions avec et sans élimination
  1451. * On choisit nit=1 car sinon on dégrade à nouveau la précision
  1452. *
  1453. nit = 1 ;
  1454. sol7 = sol1 ;
  1455. 'REPETER' it nit ;
  1456.    smb  = '-' ftot ('*' mtot sol7) ;
  1457.    dsol = 'KRES' mtot1 smb 'LDEPE' FAUX ;
  1458.    sol7 = '+' sol7 dsol ;
  1459. 'FIN' it ;
  1460. sol8 = sol2 ;
  1461. 'REPETER' it nit ;
  1462.    smb  = '-' ftot ('*' mtot sol8) ;
  1463.    dsol = 'KRES' mtot2 smb 'LDEPE' VRAI ;
  1464.    sol8 = '+' sol8 dsol ;
  1465. 'FIN' it ;
  1466. errv = 2.D-14 ; errp = 2.D-12 ; errlx = 8.D-14 ;
  1467. err1 = 'MAXIMUM' ('-' sol7 sol8) 'ABS' 'AVEC' nomvp ;
  1468. err2 = 'MAXIMUM' ('-' sol7 sol8) 'ABS' 'AVEC' nompp ;
  1469. err3 = 'MAXIMUM' ('-' sol7 sol8) 'ABS' 'SANS' ('ET' nomvp nompp) ;
  1470. tst1 = '<' err1 errv ;
  1471. tst2 = '<' err2 errp ;
  1472. tst3 = '<' err3 errlx ;
  1473. tst = tst1 'ET' tst2 'ET' tst3 ;
  1474. 'MESSAGE' ('CHAINE' 'Test 4 sur la solution avec OU sans elimination'
  1475.    ' 'ET' raffinement iteratif')  ;
  1476. 'MESSAGE' '         err UX  = ' err1 ' tol=' errv  ;
  1477. 'MESSAGE' '         err P   = ' err2 ' tol=' errp  ;
  1478. 'MESSAGE' '         err LX  = ' err3 ' tol=' errlx ;
  1479. 'SI' tst ;
  1480.    'MESSAGE' 'Test 4 OK' ;
  1481. 'SINON' ;   
  1482.    'MESSAGE' '!!! Test 4 not passed  ' ;
  1483. 'FINSI' ;
  1484. lpass = lpass 'ET' tst ;
  1485. * On choisit ensuite sol8 comme solution de référence
  1486. solref = sol8 ;
  1487. *
  1488. * Méthode de lagrangien augmenté itérée
  1489. *
  1490. * Eliminons d'abord le plus de contraintes possible
  1491. mtotd = '*' mtot 1. ;
  1492. mtotc ftotc ftot1 = 'KOPS' 'CONDENSE' mtotd ftot ;
  1493. * On vérifie que les 'FLX' ont disparus dans la matrice condensée et de
  1494. * la force (en fait elle existe toujours dans la force)
  1495. * MAJ: suite à la fiche ano 8186, on peut rétablir le test
  1496. lp1 = 'EXTRAIRE' mtotc 'COMP' 'DUAL' ;
  1497. lp2 = 'EXTRAIRE' ftotc 'COMP' ;
  1498. tst1 = 'NON' ('EXISTE' lp1 'FLX') ;
  1499. tst2 = 'NON' ('EXISTE' lp2 'FLX') ;
  1500. tst = tst1 'ET' tst2 ;
  1501. *tst = tst1 ;
  1502. 'MESSAGE' ('CHAINE' 'Test 5 sur condensation des LX et FLX')  ;
  1503. 'MESSAGE' ('CHAINE' 'lp1=') ; 'LISTE' lp1 ;
  1504. 'MESSAGE' ('CHAINE' 'lp2=') ; 'LISTE' lp2 ;
  1505. 'MESSAGE' ('CHAINE' '         tst1 = ' tst1)  ;
  1506. 'MESSAGE' ('CHAINE' '         tst2 = ' tst2)  ;
  1507. 'SI' tst ;
  1508.    'MESSAGE' 'Test 5 OK' ;
  1509. 'SINON' ;   
  1510.    'MESSAGE' '!!! Test 5 not passed  ' ;
  1511. 'FINSI' ;
  1512. lpass = lpass 'ET' tst ;
  1513. * Ensuite extraction rigidité et contraintes
  1514. krig  = 'EXTRAIRE' mtotc nomvp nomvd ;
  1515. kcnt  = 'EXTRAIRE' mtotc nomvp ('MOTS' 'LXP') ;
  1516. kcnt2 = 'EXTRAIRE' mtotc ('MOTS' 'LXP') nomvd ;
  1517. * Matrice masse (+ diagonalisé pour la pression)
  1518. kmasp  = GMASS2 _mt tmode 'NPRI' 'PN' 'NDUA' 'PN' ;
  1519. kmaspd = GMASS2 _mt tmode 'NPRI' 'PN' 'FPRI' 1.D0 'NDUA' 'PN' ;
  1520. mas    = 'MAXIMUM' ('RESULT' kmaspd) ;
  1521. dkmasp = 'EXTRAIRE' kmasp 'DIAG' ;
  1522. masd   = 'MAXIMUM' ('RESULT' dkmasp) ;
  1523. 'MESSAGE' ('CHAINE' 'mas=' mas) ;
  1524. 'MESSAGE' ('CHAINE' 'masd=' masd) ;
  1525. * Cela ne paraît pas forcément utile
  1526. *dkmass = '*' dkmasp ('/' mas masd) ;
  1527. dkmass = dkmasp ;
  1528. dkmassi = 'INVERSE' dkmass ;
  1529. *
  1530. * Matrice de pénalisation
  1531. *
  1532. mlagd = 'KOPS' 'CMCT' kcnt2 dkmassi kcnt2 ;
  1533. * Scaling
  1534. drig  = 'EXTRAIRE' krig  'DIAG' ;
  1535. dlagd = 'EXTRAIRE' mlagd 'DIAG' ;
  1536. dsca  = '/' drig dlagd nomvp nomvp nomvp ;
  1537. * Petite vérif de la division nouvellement créée
  1538. drrig = '*' dsca dlagd nomvp nomvp nomvp ;
  1539. errv = 1.D-14 ;
  1540. err1 = 'MAXIMUM' ('-' drrig drig 'ABS') ;
  1541. tst1 = '<' err1 errv ;
  1542. tst = tst1 ;
  1543. 'MESSAGE' ('CHAINE' 'Test 6 sur /')  ;
  1544. 'MESSAGE' '         err1 = ' err1  ;
  1545. 'SI' tst ;
  1546.    'MESSAGE' 'Test 6 OK' ;
  1547. 'SINON' ;   
  1548.    'MESSAGE' '!!! Test 6 not passed  ' ;
  1549. 'FINSI' ;
  1550. lpass = lpass 'ET' tst ;
  1551. *'OPTI' 'ECHO' 1 ;
  1552. *'LISTE' ('MAXIMUM' dsca) ;
  1553. *'LISTE' ('MINIMUM' dsca) ;
  1554. *'OPTION' 'DONN' 5 ;
  1555. * Matrice de pénalisation scalée : on scale kcnt
  1556. **mdsca = 'MANUEL' 'RIGI' ('**' dsca 0.5) ;
  1557. **kcnts = 'KOPS' 'CMCT' kcnt mdsca ;
  1558. **kcnt2s = 'KOPS' 'TRANSPOS' kcnts ;
  1559. *mlagds = '*' mlagd ('MINIMUM' dsca) ;
  1560. *dlagds = 'EXTRAIRE' mlagds 'DIAG' ;
  1561. *dscas  = '/' drig dlagds nomvp nomvp nomvp ;
  1562. *'LISTE' ('MAXIMUM' dscas) ;
  1563. *'LISTE' ('MINIMUM' dscas) ;
  1564. *
  1565. lalphad = 'PROG' 1. 10. 100. ;
  1566. lerrv   = 'PROG' 9.D-2 9.D-4 9.D-6 ;
  1567. lerrp   = 'PROG' 6.D-1 6.D-2 6.D-4 ;
  1568. lerrlx  = 'PROG' 5.D-2 5.D-4 5.D-6 ;
  1569. ncas = 'DIME' lalphad ;
  1570. 'REPETER' iicas ncas ;
  1571.    icas = &iicas ;
  1572.    alphadd = 'EXTRAIRE' lalphad icas ;
  1573.    alphad = '*' alphadd ('MINIMUM' dsca) ;
  1574.    mlagds = '*' mlagd alphad  ;
  1575. * Directe
  1576.    solc = 'KRES' mtotc ftotc 'LDEPE' FAUX ;
  1577. *
  1578.    nit = 2 ;
  1579.    sol9c  = '*' solc 0. ;
  1580.    ftot9c = ftotc ;
  1581.    mtot9c = mtotc ;
  1582.    maug = krig 'ET' mlagds ;
  1583.    'REPETER' it nit ;
  1584.       smb    = '-' ftot9c ('*' mtot9c sol9c) ;
  1585.       smbf = 'EXCO' nomvd smb ;
  1586.       smbg = 'EXCO' nompd smb ;
  1587.       smbf2 = '*' ('*' ('*' dkmassi smbg nompp nompd nompp) kcnt2)
  1588.          alphad ;
  1589.       smbaug = '+' smbf smbf2 ;
  1590.       dsolu tt = 'KRES' maug smbaug 'LDEPE' FAUX
  1591.          'LTIME' VRAI 'IMPR' 1 'TYPINV' 1 ;
  1592. *   'LISTE' tt ;
  1593.       dsolp = '*' ('*' dkmassi ('-' ('*' kcnt dsolu) smbg)
  1594.                        nompp nompd nompp)
  1595.                   ('*' alphad +1.) ;
  1596.       dsol = dsolu '+' dsolp ;
  1597.       sol9c = sol9c '+' dsol ;
  1598.    'FIN' it ;
  1599.    errv  = 'EXTRAIRE' lerrv  icas ;
  1600.    errp  = 'EXTRAIRE' lerrp  icas ;
  1601.    errlx = 'EXTRAIRE' lerrlx icas ;
  1602.    sol9 = 'KOPS' 'EVAPORE' sol9c mtotd ftot ftot1 ;
  1603.    'SI' graph ;
  1604.       vvit2 = 'VECT' sol9 'DEPL' 'ROUG' ;
  1605.       tit = 'CHAINE' 'Vitesses' ' jaun=direct' ' roug=lagrangien alpha='
  1606.          alphadd ;
  1607.       'TRACER' ('ET' vvit vvit2) cmt 'TITR' tit ;
  1608.    'FINSI' ;   
  1609.    err1 = 'MAXIMUM' ('-' sol8 sol9) 'ABS' 'AVEC' nomvp ;
  1610.    err2 = 'MAXIMUM' ('-' sol8 sol9) 'ABS' 'AVEC' nompp ;
  1611.    err3 = 'MAXIMUM' ('-' sol8 sol9) 'ABS' 'SANS' ('ET' nomvp nompp) ;
  1612.    tst1 = '<' err1 errv ;
  1613.    tst2 = '<' err2 errp ;
  1614.    tst3 = '<' err3 errlx ;
  1615.    tst = tst1 'ET' tst2 'ET' tst3 ;
  1616.    'MESSAGE' ('CHAINE' 'Test 7 Lagrangien augmenté' ' '
  1617.       nit ' iterations' ' alpha=' alphadd) ;
  1618.    'MESSAGE' ('CHAINE' '         err UX  = ' err1 ' tol=' errv )  ;
  1619.    'MESSAGE' ('CHAINE' '         err P   = ' err2 ' tol=' errp )  ;
  1620.    'MESSAGE' ('CHAINE' '         err LX  = ' err3 ' tol=' errlx)  ;
  1621.    'SI' tst ;
  1622.       'MESSAGE' 'Test 7 OK' ;
  1623.    'SINON' ;   
  1624.       'MESSAGE' '!!! Test 7 not passed  ' ;
  1625.    'FINSI' ;
  1626.    lpass = lpass 'ET' tst ;
  1627. 'FIN' iicas ;
  1628. *
  1629. 'SAUTER' 2 'LIGNE' ;
  1630. 'SI' lpass ;
  1631.    'MESSAGE' 'Tout sest bien passe' ;
  1632. 'SINON' ;
  1633.    'MESSAGE' 'Il y a eu des erreurs' ;
  1634. 'FINSI' ;   
  1635. 'SAUTER' 2 'LIGNE' ;
  1636. 'SI' ('NON' lpass) ;
  1637.    'ERREUR' 5 ;
  1638. 'FINSI' ;   
  1639. ************************************************************************
  1640. *
  1641. * END OF TEST PART
  1642. * 
  1643. ************************************************************************
  1644. 'FINSI' ;
  1645. 'SI' interact ;
  1646.    'OPTION' 'DONN' 5 'ECHO' 1 ;
  1647. 'FINSI' ;   
  1648. *
  1649.  'FIN' ; 
  1650.  
  1651.  
  1652.  
  1653.  
  1654.  
  1655.  
  1656.  
  1657.  
  1658.  
  1659.  
  1660.  
  1661.  
  1662.  
  1663.  
  1664.  
  1665.  
  1666.  
  1667.  

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