1 : C THNUMAC SOURCE FANDEUR 11/07/19 21:17:52 7042 2 : 3 : C======================================================================= 4 : C= T H N U M A C = 5 : C= ----------- = 6 : C= (TNUMAC dans le cas de la thermique) = 7 : C= Fonction : = 8 : C= ---------- = 9 : C= Calcul de la matrice de CONDUCTIVITE THERMOHYDRIQUE pour les = 10 : C= elements finis MASSIFs a integration NUMERIQUE = 11 : C= = 12 : C= Parametres : (E)=Entree (S)=Sortie = 13 : C= ------------ = 14 : C= NEF (E) Numero de l'ELEMENT FINI dans NOMTP (cf. CCHAMP) = 15 : C= IMAIL (E) Numero du segment IMODEL dans le segment MMODEL = 16 : C= IPCHEM (E) Pointeur sur un segment MCHELM de CARACTERISTIQUES = 17 : C= IPRIGI (E/S) Pointeur sur l'objet RIGIDITE (CONDUCTIVITE) = 18 : C= = 19 : C= Zakaria HABIBI le 30 juin 2008. = 20 : C======================================================================= 21 : 22 : SUBROUTINE THNUMAC (NEF,ipmail,ipinte,ipint1,IVAMAT,NMATT, 23 : & ipmatr,LRE) 24 : 25 : IMPLICIT INTEGER(I-N) 26 : IMPLICIT REAL*8 (A-H,O-Z) 27 : 28 : -INC CCOPTIO 29 : -INC CCREEL 30 : 31 : -INC SMCHAML 32 : -INC SMCOORD 33 : -INC SMELEME 34 : -INC SMINTE 35 : -INC SMRIGID 36 : 37 : SEGMENT MMAT1 38 : REAL*8 VALMAT(NV1) 39 : REAL*8 CEL1(NBNN,NBNN),CEL2(NBNN,NBNN),CEL3(NBNN,NBNN) 40 : REAL*8 CEL4(NBNN,NBNN),CEL5(NBNN,NBNN),CEL6(NBNN,NBNN) 41 : REAL*8 CEL7(NBNN,NBNN),CEL8(NBNN,NBNN),CEL9(NBNN,NBNN) 42 : REAL*8 CE10(NBNN,NBNN),CE11(NBNN,NBNN),CE12(NBNN,NBNN) 43 : REAL*8 XE(3,NBNN),GDT1(IDIM),GDT2(IDIM),GDT3(IDIM) 44 : REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN),FORME(NBNN) 45 : REAL*8 CMAT(NDIM,NDIM),CMAT1(IDIM,IDIM),CMAT2(IDIM,IDIM) 46 : ENDSEGMENT 47 : 48 : SEGMENT MAXE 49 : REAL*8 TXR(IDIM,IDIM),XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM) 50 : ENDSEGMENT 51 : 52 : SEGMENT MPTVAL 53 : INTEGER IPOS(NS),NSOF(NS),IVAL(NCOSOU) 54 : CHARACTER*16 TYVAL(NCOSOU) 55 : ENDSEGMENT 56 : 57 : C 1 - INITIALISATIONS ET VERIFICATIONS 58 : C ====================================== 59 : C Recuperation d'informations sur le maillage elementaire 60 : C ===== 61 : MELEME = ipmail 62 : c* SEGACT,MELEME 63 : NBNN = NUM(/1) 64 : NBELEM = NUM(/2) 65 : C ===== 66 : C Recuperation d'informations sur l'element fini 67 : C ===== 68 : MINTE = IPINTE 69 : c* SEGACT,MINTE 70 : NBPGAU = POIGAU(/1) 71 : C ===== 72 : C Recuperation des fonctions de forme et de leurs derivees au 73 : C centre de gravite de l'element pour le calcul des axes locaux 74 : C d'orthotropie ou d'anisotropie 75 : C ===== 76 : IF (ipint1.NE.0) THEN 77 : MINTE1 = IPINT1 78 : c* SEGACT,MINTE1 79 : NBSH = MINTE1.SHPTOT(/2) 80 : ENDIF 81 : C ===== 82 : C Initialisation des segments de travail 83 : C ===== 84 : MPTVAL = IVAMAT 85 : IF (IFOMOD.EQ.1) THEN 86 : NDIM=3 87 : ELSE 88 : NDIM=IDIM 89 : ENDIF 90 : NV1 = NMATT 91 : SEGINI,MMAT1 92 : MAXE = 0 93 : IF (ipint1.GT.0) THEN 94 : SEGINI,MAXE 95 : ENDIF 96 : 97 : C ===== 98 : C Matrice CONDUCTIVITE GLOBALE a calculer 99 : C ===== 100 : XMATRI = ipmatr 101 : C* SEGACT,xMATRI*MOD 102 : NLIGRP = LRE 103 : NLIGRD = LRE 104 : 105 : C 2 - BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL/IMAMOD 106 : C ================================================================== 107 : DO iElt = 1, NBELEM 108 : C === 109 : C - Mise a zero de la matrice de CONDUCTIVITE de l'element iElt 110 : C === 110 : CALL ZERO(CEL1,NBNN,NBNN) 110 : CALL ZERO(CEL2,NBNN,NBNN) 110 : CALL ZERO(CEL3,NBNN,NBNN) 110 : CALL ZERO(CEL4,NBNN,NBNN) 110 : CALL ZERO(CEL5,NBNN,NBNN) 110 : CALL ZERO(CEL6,NBNN,NBNN) 110 : CALL ZERO(CEL7,NBNN,NBNN) 110 : CALL ZERO(CEL8,NBNN,NBNN) 110 : CALL ZERO(CEL9,NBNN,NBNN) 110 : CALL ZERO(CE10,NBNN,NBNN) 110 : CALL ZERO(CE11,NBNN,NBNN) 110 : CALL ZERO(CE12,NBNN,NBNN) 111 : C === 112 : C - Recuperation des coordonnees GLOBALES des noeuds de l'element 113 : C === 113 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,iElt,XE) 114 : C === 115 : C - Calculs des axes locaux d'orthotropie ou d'anisotropie 116 : C === 117 : IF (ipint1.NE.0) THEN 117 : CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR) 118 : IF (nbsh.EQ.-1) THEN 118 : CALL ERREUR(525) 119 : GOTO 9990 120 : ENDIF 121 : ENDIF 122 : C === 123 : C - Boucle sur les points de Gauss de l'element iElt 124 : C === 125 : iFois=0 126 : DO iGau = 1, NbPGau 127 : C ======= 128 : C 2.4.1 - Calcul du jacobien, des fonctions de forme et de leurs 129 : C derivees au point de Gauss iGau 130 : C ======= 130 : CALL TCOND5(iGau,NBNN,NDIM,XE,SHPTOT,SHP,GRAD,DJAC) 131 : IF (IERR.NE.0) GOTO 9990 132 : IF (DJAC.LT.XZERO) iFois=iFois+1 133 : DJAC=ABS(DJAC) 134 : C ======= 135 : C 2.4.2 - Erreur si le jacobien est nul en ce point de Gauss 136 : C ======= 137 : IF (DJAC.LT.XPetit) THEN 138 : INTERR(1)=iElt 138 : CALL ERREUR(259) 139 : GOTO 9990 140 : ENDIF 141 : 142 : C* VALEURS DES COMPOSANTES DES GRADIENTS POUR G 143 : C MELVAL=IVAL(19) 144 : C IBMN=MIN(iElt,VELCHE(/2)) 145 : C IGMN=MIN(iGau,VELCHE(/1)) 146 : C GDT1(1)=VELCHE(IGMN,IBMN) 147 : C MELVAL=IVAL(20) 148 : C IBMN=MIN(iElt,VELCHE(/2)) 149 : C IGMN=MIN(iGau,VELCHE(/1)) 150 : C GDT1(2)=VELCHE(IGMN,IBMN) 151 : C* VALEURS DES COMPOSANTES DES GRADIENTS POUR C 152 : C MELVAL=IVAL(19) 153 : C IBMN=MIN(iElt,VELCHE(/2)) 154 : C IGMN=MIN(iGau,VELCHE(/1)) 155 : C GDT2(1)=VELCHE(IGMN,IBMN) 156 : C MELVAL=IVAL(20) 157 : C IBMN=MIN(iElt,VELCHE(/2)) 158 : C IGMN=MIN(iGau,VELCHE(/1)) 159 : C GDT2(2)=VELCHE(IGMN,IBMN) 160 : C* VALEURS DES COMPOSANTES DES GRADIENTS POUR T 161 : C MELVAL=IVAL(19) 162 : C IBMN=MIN(iElt,VELCHE(/2)) 163 : C IGMN=MIN(iGau,VELCHE(/1)) 164 : C GDT3(1)=VELCHE(IGMN,IBMN) 165 : C MELVAL=IVAL(20) 166 : C IBMN=MIN(iElt,VELCHE(/2)) 167 : C IGMN=MIN(iGau,VELCHE(/1)) 168 : C GDT3(2)=VELCHE(IGMN,IBMN) 169 : 170 : C ======= 171 : C 2.4.3 - Recuperation de la ou des valeurs de conductibilite au point 172 : C de Gauss iGau (tableau VALMAT) 173 : C ======= 174 : DJAC=DJAC*POIGAU(iGau) 175 : DO i = 1, NMATT 176 : IF (IVAL(i).NE.0) THEN 177 : MELVAL=IVAL(i) 178 : IBMN=MIN(iElt,VELCHE(/2)) 179 : IGMN=MIN(iGau,VELCHE(/1)) 180 : VALMAT(i)=VELCHE(IGMN,IBMN) 181 : ELSE 182 : VALMAT(i)=XZERO 183 : ENDIF 184 : ENDDO 185 : C ======= 186 : C 2.4.4 - Cas d'un materiau ISOTROPE de conductibilite K 187 : C Calcul de la contribution du point de Gauss a la matrice 188 : C CONDUCTIVITE elementaire de cet element fini 189 : C Ajout du terme XK*transposee(GRAD)*GRAD 190 : C Seul cas valide en dimension 1 191 : C ======= 192 : XK=VALMAT(1)*DJAC 192 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL1) 193 : XK=VALMAT(2)*DJAC 193 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL2) 194 : XK=VALMAT(3)*DJAC 194 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL3) 195 : XK=VALMAT(4)*DJAC 195 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL4) 196 : XK=VALMAT(5)*DJAC 196 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL5) 197 : XK=VALMAT(6)*DJAC 197 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL6) 198 : XK=VALMAT(7)*DJAC 198 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL7) 199 : XK=VALMAT(8)*DJAC 199 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL8) 200 : XK=VALMAT(9)*DJAC 200 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL9) 201 : 202 : CC Calcul des termes en GRAD(PG,PC,T) 203 : C RK =DJAC 204 : C DO 700 K=1,IDIM 205 : C XK1 = GDT1(K)*RK 206 : C DO 300 I=1,NBNN 207 : C DO 400 J=1,NBNN 208 : C CE10(J,I) = CE10(J,I) + SHP(1,I)*GRAD(K,J)*XK1 209 : C 400 CONTINUE 210 : C 300 CONTINUE 211 : C 700 CONTINUE 212 : C 213 : C DO 701 K=1,IDIM 214 : C XK2 = GDT2(K)*RK 215 : C DO 301 I=1,NBNN 216 : C DO 401 J=1,NBNN 217 : C CE11(J,I) = CE11(J,I) + SHP(1,I)*GRAD(K,J)*XK2 218 : C 401 CONTINUE 219 : C 301 CONTINUE 220 : C 701 CONTINUE 221 : C 222 : C 223 : C DO 702 K=1,IDIM 224 : C XK3 = GDT3(K)*RK 225 : C DO 302 I=1,NBNN 226 : C DO 402 J=1,NBNN 227 : C CE12(J,I) = CE12(J,I) + SHP(1,I)*GRAD(K,J)*XK3 228 : C 402 CONTINUE 229 : C 302 CONTINUE 230 : C 702 CONTINUE 231 : ENDDO 232 : 233 : C ===== 234 : C 2.5 - Erreur si, en un point de Gauss, le jacobien change de signe. 235 : C ===== 236 : IF (iFois.NE.0.AND.iFois.NE.NbPGau) THEN 237 : INTERR(1)=iElt 237 : CALL ERREUR(195) 238 : GOTO 9990 239 : ENDIF 240 : C ===== 241 : C 2.6 - Stockage de la matrice de CONDUCTIVITE pour cet element fini 242 : C ===== (remplissage de XMATRI) 243 : C Note : CE10,CE11,CE12 servaient a calculer les termes en Grad(PG,PC,T) 244 : C du couplage. Desormais laisses a 0. en attendant de voir com- 245 : C -ment traiter le couplage (champ (PG,PC,T) en entree de COND ?) 245 : CALL REMPTH 246 : & (CEL1,CEL2,CEL3,CEL4,CEL5,CEL6, 247 : & CEL7,CEL8,CEL9,CE10,CE11,CE12,NBNN,RE(1,1,ielt)) 248 : 249 : ENDDO 250 : 251 : C 3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS 252 : C ==================================================== 253 : 9990 CONTINUE 254 : SEGSUP,MMAT1 255 : IF (IPINT1.GT.0) THEN 256 : SEGSUP,MAXE 257 : ENDIF 258 : 259 : RETURN 260 : END 261 : 262 :
© Cast3M 2003 - All rights reserved.
Disclaimer