1 : C MHYBR2 SOURCE CHAT 11/03/16 21:27:56 6902 2 : SUBROUTINE MHYBR2(IMAIL,IPMODE,IPCHEM,IPRIGI,IPGEOS,ILUMP) 3 : C----------------------------------------------------------------------- 4 : C Calcul de l'inverse de la matrice masse hybride. 5 : C Traitement du cas des elements finis massifs a integration numerique 6 : C pour un maillage elementaire et une formulation hybride. 7 : C----------------------------------------------------------------------- 8 : C 9 : C--------------------------- 10 : C Parametres Entree/Sortie : 11 : C--------------------------- 12 : C 13 : C E/ IMAIL : Numero du maillage elementaire considere, 14 : C dans l'objet modele. 15 : C E/ IPMODE : Pointeur sur un segment IMODEL. 16 : C E/ IPCHEM : Pointeur sur un chamelem de caracteristiques. 17 : C E/ IPGEOS : Pointeur sur le maillage sommet 18 : C E/S IPRIGI : Pointeur sur l'objet rigidite resultat. 19 : C 20 : C---------------------- 21 : C Variables en COMMON : 22 : C---------------------- 23 : C 24 : C E/ XCOOR : VOIR SMCOORD 25 : C E/ IERR : VOIR CCOPTIO 26 : C E/ IDIM : VOIR CCOPTIO 27 : C E/ INTERR : VOIR CCOPTIO 28 : C E/ IFOMOD : VOIR CCOPTIO 29 : C E/ XPETIT : VOIR CCREEL 30 : C 31 : C---------------------- 32 : C Tableaux de travail : 33 : C---------------------- 34 : C 35 : C NBNN : Nombre de noeuds dans l'element considere 36 : C NEFHYB : Numero de l'element fini dans NOMTP. 37 : C NEF : Numero de l'element fini support geometrique 38 : C dans NOMTP (voir CCHAMP) 39 : C NBELEM : Nombre d'element dans le maillage elementaire 40 : C NBPGAU : Nombre de points de gauss pour l'element fini NEF 41 : C CEL : Matrice de conductivite elementaire 42 : C XE : Coordonnees des noeuds dans le repere global 43 : C CMAT : Matrice de permeabilite dans le repere global 44 : C SHP : Tableau de travail contenant les fonctions de forme au 45 : C point de gauss utilise + derivees 46 : C SHY : Contient les fonctions de base hybride en un point 47 : C mais pas les derivees de la fonction de base. 48 : C VALMAT : Valeurs des coeff. de la matrice CMAT et des 49 : C cosinus directeurs des axes d'ortho. / repere local 50 : C XGLOB : Cosinus directeurs des axes d'ortho. / repere global 51 : C XLOC : Cosinus directeurs des axes d'ortho. / repere local 52 : C TXR : Cosinus directeurs des axes locaux / repere global 53 : C 54 : C 55 : C----------------------------------------------------------------------- 56 : C 57 : C Langage : ESOPE + FORTRAN77 58 : C 59 : C Auteurs : F.DABBENE 08/93 60 : C 61 : C----------------------------------------------------------------------- 62 : IMPLICIT INTEGER(I-N) 63 : IMPLICIT REAL*8(A-H,O-Z) 64 : * 65 : -INC CCHAMP 66 : -INC CCOPTIO 67 : -INC CCREEL 68 : -INC SMCOORD 69 : -INC SMINTE 70 : -INC SMMODEL 71 : -INC SMRIGID 72 : -INC SMELEME 73 : -INC SMCHAML 74 : * 75 : SEGMENT MAXE 76 : REAL*8 TXR(IDIM,IDIM),XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM) 77 : ENDSEGMENT 78 : * 79 : SEGMENT MMAT1 80 : REAL*8 CEL(NBDDL,NBDDL),CEL1(NBDDL,NBDDL),XE(3,NBNN) 81 : REAL*8 SHP(6,NBNN),SHY(IDIM,NBDDL),ZJAC(IDIM,IDIM) 82 : REAL*8 CMAT(IDIM,IDIM),CMAT1(IDIM,IDIM),CMAT2(IDIM,IDIM) 83 : INTEGER ICSTO(NBDDL) 84 : ENDSEGMENT 85 : * 86 : SEGMENT NOTYPE 87 : CHARACTER*16 TYPE(NBTYPE) 88 : ENDSEGMENT 89 : * 90 : SEGMENT MVELCH 91 : REAL*8 VALMAT(NV1) 92 : ENDSEGMENT 93 : * 94 : SEGMENT MPTVAL 95 : INTEGER IPOS(NS) , NSOF(NS) 96 : INTEGER IVAL(NCOSOU) 97 : CHARACTER*16 TYVAL(NCOSOU) 98 : ENDSEGMENT 99 : * 100 : SEGMENT HYBSTO 101 : REAL*8 HYBASE(NDIM,NBDDL,NBPP) 102 : ENDSEGMENT 103 : * 104 : CHARACTER*8 CNM 105 : CHARACTER*(NCONCH) CONM 106 : PARAMETER(NINF=3) 107 : INTEGER INFOS(NINF) 108 : LOGICAL lsupma 109 : * 110 : * Recup. des caracteristiques geometriques du maillage elementaire 111 : * et du maillage hybride dual 112 : * 113 : IMODEL = IPMODE 114 : SEGACT IMODEL 115 : CONM = CONMOD 116 : IPMAIL = IMAMOD 117 : MELEME = IPGEOS 118 : SEGACT MELEME 119 : NBNN = NUM(/1) 120 : NBELEM = NUM(/2) 121 : NEFHYB = NEFMOD 122 : NEF = NUMGEO(NEFHYB) 123 : MFR = NUMMFR(NEFHYB) 124 : * 125 : MRIGID = IPRIGI 126 : SEGACT MRIGID 127 : IPT1 = IRIGEL(1,IMAIL) 128 : SEGACT IPT1 129 : NBDDL = IPT1.NUM(/1) 130 : SEGDES IPT1 131 : * 132 : * Recup. des caracteristiques d'integration de l'EF support geometrique 133 : * de l'EF hybride 134 : * 134 : CALL TSHAPE(NEF,'GAUSS',IPINTE) 135 : IF (IERR.NE.0) THEN 136 : SEGDES IMODEL , MELEME 137 : RETURN 138 : ENDIF 139 : * 140 : * Recup. des fonctions de bases hybrides 141 : * 141 : CALL HYSHPT(NEFHYB,NBDDL,IPINTE,IPTHYB) 142 : IF (IERR.NE.0) THEN 143 : SEGDES IMODEL , MELEME 144 : RETURN 145 : ENDIF 146 : * 147 : * Activation des segments "d'integration" 148 : * 149 : MINTE = IPINTE 150 : SEGACT MINTE 151 : NBPGAU = POIGAU(/1) 152 : HYBSTO = IPTHYB 153 : SEGACT HYBSTO 154 : * 155 : * Recup. des caracteristiques d'integration au centre de l'EF 156 : * 156 : CALL RESHPT(1,NBNN,NEF,NEF,0,IPT1,IRT1) 157 : MINTE1 = IPT1 158 : SEGACT MINTE1 159 : * 160 : * Initialisation des chapeaux de l'objet rigidité 161 : * 162 : xMATRI = IRIGEL(4,IMAIL) 163 : SEGACT xMATRI*MOD 164 : NLIGRP = NBDDL 165 : NLIGRD = NBDDL 166 : * 167 : * Recherche du nom du materiau ...trope, de son numero, de sa nature 168 : * 169 : NFOR = FORMOD(/2) 170 : NMAT = MATMOD(/2) 170 : CALL NOMATE(FORMOD,NFOR,MATMOD,NMAT,CNM,INM,INT) 171 : IF (IERR.NE.0) THEN 172 : SEGDES IMODEL , MELEME 173 : SEGDES MINTE , MINTE1 174 : SEGSUP xMATRI , MRIGID , HYBSTO 175 : RETURN 176 : ENDIF 177 : * 178 : * Remplissage du tableau INFOS (informations sur element). 179 : * 180 : INFOS(1) = 0 181 : INFOS(2) = 0 182 : INFOS(3) = NIFOUR 183 : * 184 : * Remplissage de MOMATR : Noms des composantes obligatoires et 185 : * facultatives contenus dans IPCHEM. 186 : * 187 : if(lnomid(6).ne.0) then 188 : nomid=lnomid(6) 189 : segact nomid 190 : momatr=nomid 191 : nmatr=lesobl(/2) 192 : nmatf=lesfac(/2) 193 : lsupma=.false. 194 : else 195 : lsupma=.true. 195 : CALL IDMATR(MFR,IMODEL,MOMATR,NMATR,NMATF) 196 : endif 197 : NMATT = NMATR 198 : NV1 = NMATT 199 : * 200 : * Initialisation de NOTYPE : Type des infos contenus dans IPCHEM. 201 : * 202 : NBTYPE = 1 203 : SEGINI NOTYPE 204 : MOTYPE = NOTYPE 205 : TYPE(1) ='REAL*8' 206 : * 207 : * Verification des informations transmises et recuperation de IVAMAT, 208 : * pointeur vers le segment MPTVAL contenant les pointeurs vers les 209 : * composantes obligatoires du MCHAML de caracteristiques. 210 : * 210 : CALL KOMCHA(IPCHEM,IPGEOS,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT) 211 : SEGSUP NOTYPE 212 : SEGACT MELEME 213 : IF (IERR.NE.0) THEN 214 : SEGDES MELEME 215 : SEGDES IMODEL 216 : SEGDES MINTE , MINTE1 217 : SEGSUP xMATRI , MRIGID , HYBSTO 218 : RETURN 219 : ENDIF 220 : * 221 : * Initialisation des tableaux de travail 222 : * 223 : NDIM = IDIM * (IDIM+1) 224 : SEGINI MMAT1 , MVELCH , MAXE 225 : * 226 : *------------------------------------------------------- 227 : * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL 228 : *------------------------------------------------------- 229 : * 230 : DO 30 IEL=1,NBELEM 231 : * 232 : * Initialisations 233 : * 234 : IFOIS = 0 235 : IFOI2 = 0 235 : CALL ZERO(CEL,NBDDL,NBDDL) 236 : * 237 : * Recuperation des coordonnees globales des noeuds de l'element IEL 238 : * 238 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE) 239 : * 240 : * Calcul des axes locaux dans le cas orthotrope ou anisotrope 241 : * 242 : IF (INM.EQ.2.OR.INM.EQ.3) THEN 243 : NBSH = MINTE1.SHPTOT(/2) 243 : CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR) 244 : if (nbsh.eq.-1) then 245 : 245 : call erreur(525) 246 : return 247 : endif 248 : ENDIF 249 : 250 : *-------------------------------- 251 : * BOUCLE SUR LES POINTS DE GAUSS 252 : *-------------------------------- 253 : DO 20 IGAU=1,NBPGAU 254 : * 255 : * Calcul de la matrice jacobienne de la fonction de passage du 256 : * repere local au repere global, de son determinant ET recup. 257 : * des fonctions de base hybride au point de gauss IGAU. 258 : * 258 : CALL MHYBR3(IGAU,NBNN,NBDDL,NDIM,IDIM,IDIM,XE,HYBASE, 259 : S SHPTOT,SHY,SHP,ZJAC,DJAC) 260 : * 261 : * Controle du maillage 262 : * 263 : IF (DJAC.LT.0.D0) IFOIS = IFOIS + 1 264 : IF (ABS(DJAC).LT.XPETIT) THEN 265 : IFOI2 = IFOI2 + 1 266 : DJAC = XPETIT 267 : ENDIF 268 : * 269 : * Calcul du poids d'integration global affecte dans DJAC. 270 : * 271 : DJAC = POIGAU(IGAU) / ABS(DJAC) 272 : * 273 : * Recherche des valeurs des composantes de la matrice de permeabilite : 274 : * VALMAT(im) contient la im° composante au point de gauss IGAU. 275 : * 276 : MPTVAL = IVAMAT 277 : DO 10 IM=1,NMATT 278 : IF (IVAL(IM).NE.0) THEN 279 : MELVAL = IVAL(IM) 280 : IBMN = MIN(IEL,VELCHE(/2)) 281 : IGMN = MIN(IGAU,VELCHE(/1)) 282 : VALMAT(IM) = VELCHE(IGMN,IBMN) 283 : ELSE 284 : VALMAT(IM) = 0.D0 285 : ENDIF 286 : 10 CONTINUE 287 : * 288 : *= Passage de la matrice de permeabilite du repere local au global 289 : * 289 : CALL MATGLO(CMAT,CMAT1,CMAT2,TXR,XLOC,XGLOB,VALMAT, 290 : S IDIM,IDIM,INM,IFOMOD) 291 : * 292 : *- Calcul de la contribution du point de gauss IGAU a la matrice 293 : *- elementaire CEL de l'element IEL : 294 : *- POIGAU/DJAC * transpose( ZJAC*SHY ) *inv(CMAT)* ( ZJAC*SHY ) 295 : *- On calcule CMAT2=inv(CMAT) avec INVRS puis 296 : *- on calcule CMAT1=transpose(ZJAC)*CMAT2*ZJAC avec PRODT puis 297 : *- on somme POIGAU/DJAC * transp.(SHY)*CMAT1*SHY avec BDBST. 298 : * 298 : CALL INVRS(CMAT,IDIM,CMAT2,CJAC) 299 : IF (CJAC.EQ.0.D0) THEN 299 : CALL ERREUR(406) 300 : INTERR(1) = IEL 300 : CALL ERREUR(259) 301 : ENDIF 302 : IF (IERR.NE.0) THEN 303 : SEGSUP xMATRI , MRIGID 304 : GOTO 40 305 : ENDIF 305 : CALL PRODT(CMAT1,CMAT2,ZJAC,IDIM,IDIM) 305 : CALL BDBST(SHY,DJAC,CMAT1,NBDDL,IDIM,CEL) 306 : 20 CONTINUE 307 : * 308 : * Le jacobien est negatif --> MAILLAGE INCORRECT 309 : * 310 : IF (IFOIS.NE.0.AND.IFOIS.NE.NBPGAU) THEN 311 : INTERR(1) = IEL 311 : CALL ERREUR(195) 312 : SEGSUP xMATRI , MRIGID 313 : GOTO 40 314 : ENDIF 315 : * 316 : * Le jacobien est tres petit --> MAILLAGE INCORRECT 317 : * 318 : IF (IFOI2.EQ.NBPGAU) THEN 319 : INTERR(1) = IEL 319 : CALL ERREUR(259) 320 : SEGSUP xMATRI , MRIGID 321 : GOTO 40 322 : ENDIF 323 : * 324 : * Lump pour les carrés 325 : IF(ILUMP.GT.0)THEN 326 : * 327 : DO 60 II= 1,NBDDL 328 : SCLI=0.D0 329 : DO 65 JJ=1,NBDDL 330 : SCLI= ABS(CEL(II,JJ))+SCLI 331 : CEL(II,JJ)=0.D0 332 : 65 CONTINUE 333 : CEL(II,II)=SCLI 334 : 60 CONTINUE 335 : ENDIF 336 : * 337 : * Inversion de la matrice masse hybride 338 : * 338 : CALL INVER(CEL,NBDDL,ICRIT,CEL1,ICSTO,XPETIT) 339 : IF (ICRIT.EQ.1) THEN 340 : MOTERR(1:8) = 'DARCY' 340 : CALL ERREUR(700) 341 : SEGSUP xMATRI , MRIGID 342 : GOTO 40 343 : ENDIF 344 : * 345 : * Remplissage de XMATRI 346 : * 347 : * SEGINI XMATRI 348 : * IMATTT(IEL) = XMATRI 348 : CALL REMPMT(CEL,NBDDL,RE(1,1,iel)) 349 : * SEGDES XMATRI 350 : 30 CONTINUE 351 : * 352 : * Desactivation des segments 353 : * 354 : SEGDES xMATRI , MRIGID 355 : 40 CONTINUE 356 : SEGSUP MMAT1 , MVELCH , MAXE , HYBSTO 357 : SEGDES MELEME 358 : SEGDES IMODEL 359 : SEGDES MINTE , MINTE1 360 : * 361 : MPTVAL = IVAMAT 362 : DO 50 I=1,NMATT 363 : MELVAL = IVAL(I) 364 : SEGDES MELVAL 365 : 50 CONTINUE 366 : SEGSUP MPTVAL 367 : NOMID = MOMATR 368 : if(lsupma)SEGSUP NOMID 369 : * 370 : RETURN 371 : END 372 : 373 : 374 : 375 :
© Cast3M 2003 - All rights reserved.
Disclaimer