1 : C EF0KON SOURCE MAGN 10/05/18 21:15:26 6675 2 : C************************************************************** 3 : C============================================================== 4 : C Cas de l'equation convection diffusion C 5 : C formulation EFM (elements finis mixte) C 6 : C Auteur : M. MARIN C 7 : C C 8 : C IERRKON : indicatif d'erreur = 1 si il y a une erreur C 9 : C = 0 sinon C 10 : C Cette routine cree les matrices elementaires pour un C 11 : C terme u*teta ou : C 12 : C u : CHPOINT vectoriel au sommet C 13 : C teta : CHPOINT scalaire au centre C 14 : C C 15 : C============================================================== 16 : C************************************************************** 17 : 18 : SUBROUTINE EF0KON(KIZX,IZTU1,IZTGG1,IZTGG2,IZTGG3,TYPC,IERRKON, 19 : & MELEME,MTABZ,NOMI,IKR,MLENTI,IAXI,NOMII,MZPHI,TYPCFI, 20 : & IZTCO,NELZ,IKU,IKM,AIMPL,IDCEN,MLENT1,DT) 21 : IMPLICIT INTEGER(I-N) 22 : IMPLICIT REAL*8 (A-H,O-Z) 23 : 24 : 25 : -INC CCVQUA4 26 : -INC CCOPTIO 27 : -INC CCGEOME 28 : -INC SMCOORD 29 : -INC SMLENTI 30 : -INC SMELEME 31 : POINTEUR MELEMC.MELEME,MELEMS.MELEME 32 : -INC SMCHAML 33 : -INC SMCHPOI 34 : POINTEUR IZTU1.MPOVAL,MZPHI.MPOVAL,IZTCO.MPOVAL 35 : POINTEUR IZTGG1.MPOVAL,IZTGG2.MPOVAL,IZTGG3.MPOVAL 36 : -INC SIZFFB 37 : -INC SMMATRIK 38 : -INC SMLMOTS 39 : -INC SMTABLE 40 : POINTEUR KIZX.MTABLE,MTABZ.MTABLE 41 : 42 : C=========================================== 43 : C Variables locale ou argument C 44 : CHARACTER*8 TYPC,NOMI,NOM,NOM0,NOMII,TYPCFI 45 : INTEGER IERRKON,IKR,IKU,IKM 46 : C=========================================== 47 : 48 : c WRITE(6,*) KIZX,IZTU1,IZTGG1,IZTGG2,IZTGG3,TYPC,IERRKON, 49 : c & MELEM,MTABZ,NOMI,IKR,MLENTI,IAXI,NOMII,MZPHI,TYPCFI, 50 : c & IZTCO,NELZ,IKU,IKM,AIMPL 51 : C Validation du passage d'argument 52 : 53 : IERRKON=0 54 : 55 : C Initialisation pour le nombre total final d'elements (TRI3+QUA4+...) 56 : C sur le domaine ou s'applique KONV 57 : NUTOEL=0 58 : c WRITE(6,*) 'MELEME=',MELEME 59 : 60 : C WRITE(6,*) 'Type du Coef 1:',IKR 61 : 62 : C On teste si la variable primale est de type CENTRE 63 : IF (TYPC.NE.'CENTRE ') THEN 63 : CALL ERRKON(1,IERRKON) 64 : RETURN 65 : END IF 66 : 67 : C Test si la variable primale est de type scalaire 68 : IF (IZTU1.VPOCHA(/2).NE.1) THEN 68 : CALL ERRKON(5,IERRKON) 69 : RETURN 70 : END IF 71 : 72 : C On teste si la variable duale est de type SOMMET 73 : IF (TYPCFI.NE.'SOMMET ') THEN 73 : CALL ERRKON(4,IERRKON) 74 : RETURN 75 : END IF 76 : 77 : C On teste si la variable duale est un vecteur 78 : IF ((MZPHI.VPOCHA(/2).EQ.1).AND.(IDIM.NE.1)) THEN 78 : CALL ERRKON(6,IERRKON) 79 : RETURN 80 : END IF 81 : 82 : C On teste pour savoir si RO est bien un scalaire 83 : IF (IKR.NE.1) THEN 83 : CALL ERRKON(7,IERRKON) 84 : RETURN 85 : END IF 86 : C On recupere le nombre de composante de la vitesse 87 : NBME=IDIM 88 : c WRITE(6,*) 'Nombre de composante de UN :',NBME 89 : 90 : C On active le segment contenant le maillage de la zone ou 91 : C KONV s'applique et on initialise le nombre de sous-objet. 92 : C ------------------------------------------------------ 93 : C NBSOUS et NBME servent a initialiser le segment IMATRI 94 : SEGACT MELEME 95 : NBSOUS=MELEME.LISOUS(/1) 96 : c WRITE(6,*) 'NBSOUS=',NBSOUS 97 : IF (NBSOUS.EQ.0) NBSOUS=1 98 : 99 : C Initialisation de la matrice elementaire de l'operateur 100 : NRIGE=7 101 : NKID=9 102 : NKMT=7 103 : NMATRI=1 104 : SEGINI MATRIK 105 : 106 : 107 : C On recupere le maillage des centres 107 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 108 : 109 : C On recupere le maillage des sommets 109 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 110 : 111 : C Inconnue primale 112 : IRIGEL(1,1)=MELEMC 113 : 114 : C Inconnue duale 115 : IRIGEL(2,1)=MELEME 116 : 117 : C Matrice rectangulaire 118 : IRIGEL(7,1)=3 119 : 120 : SEGINI IMATRI 121 : IRIGEL(4,1)= IMATRI 122 : 123 : C Intialisation des supports inconnue primale et duale 124 : C pour le segment IMATRI 125 : KSPGP=MELEMC 126 : KSPGD=MELEMS 127 : 128 : DO I=1,NBME 129 : WRITE(NOM,FMT='(I1,A3)')I,NOMII(1:3) 130 : LISPRI(I)=NOMI(1:4)//' ' 131 : LISDUA(I)=NOM(1:4)//' ' 132 : END DO 133 : 134 : C On recupere le segment qui contient int(Ni) 134 : CALL LEKTAB(MTABZ,'XXPSOML',MCHELM) 135 : 136 : IF (MCHELM.EQ.0) THEN 136 : CALL ERRKON(2,IERRKON) 137 : RETURN 138 : END IF 139 : 140 : SEGACT MCHELM 141 : 142 : C==========Boucle principale==============C 143 : C On effectue cette boucle pour tous les 144 : C sous-objets 145 : 146 : 147 : DO L=1,NBSOUS 148 : IPT1=MELEME 149 : SEGACT MELEME 150 : 151 : C Si NBSOUS > 1 on agit sur les sous-objets LISOUS 152 : IF (NBSOUS.NE.1) IPT1=LISOUS(L) 153 : SEGACT IPT1 154 : 155 : c WRITE(6,*) 'IPT1=',IPT1,'NBSOUS=',NBSOUS 156 : C On rempli NOM0 du veritable nom des elements traites (TRI3, QUA4,...) 157 : NOM0=NOMS(IPT1.ITYPEL)//' ' 158 : 159 : C Cette routine cree des objets contenant les fonctions 160 : C de forme lie aux elements traites 161 : C Le segment IZFFM contient les fonctions de forme et les gradients 162 : C a chaque point de Gauss 162 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 163 : 164 : SEGACT IZFFM 165 : IZHR=KZHR(1) 166 : SEGACT IZHR*MOD 167 : 168 : C Dimension de l'element de reference 169 : NES=GR(/1) 170 : 171 : C Nombre de points de Gauss 172 : NPG=GR(/3) 173 : c WRITE(6,*) 'Nbre de points de Gauss :' , NPG 174 : 175 : C NP est le nombre de points de l'element de reference 176 : C MP est la troisieme dimension du tableau contenant les matrices 177 : C elementaire de chaque elements (segment IZAFM). C'est le 178 : C nombre de point par element 179 : NP = 1 180 : MP = IPT1.NUM(/1) 181 : C NBEL est le nombre d'element de type NOM0 182 : NBEL=IPT1.NUM(/2) 183 : 184 : C On declare un segment IZAFM par composante 185 : SEGINI IPM1 186 : LIZAFM(L,1)=IPM1 187 : 188 : IF (IDIM.EQ.2) THEN 189 : SEGINI IPM2 190 : LIZAFM(L,2)=IPM2 191 : ELSE IF (IDIM.EQ.3) THEN 192 : SEGINI IPM2 193 : SEGINI IPM3 194 : LIZAFM(L,2)=IPM2 195 : LIZAFM(L,3)=IPM3 196 : END IF 197 : 198 : C Nombre total de points NTP 199 : NPT=IZTGG2.VPOCHA(/1) 200 : 201 : C On recupere le tableau contenant les integrales 202 : C de fonction test 203 : MCHAML=ICHAML(L) 204 : SEGACT MCHAML 205 : MELVAL=IELVAL(1) 206 : SEGACT MELVAL 207 : SEGACT MELEMC 208 : C ****** 208 : CALL XCNEF0(NBEL,NP,MP,IPM1.AM,IZTGG1.VPOCHA,IZTGG2 209 : & .VPOCHA,IZTGG3.VPOCHA,NPT,IDIM,IDCEN,XYZ, 210 : & NUTOEL,XCOOR,IPT1.NUM,MLENTI.LECT,IPM2.AM,IPM3.AM, 211 : & FN,GR,PG,HR,PGSQ,RPG,NES,NPG,IAXI,VELCHE, 212 : & NBME,IZTGG3.VPOCHA,IZTCO.VPOCHA,NELZ, 213 : & IKR,IKU,IKM,IZTU1.VPOCHA,AIMPL,MLENT1.LECT,DT, 214 : & MELEMC) 215 : C ****** 216 : SEGDES MELEMC 217 : SEGDES IPT1,IPM1,IPM2 218 : IF (IDIM.EQ.3) THEN 219 : SEGDES IPM3 220 : END IF 221 : SEGDES MCHAML,MELVAL 222 : SEGSUP IZHR,IZFFM 223 : 224 : 225 : C On met a jour le nombre d'element total (quelquesoit son type) 226 : NUTOEL=NUTOEL+NBEL 227 : END DO 228 : 229 : C==========Fin Boucle principale==========C 230 : 231 : C On desactive les segments 232 : 233 : SEGDES MELEME 234 : SEGDES MCHELM 235 : 236 : C Si IKR = 0 alors IZTGG1 (coef 1) est un CHPOINT 237 : C que l'on desactive. 238 : IF(IKR.EQ.0)THEN 239 : SEGDES IZTGG1 240 : ENDIF 241 : 242 : C On ecrit dans la table KONV la matrice elementaire 243 : C de l'operateur 243 : CALL ECMO(KIZX,'MATELM','MATRIK',MATRIK) 244 : 245 : SEGDES IMATRI,MATRIK 246 : RETURN 247 : 248 : END 249 : 250 : 251 : 252 : 253 : 254 : 255 : 256 : 257 : 258 :
© Cast3M 2003 - All rights reserved.
Disclaimer