$$$$ FCOURANT * FCOURANT PROCEDUR GOUNAND 26/01/12 21:15:03 12448 ************************************************************************ * NOM : FCOURANT * DESCRIPTION : Calcul la fonction de courant en 2D et 2D Axi * par une méthode d'éléments finis moindres carrés * * LANGAGE : GIBIANE-CAST3M * AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF) * mél : gounand@semt2.smts.cea.fr ********************************************************************** * VERSION : v1, 22/11/2007, version initiale * HISTORIQUE : v1, 22/11/2007, création * HISTORIQUE : 2018/10/15 : 2eme forme en axi + (fiche 9973) * composantes UR et UZ par defaut * recodage, utilisation RESO par defaut, utilisation GIBI.ERREUR * HISTORIQUE : 2019/07/10 : choix de la forme normale par défaut * sinon sur hsu5 compressible mail 11 pour rho u, on a des * recirculations bizarres près de l'axe et des valeurs * négatives de flux. (fiche 10256) * HISTORIQUE : 2025/12/29, simplification donnees pour NLIN * HISTORIQUE : * HISTORIQUE : ************************************************************************ * * Procedure fonction de courant 2D et 2D axi * * psi est la fonction de courant (streamfunction) * en 2D, on résout : dpsi/dx = u_y * dpsi/dy = - u_x * * psi est la fonction de courant (Stokes streamfunction) * B_\theta = psi / r est le potentiel vecteur (B_r = B_z =0) * * * en 2D axi , on résout soit la forme "regularisee" : * dpsi/dr = 2 pi r u_z * dpsi/dz = - 2 pi r u_r * soit la forme "normale" : * ou : 1/(2 pi r) dpsi/dr = u_z * 1/(2 pi r) dpsi/dz = - u_r * * Par defaut on resout la forme normale. On ne le met pas dans la * notice par souci de simplicite. * Il y a encore une autre forme intermédiaire avec des racines de * (2 pi r) * * 'DEBPROC' FCOURANT ; * 'ARGUMENT' mail*'MAILLAGE' ; disc = DEADUTIL 'TYPM' mail ; 'ARGUMENT' un*'CHPOINT' ; * * Verif que un est bien defini sur mail * mun = 'EXTR' un 'MAIL' ; pmail = 'CHAN' mail 'POI1' ; pmun = 'CHAN' mun 'POI1' ; minter = 'INTE' pmail pmun ; mnoun = 'DIFF' pmail minter ; 'SI' ('>' ('NBEL' mnoun) 0) ; nnoex = 'POIN' mnoun 1 ; * 771 0 ATTENTION !!! Le CHPOINT ne definit pas de valeur au noeud %i1 'ERRE' 771 'AVEC' nnoex ; 'FINS' ; * 'ARGUMENT' rigblo/'RIGIDITE' ; lrb = 'EXISTE' rigblo ; 'SI' lrb ; 'ARGUMENT' chblo/'CHPOINT' ; lcb = 'EXISTE' chblo ; 'SINON' ; lcb = FAUX ; 'FINSI' ; * 'SI' ('NON' ('EXISTE' rigblo)) ; pref = 'POIN' ('CHANGER' mail 'POI1') 1 ; rigblo = 'BLOQ' 'T' pref ; 'FINSI' ; * 'ARGUMENT' rvm/'TABLE' ; 'SI' ('NON' ('EXISTE' rvm)) ; typsolv = 0 ; 'SINO' ; typsolv = 1 ; 'FINSI' ; * dim = 'VALEUR' 'DIME' ; * 'SI' ('NEG' dim 2) ; * 709 2 * Fonction indisponible en dimension %i1. 'ERREUR' 709 'AVEC' dim ; 'FINSI' ; * laxi = 'EGA' ('VALEUR' 'MODE') 'AXIS' ; * 'SI' laxi ; ncr = 'MOTS' 'UR' ; ncz = 'MOTS' 'UZ' ; * Forme regularisee ou normale 'ARGU' mform/'MOT' ; 'SI' ('NON' ('EXIS' mform)) ; mform = 'CHAI' 'NORM' ; 'FINS' ; lform = 'MOTS' 'REGU' 'NORM' ; iform = 'POSI' mform 'DANS' lform ; *1052 2 *Mot-cle incorrect "%M1:4". Voici la liste des valeurs admises : *1052 2 *%M5:40 'SI' ('EGA' iform 0) ; 'ERRE' 1052 'AVEC' mform 'REGU NORM' ; 'FINS' ; * iform = 0 : régularisée * iform = 1 : standard *dbg 'SI' ('EGA' iform 1) ; *dbg 'MESS' 'FCOURANT Axi : forme regularisee' ; *dbg 'SINO' ; *dbg 'MESS' 'FCOURANT Axi : forme normale' ; *dbg 'FINS' ; * Si on ne trouve pas de composantes UR ou UZ dans le champ un en entree * mais des composantes UX,UY alors on les autorise car c'est la * convention de nommage avec EQEX et EXEC... ncun = 'EXTR' un 'COMP' ; 'SI' ('NON' ('EXIS' (ncr 'ET' ncz) ncun 'OU')) ; 'SI' ('EXIS' ('MOTS' 'UX' 'UY') ncun 'ET') ; ncr = 'UX' ; ncz = 'UY' ; 'FINS' ; 'FINS' ; * cdpr = '*' ('COORDONNEE' 1 mail) ('*' PI 2.D0) ; * 'SI' ('EGA' iform 1) ; * Cas 2D Axi Forme regularisee numop = 2 ; numvar = 1 ; numder = 2 ; numdat = 1 ; numcof = 1 ; * A = ININLIN numop numvar numdat numcof numder ; A . 'VAR' . 1 . 'NOMDDL' = 'T' ; A . 'VAR' . 1 . 'DISC' = disc ; * A . 'DAT' . 1 . 'VALEUR' = -1. ; A . 'COF' . 1 . 'LDAT' = 1 ; * A . 1 . 1 . 1 = 0 ; A . 2 . 1 . 2 = 1 ; * numdat = 2 ; numcof = 2 ; B = ININLIN numop numvar numdat numcof numder ; B . 'VAR' . 1 . 'NOMDDL' = 'Q' ; B . 'VAR' . 1 . 'DISC' = disc ; * B . 'DAT' . 1 . 'VALEUR' = -1. ; B . 'COF' . 1 . 'LDAT' = 1 ; * B . 'DAT' . 2 . 'NOMDDL' = 'SCAL' ; B . 'DAT' . 2 . 'DISC' = disc ; B . 'DAT' . 2 . 'VALEUR' = cdpr ; B . 'COF' . 2 . 'LDAT' = 2 ; * B . 1 . 1 . 1 = 2 ; B . 2 . 1 . 2 = 'LECT' 1 2 ; * numvar = 2 ; numdat = 1 ; numcof = 1 ; C = ININLIN numop numvar numdat numcof numder ; C . 'VAR' . 1 . 'NOMDDL' = ncr ; C . 'VAR' . 1 . 'DISC' = disc ; C . 'VAR' . 1 . 'VALEUR' = un ; C . 'VAR' . 2 . 'NOMDDL' = ncz ; C . 'VAR' . 2 . 'DISC' = disc ; C . 'VAR' . 2 . 'VALEUR' = un ; C . 'DAT' . 1 . 'NOMDDL' = 'SCAL' ; C . 'DAT' . 1 . 'DISC' = disc ; C . 'DAT' . 1 . 'VALEUR' = cdpr ; C . 'COF' . 1 . 'LDAT' = 1 ; * C . 1 . 2 . 0 = 1 ; C . 2 . 1 . 0 = 1 ; 'SINO' ; * Cas 2D Axi Forme normale numop = 2 ; numvar = 1 ; numder = 2 ; numdat = 2 ; numcof = 2 ; * A = ININLIN numop numvar numdat numcof numder ; A . 'VAR' . 1 . 'NOMDDL' = 'T' ; A . 'VAR' . 1 . 'DISC' = disc ; * A . 'DAT' . 1 . 'VALEUR' = -1. ; A . 'COF' . 1 . 'LDAT' = 1 ; * A . 'DAT' . 2 . 'NOMDDL' = 'SCAL' ; A . 'DAT' . 2 . 'DISC' = disc ; A . 'DAT' . 2 . 'VALEUR' = cdpr ; A . 'COF' . 2 . 'LDAT' = 2 ; A . 1 . 1 . 1 = -2 ; A . 2 . 1 . 2 = 'LECT' 1 -2 ; * numdat = 1 ; numcof = 1 ; B = ININLIN numop numvar numdat numcof numder ; B . 'VAR' . 1 . 'NOMDDL' = 'Q' ; B . 'VAR' . 1 . 'DISC' = disc ; * B . 'DAT' . 1 . 'VALEUR' = -1. ; B . 'COF' . 1 . 'LDAT' = 1 ; * B . 1 . 1 . 1 = 0 ; B . 2 . 1 . 2 = 1 ; * numvar = 2 ; numdat = 0 ; numcof = 0 ; C = ININLIN numop numvar numdat numcof numder ; C . 'VAR' . 1 . 'NOMDDL' = ncr ; C . 'VAR' . 1 . 'DISC' = disc ; C . 'VAR' . 1 . 'VALEUR' = un ; C . 'VAR' . 2 . 'NOMDDL' = ncz ; C . 'VAR' . 2 . 'DISC' = disc ; C . 'VAR' . 2 . 'VALEUR' = un ; * C . 1 . 2 . 0 = 0 ; C . 2 . 1 . 0 = 0 ; 'FINS' ; * 'SINON' ; * Cas 2D PLAN ncx = 'UX' ; ncy = 'UY' ; * numop = 2 ; numvar = 1 ; numder = 2 ; numdat = 1 ; numcof = 1 ; * A = ININLIN numop numvar numdat numcof numder ; A . 'VAR' . 1 . 'NOMDDL' = 'T' ; A . 'VAR' . 1 . 'DISC' = disc ; A . 'DAT' . 1 . 'VALEUR' = -1 ; A . 'COF' . 1 . 'LDAT' = 1 ; A . 1 . 1 . 1 = 0 ; A . 2 . 1 . 2 = 1 ; * B = ININLIN numop numvar numdat numcof numder ; B . 'VAR' . 1 . 'NOMDDL' = 'Q' ; B . 'VAR' . 1 . 'DISC' = disc ; B . 'DAT' . 1 . 'VALEUR' = -1 ; B . 'COF' . 1 . 'LDAT' = 1 ; B . 1 . 1 . 1 = 0 ; B . 2 . 1 . 2 = 1 ; * numvar = 2 ; numdat = 0 ; numcof = 0 ; C = ININLIN numop numvar numdat numcof numder ; C . 'VAR' . 1 . 'NOMDDL' = ncx ; C . 'VAR' . 1 . 'DISC' = disc ; C . 'VAR' . 1 . 'VALEUR' = un ; C . 'VAR' . 2 . 'NOMDDL' = ncy ; C . 'VAR' . 2 . 'DISC' = disc ; C . 'VAR' . 2 . 'VALEUR' = un ; * C . 1 . 2 . 0 = 0 ; C . 2 . 1 . 0 = 0 ; 'FINSI' ; * mat = NLINP mail disc A B 'GAU7' ; * Remettre la symetrie si necessaire * Au 2018/10/19 : le test de symetrie de RESO est assez contraignant * pour que cela ne marche pas ! *lpri = 'EXTR' mat 'COMP' ; *ldua = 'EXTR' mat 'COMP' 'DUAL' ; *mat = 'CHAN' 'INCO' mat lpri lpri ldua ldua 'SYME' ; *old mat = 'KOPS' 'RIMA' mat ; smb = NLINP mail disc C B 'GAU7' ; mtot = mat 'ET' rigblo ; ftot = smb ; 'SI' lcb ; ftot = 'ET' ftot chblo ; 'FINS' ; * 'SI' ('EGA' typsolv 0) ; solt = 'RESO' mtot ftot ; 'SINO' ; solt = 'KRES' mtot 'TYPI' rvm 'SMBR' ftot ; 'FINS' ; psi = 'EXCO' 'T' solt 'PSI' ; 'RESPRO' psi ; * 'FINPROC' ;