Télécharger la plaquette Cast3M


fr en

Télécharger fcourant.procedur

Retour à la liste

Numérotation des lignes : 

  1. * FCOURANT  PROCEDUR  GOUNAND   19/07/12    21:15:09     10256          
  2. ************************************************************************
  3. * NOM         : FCOURANT
  4. * DESCRIPTION : Calcul la fonction de courant en 2D et 2D Axi
  5. *               par une méthode d'éléments finis moindres carrés
  6. *
  7. * LANGAGE     : GIBIANE-CAST3M
  8. * AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  9. *               mél : gounand@semt2.smts.cea.fr
  10. **********************************************************************
  11. * VERSION    : v1, 22/11/2007, version initiale
  12. * HISTORIQUE : v1, 22/11/2007, création
  13. * HISTORIQUE : 2018/10/15 : 2eme forme en axi +  (fiche 9973)
  14. *              composantes UR et UZ par defaut
  15. *              recodage, utilisation RESO par defaut, utilisation GIBI.ERREUR 
  16. * HISTORIQUE : 2019/07/10 :  choix de la forme normale par défaut
  17. *              sinon sur hsu5 compressible mail 11 pour rho u, on a des
  18. *              recirculations bizarres près de l'axe et des valeurs
  19. *              négatives de flux. (fiche 10256)
  20. * HISTORIQUE :
  21. ************************************************************************
  22. *
  23. * Procedure fonction de courant 2D et 2D axi
  24. *
  25. *             psi est la fonction de courant (streamfunction)
  26. * en 2D, on résout :   dpsi/dx =   u_y
  27. *                      dpsi/dy = - u_x
  28. * 
  29. *             psi est la fonction de courant (Stokes streamfunction)
  30. *             B_\theta = psi / r est le potentiel vecteur (B_r = B_z =0)
  31. *
  32. *
  33. * en 2D axi , on résout soit la forme "regularisee" :
  34. *                         dpsi/dr  =   2 pi r u_z
  35. *                         dpsi/dz  = - 2 pi r u_r
  36. *                       soit la forme "normale" :
  37. *                 ou    :  1/(2 pi r)  dpsi/dr =   u_z
  38. *                          1/(2 pi r)  dpsi/dz = - u_r
  39. *
  40. * Par defaut on resout la forme normale. On ne le met pas dans la
  41. * notice par souci de simplicite.
  42. * Il y a encore une autre forme intermédiaire avec des racines de
  43. * (2 pi r)
  44. *
  45. *
  46. 'DEBPROC' FCOURANT ;
  47. *
  48. 'ARGUMENT' mail*'MAILLAGE' ;
  49. disc = DEADUTIL 'TYPM' mail ;
  50. mquaf = 'CHAN' mail 'QUAF' ;
  51. 'ARGUMENT' un*'CHPOINT' ;
  52. *
  53. * Verif que un est bien defini sur mail
  54. *
  55. mun = 'EXTR' un 'MAIL' ;
  56. pmail = 'CHAN' mail 'POI1' ; pmun = 'CHAN' mun 'POI1' ;
  57. minter = 'INTE' pmail pmun ;
  58. mnoun = 'DIFF' pmail minter ;
  59. 'SI' ('>' ('NBEL' mnoun) 0) ;
  60.    nnoex = 'POIN' mnoun 1 ;
  61. * 771 0 ATTENTION !!! Le CHPOINT ne definit pas de valeur au noeud %i1
  62.    'ERRE' 771 'AVEC' nnoex ;
  63. 'FINS' ;
  64. *
  65. 'ARGUMENT' rigblo/'RIGIDITE' ;
  66. lrb = 'EXISTE' rigblo ;
  67. 'SI' lrb ;
  68.    'ARGUMENT' chblo/'CHPOINT' ;
  69.    lcb = 'EXISTE' chblo ;
  70. 'SINON' ;
  71.    lcb = FAUX ;   
  72. 'FINSI' ;
  73. *
  74. 'SI' ('NON' ('EXISTE' rigblo)) ;
  75.    pref  = 'POIN' ('CHANGER' mail 'POI1') 1 ;
  76.    rigblo  = 'BLOQ' 'T' pref ;
  77. 'FINSI' ;
  78. *
  79. 'ARGUMENT' rvm/'TABLE' ;
  80. 'SI' ('NON' ('EXISTE' rvm)) ;
  81.    typsolv = 0 ;
  82. 'SINO' ;   
  83.    typsolv = 1 ;
  84. 'FINSI' ;   
  85. *
  86. dim = 'VALEUR' 'DIME' ;
  87. *
  88. 'SI' ('NEG' dim 2) ;
  89. * 709 2
  90. * Fonction indisponible en dimension %i1.
  91.    'ERREUR' 709 'AVEC' dim ;
  92. 'FINSI' ;
  93. *
  94. laxi = 'EGA' ('VALEUR' 'MODE') 'AXIS' ;
  95. *
  96. 'SI' laxi ;
  97.    ncr = 'MOTS' 'UR' ; ncz = 'MOTS' 'UZ' ;
  98. * Forme regularisee ou normale    
  99.    'ARGU' mform/'MOT' ;
  100.    'SI' ('NON' ('EXIS' mform)) ;
  101.       mform = 'CHAI' 'NORM' ;
  102.    'FINS' ;
  103.    lform = 'MOTS' 'REGU' 'NORM' ;
  104.    iform = 'POSI' mform 'DANS' lform ;
  105. *1052 2
  106. *Mot-cle incorrect "%M1:4". Voici la liste des valeurs admises :
  107. *1052 2
  108. *%M5:40
  109.    'SI' ('EGA' iform 0) ;
  110.       'ERRE' 1052 'AVEC' ('CHAI' mform 'REGU NORM') ;
  111.    'FINS' ;
  112. * iform = 0 : régularisée   
  113. * iform = 1 : standard
  114. *dbg   'SI' ('EGA' iform 1) ;
  115. *dbg      'MESS' 'FCOURANT Axi : forme regularisee' ;
  116. *dbg   'SINO' ;      
  117. *dbg      'MESS' 'FCOURANT Axi : forme normale' ;
  118. *dbg   'FINS' ;      
  119. * Si on ne trouve pas de composantes UR ou UZ dans le champ un en entree
  120. * mais des composantes UX,UY alors on les autorise car c'est la
  121. * convention de nommage avec EQEX et EXEC...
  122.    ncun = 'EXTR' un 'COMP' ;
  123.    'SI' ('NON' ('EXIS' (ncr 'ET' ncz) ncun 'OU')) ;
  124.       'SI' ('EXIS' ('MOTS' 'UX' 'UY') ncun 'ET') ;
  125.          ncr = 'MOTS' 'UX' ; ncz = 'MOTS' 'UY' ;
  126.       'FINS' ;
  127.    'FINS' ;      
  128. *
  129.    cdpr  = '*' ('COORDONNEE' 1 mail) ('*' PI 2.D0) ;
  130. *
  131.    'SI' ('EGA' iform 1) ;
  132. * Cas 2D Axi Forme regularisee
  133.       numop  = 2 ; numvar = 1 ; numder = 2 ;
  134.       numdat = 1 ; numcof = 1 ;
  135. *
  136.       A = ININLIN numop numvar numdat numcof numder ;
  137.       A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'T' ;
  138.       A . 'VAR' . 1 . 'DISC'   = disc ;
  139. *
  140.       A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'DUMM' ;
  141.       A . 'DAT' . 1 . 'DISC'   = 'CSTE' ;
  142.       A . 'DAT' . 1 . 'VALEUR' = -1. ;
  143.       A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  144.       A . 'COF' . 1 . 'LDAT'   = 'LECT' 1 ;
  145. *      
  146.       A . 1 . 1 . 1 = 'LECT' ;
  147.       A . 2 . 1 . 2 = 'LECT' 1 ;
  148. *
  149.       numdat = 2 ; numcof = 2 ;
  150.       B = ININLIN numop numvar numdat numcof numder ;
  151.       B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'Q' ;
  152.       B . 'VAR' . 1 . 'DISC'   = disc ;
  153. *     
  154.       B . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'DUMM' ;
  155.       B . 'DAT' . 1 . 'DISC'   = 'CSTE' ;
  156.       B . 'DAT' . 1 . 'VALEUR' = -1. ;
  157.       B . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  158.       B . 'COF' . 1 . 'LDAT'   = 'LECT' 1 ;
  159. *      
  160.       B . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  161.       B . 'DAT' . 2 . 'DISC'   = disc ;
  162.       B . 'DAT' . 2 . 'VALEUR' = cdpr ;
  163.       B . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
  164.       B . 'COF' . 2 . 'LDAT'   = 'LECT' 2 ;
  165. *   
  166.       B . 1 . 1 . 1 = 'LECT' 2 ;
  167.       B . 2 . 1 . 2 = 'LECT' 1 2 ;
  168. *
  169.       numvar = 2 ; numdat = 1 ; numcof = 1 ;
  170.       C = ININLIN numop numvar numdat numcof numder ;
  171.       C . 'VAR' . 1 . 'NOMDDL' = ncr ;
  172.       C . 'VAR' . 1 . 'DISC'   = disc ;
  173.       C . 'VAR' . 1 . 'VALEUR' = un ;
  174.       C . 'VAR' . 2 . 'NOMDDL' = ncz ;
  175.       C . 'VAR' . 2 . 'DISC'   = disc ;
  176.       C . 'VAR' . 2 . 'VALEUR' = un ;
  177.       C . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  178.       C . 'DAT' . 1 . 'DISC'   = disc ;
  179.       C . 'DAT' . 1 . 'VALEUR' = cdpr ;
  180.       C . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  181.       C . 'COF' . 1 . 'LDAT'   = 'LECT' 1 ;
  182. *
  183.       C . 1 . 2 . 0 = 'LECT' 1 ;
  184.       C . 2 . 1 . 0 = 'LECT' 1 ;
  185.    'SINO' ;
  186. * Cas 2D Axi Forme normale
  187.       numop  = 2 ; numvar = 1 ; numder = 2 ;
  188.       numdat = 2 ; numcof = 2 ;
  189. *
  190.       A = ININLIN numop numvar numdat numcof numder ;
  191.       A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'T' ;
  192.       A . 'VAR' . 1 . 'DISC'   = disc ;
  193. *      
  194.       A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'DUMM' ;
  195.       A . 'DAT' . 1 . 'DISC'   = 'CSTE' ;
  196.       A . 'DAT' . 1 . 'VALEUR' = -1. ;
  197.       A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  198.       A . 'COF' . 1 . 'LDAT'   = 'LECT' 1 ;
  199. *
  200.       A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  201.       A . 'DAT' . 2 . 'DISC'   = disc ;
  202.       A . 'DAT' . 2 . 'VALEUR' = cdpr ;
  203.       A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
  204.       A . 'COF' . 2 . 'LDAT'   = 'LECT' 2 ;
  205.       A . 1 . 1 . 1 = 'LECT' -2 ;
  206.       A . 2 . 1 . 2 = 'LECT' 1 -2 ;
  207. *
  208.       numdat = 1 ; numcof = 1 ;
  209.       B = ININLIN numop numvar numdat numcof numder ;
  210.       B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'Q' ;
  211.       B . 'VAR' . 1 . 'DISC'   = disc ;
  212. *     
  213.       B . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'DUMM' ;
  214.       B . 'DAT' . 1 . 'DISC'   = 'CSTE' ;
  215.       B . 'DAT' . 1 . 'VALEUR' = -1. ;
  216.       B . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  217.       B . 'COF' . 1 . 'LDAT'   = 'LECT' 1 ;
  218. *   
  219.       B . 1 . 1 . 1 = 'LECT'   ;
  220.       B . 2 . 1 . 2 = 'LECT' 1 ;
  221. *
  222.       numvar = 2 ; numdat = 0 ; numcof = 0 ;
  223.       C = ININLIN numop numvar numdat numcof numder ;
  224.       C . 'VAR' . 1 . 'NOMDDL' = ncr ;
  225.       C . 'VAR' . 1 . 'DISC'   = disc ;
  226.       C . 'VAR' . 1 . 'VALEUR' = un ;
  227.       C . 'VAR' . 2 . 'NOMDDL' = ncz ;
  228.       C . 'VAR' . 2 . 'DISC'   = disc ;
  229.       C . 'VAR' . 2 . 'VALEUR' = un ;
  230. *
  231.       C . 1 . 2 . 0 = 'LECT' ;
  232.       C . 2 . 1 . 0 = 'LECT' ;
  233.    'FINS' ;
  234. *
  235. 'SINON'  ;
  236. * Cas 2D PLAN
  237.    ncx = 'MOTS' 'UX' ; ncy = 'MOTS' 'UY' ;
  238. *
  239.    numop  = 2 ; numvar = 1 ; numder = 2 ;
  240.    numdat = 1 ; numcof = 1 ;
  241. *
  242.    A = ININLIN numop numvar numdat numcof numder ;
  243.    A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'T' ;
  244.    A . 'VAR' . 1 . 'DISC'   = disc ;
  245.    A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'DUMM' ;
  246.    A . 'DAT' . 1 . 'DISC'   = 'CSTE' ;
  247.    A . 'DAT' . 1 . 'VALEUR' = -1 ;
  248.    A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  249.    A . 'COF' . 1 . 'LDAT'   = 'LECT' 1 ;
  250.  
  251.    A . 1 . 1 . 1 = 'LECT' ;
  252.    A . 2 . 1 . 2 = 'LECT' 1 ;
  253. *
  254.    B = ININLIN numop numvar numdat numcof numder ;
  255.    B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'Q' ;
  256.    B . 'VAR' . 1 . 'DISC'   = disc ;
  257.    B . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'DUMM' ;
  258.    B . 'DAT' . 1 . 'DISC'   = 'CSTE' ;
  259.    B . 'DAT' . 1 . 'VALEUR' = -1 ;
  260.    B . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  261.    B . 'COF' . 1 . 'LDAT'   = 'LECT' 1 ;
  262.  
  263.    B . 1 . 1 . 1 = 'LECT' ;
  264.    B . 2 . 1 . 2 = 'LECT' 1 ;
  265. *
  266.    numvar = 2 ; numdat = 0 ; numcof = 0 ;
  267.    C = ININLIN numop numvar numdat numcof numder ;
  268.    C . 'VAR' . 1 . 'NOMDDL' = ncx ;
  269.    C . 'VAR' . 1 . 'DISC'   = disc ;
  270.    C . 'VAR' . 1 . 'VALEUR' = un ;
  271.    C . 'VAR' . 2 . 'NOMDDL' = ncy ;
  272.    C . 'VAR' . 2 . 'DISC'   = disc ;
  273.    C . 'VAR' . 2 . 'VALEUR' = un ;
  274. *
  275.    C . 1 . 2 . 0 = 'LECT' ;
  276.    C . 2 . 1 . 0 = 'LECT' ;
  277. 'FINSI' ;
  278. *
  279. mat  = NLINP mquaf disc A B 'GAU7' ;
  280. * Remettre la symetrie si necessaire
  281. * Au 2018/10/19 : le test de symetrie de RESO est assez contraignant
  282. * pour que cela ne marche pas !
  283. *lpri = 'EXTR' mat 'COMP' ;
  284. *ldua = 'EXTR' mat 'COMP' 'DUAL' ;
  285. *mat = 'CHAN' 'INCO' mat lpri lpri ldua ldua 'SYME' ;
  286. *old mat = 'KOPS' 'RIMA' mat ;
  287. smb  = NLINP mquaf disc C B 'GAU7' ;
  288. mtot = mat 'ET' rigblo ;
  289. ftot = smb ;
  290. 'SI' lcb ;
  291.    ftot = 'ET' ftot chblo ;
  292. 'FINS' ;
  293. *
  294. 'SI' ('EGA' typsolv 0) ;
  295.    solt = 'RESO' mtot ftot ;
  296. 'SINO' ;
  297.    solt = 'KRES' mtot 'TYPI' rvm 'SMBR' ftot ;
  298. 'FINS' ;   
  299. psi = 'EXCO' 'T' solt 'PSI' ;
  300. 'RESPRO' psi ;
  301. *
  302. 'FINPROC' ;
  303.  
  304.  
  305.  
  306.  

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