Télécharger la plaquette Cast3M


fr en

Télécharger fcourant.procedur

Retour à la liste

Numérotation des lignes : 

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

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