1 : C ROT3M SOURCE CHAT 11/03/16 21:31:22 6902 2 : SUBROUTINE ROT3M(NEF,IMAIL,IPMODE,IPCHEM,IPSUPJ,XMATRI, 3 : $ ICPR,ICPR2) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8(A-H,O-Z) 6 : ************************************************************************ 7 : * 8 : * R O T 3 M 9 : * --------- 10 : * 11 : * FONCTION: 12 : * --------- 13 : * CALCUL DE LA MATRICE DE MUTUELLES POUR L'ELEMENT ROT3 14 : * 15 : * MODULES UTILISES: 16 : * ----------------- 17 : * 18 : -INC CCOPTIO 19 : -INC CCREEL 20 : *- 21 : -INC SMCOORD 22 : -INC SMINTE 23 : -INC CCHAMP 24 : -INC SMMODEL 25 : -INC SMRIGID 26 : -INC SMELEME 27 : -INC SMCHAML 28 : * 29 : * PARAMETRES: (E)=ENTREE (S)=SORTIE (+ = CONTENU DANS UN COMMUN) 30 : * ----------- 31 : * 32 : * NEF (E) NUMERO DE L'ELEMENT-FINI DANS NOMTP (VOIR CCHAMP) 33 : * IMAIL (E) NUMERO DU MAILLAGE ELEMENTAIRE CONSIDERE,DANS 34 : * L'OBJET MODELE 35 : * IPMODE (E) POINTEUR SUR UN SEGMENT IMODEL 36 : * IPCHEM (E) POINTEUR SUR LE CHAMELEM DE CARACTERISTIQUE 37 : * IPSUPJ (E) POINTEUR SUR LE MAILLAGE SUPPORT DES COURANTS DE FOUCAULT 38 : * +XCOOR (E) VOIR SMCOORD 39 : * +IDIM (E) VOIR CCOPTIO 40 : * +IFOMOD (E) VOIR CCOPTIO 41 : * +XZERO (E) VOIR CCREEL 42 : * IPRIGI (E/S) POINTEUR SUR L'OBJET RESULTAT,DE TYPE RIGIDITE 43 : * 44 : * VARIABLES: 45 : * ---------- 46 : * 47 : * NBNN NOMBRE DE NOEUDS DANS L'ELEMENT CONSIDERE 48 : * NEF NUMERO DE L'ELEMENT FINI DANS NOMTP (VOIR CCHAMP) 49 : * NBELEM NOMBRE D'ELEMENTS DANS LE MAILLAGE ELEMENTAIRE 50 : * NBPGAU NOMBRE DE POINTS DE GAUSS DANS L'ELEMENT-FINI 51 : * NDIM NOMBRE DE LIGNES DE LA MATRICE GRADIENT 52 : * XE(3,NBNN) COORDONNEES DE L'ELEMENT DANS LE REPERE GLOBAL 53 : * SHP(6,NBNN) TABLEAU DE TRAVAIL 54 : * GRAD(NDIM,NBNN) MATRICE GRADIENT DES FONCTIONS DE FORME BIDIM. 55 : * VALMAT(NMATR) TABLEAU DE TRAVAIL 56 : * 57 : SEGMENT MAXT 58 : REAL*8 RA(NBNN,NBNN) 59 : ENDSEGMENT 60 : SEGMENT,MMAT1 61 : REAL*8 VALMAT(NMATR) 62 : REAL*8 XE(3,NBNN),XE1(3,NBNN) 63 : REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN,NBPGAU) 64 : REAL*8 COSD1(3),COSD2(3) 65 : ENDSEGMENT 66 : POINTEUR MMAT2.MMAT1,MMATX.MMAT1 67 : * 68 : SEGMENT NOTYPE 69 : CHARACTER*16 TYPE(NBTYPE) 70 : ENDSEGMENT 71 : * 72 : SEGMENT MPTVAL 73 : INTEGER IPOS(NS) ,NSOF(NS) 74 : INTEGER IVAL(NCOSOU) 75 : CHARACTER*16 TYVAL(NCOSOU) 76 : ENDSEGMENT 77 : * 78 : SEGMENT INFO 79 : INTEGER INFELL(JG) 80 : ENDSEGMENT 81 : * 82 : SEGMENT SGAUSS 83 : REAL*8 XGAUSS(3,NBPGAU) 84 : REAL*8 DX(NBPGAU) 85 : ENDSEGMENT 86 : POINTEUR SGX.SGAUSS,SGY.SGAUSS 87 : SEGMENT ICPR(NA) 88 : SEGMENT ICPR2(NA) 89 : * 90 : CHARACTER*8 CNM 91 : CHARACTER*(NCONCH) CONM 92 : PARAMETER (NINF=3) 93 : INTEGER INFOS(NINF) 94 : LOGICAL SELF,NEAR 95 : * 96 : * 97 : * 98 : * AUTEUR, DATE DE CREATION: 99 : * ------------------------- 100 : * 101 : * YANN STEPHAN , FEVRIER 1997 (COPIE DE ROT3R) 102 : * 103 : * LANGAGE: 104 : * -------- 105 : * 106 : * ESOPE + FORTRAN77 107 : * 108 : ************************************************************************ 109 : * PERMEABILITE DU VIDE SUR 4PI 110 : DATA PM0S4P/1.D-7/ 111 : * 112 : * RECUPERATION DES CARACTERISTIQUES GEOMETRIQUES DU MAILLAGE 113 : * ELEMENTAIRE 114 : * 115 : IMODEL=IPMODE 116 : CONM =CONMOD 117 : IPMAIL=IMAMOD 118 : MELEME=IMAMOD 119 : SEGACT,MELEME 120 : NBNN=NUM(/1) 121 : NBELEM=NUM(/2) 122 : * 123 : * RECUPERATION DES CARACTERISTIQUES D'INTEGRATION 124 : * POUR LA MATRICE MUTUELLE (RIGIDITE) DE L'ELEMENT 125 : * FINI LIE A NOTRE MAILLAGE 126 : * 127 : if(infmod(/1).lt.5) then 127 : CALL ELQUOI(NEF,0,3,IPINF,IMODEL) 128 : IF(IERR.NE.0) THEN 129 : SEGDES IMODEL*NOMOD,MELEME*NOMOD 130 : RETURN 131 : ENDIF 132 : INFO=IPINF 133 : IPINTE=INFELL(11) 134 : segsup info 135 : else 136 : ipinte=infmod(5) 137 : endif 138 : * 139 : * INFORMATION SUR L'ELEMENT 140 : * 141 : MINTE=IPINTE 142 : SEGACT,MINTE 143 : NBPGAU=POIGAU(/1) 144 : SEGINI SGX,SGY 145 : * 146 : * RECHERCHE LES POINTEURS DES SEGMENTS MELVAL QUI CORRESPONDENT 147 : * A LA PERMEABILITE ET L'EPAISSEUR DES ELEMENT 148 : * 149 : NFOR=FORMOD(/2) 150 : NMAT=MATMOD(/2) 150 : CALL NOMATE(FORMOD,NFOR,MATMOD,NMAT,CNM,INM,INT) 151 : IF (IERR.NE.0) THEN 152 : SEGDES,MELEME*NOMOD 153 : SEGDES,MINTE*NOMOD 154 : RETURN 155 : ENDIF 156 : * 157 : * REMLIR LE TABLEAU INFOS (INFORMATIONS SUR ELEMENT) 158 : INFOS(1)=0 159 : INFOS(2)=0 160 : INFOS(3)=NIFOUR 161 : * 162 : IF(CNM.EQ.'ISOTROPE'.OR.CNM.EQ.'ORTHOTRO')THEN 163 : NBRFAC=0 164 : NBROBL=1 165 : SEGINI NOMID 166 : MOMATR=NOMID 167 : LESOBL(1)='PERM' 168 : NMATR=1 169 : NMATF=0 170 : ELSE 170 : CALL ERREUR(251) 171 : SEGDES, MELEME*NOMOD 172 : SEGDES, MINTE*NOMOD 173 : RETURN 174 : ENDIF 175 : SEGDES NOMID 176 : * 177 : NBTYPE=1 178 : SEGINI NOTYPE 179 : MOTYPE=NOTYPE 180 : TYPE(1)='REAL*8' 181 : * 181 : CALL KOMCHA(IPCHEM,IPMAIL,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT) 182 : SEGSUP NOTYPE 183 : IF(IERR.NE.0)THEN 184 : SEGDES,MELEME*NOMOD 185 : SEGDES,MINTE*NOMOD 186 : RETURN 187 : ENDIF 188 : * 189 : MPTVAL=IVAMAT 190 : * 191 : NDIM=IDIM 192 : NDIM2=IDIM-1 193 : SEGINI MAXT 194 : SEGINI,MMAT1 195 : MMATX=MMAT1 196 : SEGINI,MMAT2 197 : * 198 : * CALCUL DE LA DISTANCE DE REFERENCE 199 : * 200 : DREF=0. 201 : MIJ=0 202 : DO 20 IEL=1,NBELEM 202 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE) 202 : CALL MAXLT3(NBNN,XE,DARET) 203 : DREF=MAX(DREF,DARET) 204 : 20 CONTINUE 205 : * 206 : * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL 207 : 208 : * 209 : SEGACT ICPR,ICPR2 210 : DO 10 IEL=1,NBELEM 211 : 212 : * 213 : MMAT1=MMATX 214 : SGAUSS=SGX 215 : * 216 : * ON CHERCHE LES COORDONNEES DES NOEUDS DE L'ELEMENT IEL, 217 : * DANS LE REPERE GLOBAL 218 : * 218 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE) 219 : * 220 : * CALCUL DES COORDONNEES DES NOEUDS DANS LE REPERE LOCAL DE L' 221 : * ELEMENT COQUE 222 : * 222 : CALL COQLOC(NBNN,XE,COSD1,COSD2,XE1) 223 : * 224 : * ON CALCULE LES FONCTIONS DE FORME AUX POINTS DE GAUSS 225 : * 225 : CALL ELGAUS(MINTE,MMAT1,SGAUSS,IFOIS,IFOI2) 226 : * 227 : IF(IFOIS.NE.0.AND.IFOIS.NE.NBPGAU)THEN 228 : * 229 : * LE JACOBIEN EST NEGATIF ,MAILLAGE INCORRECT 230 : INTERR(1)=IEL 230 : CALL ERREUR(195) 231 : GO TO 999 232 : ELSEIF(IFOI2.EQ.NBPGAU)THEN 233 : * 234 : * CAS OU LE JACOBIEN EST TRES PETIT 235 : * 236 : INTERR(1)=IEL 236 : CALL ERREUR (259) 237 : GO TO 999 238 : ENDIF 239 : * 240 : * ON BOUCLE SUR LE MAILLAGE SUPPORT DE COURANTS 241 : IPT1=IPSUPJ 242 : SEGACT, IPT1 243 : NBSOUJ=IPT1.LISOUS(/1) 244 : IF(NBSOUJ.EQ.0) NBSOUJ=1 245 : DO 110 ISOUJ=1,NBSOUJ 246 : IF(NBSOUJ.EQ.1) THEN 247 : IPT2=IPT1 248 : ELSE 249 : IPT2=IPT1.LISOUS(ISOUJ) 250 : SEGACT, IPT2 251 : ENDIF 252 : NBELJ=IPT2.NUM(/2) 253 : NBNNJ=IPT2.NUM(/1) 254 : NBNNT=NBNN+NBNNJ 255 : NLIGRP=NBNN 256 : NLIGRD=NBNN 257 : DO 111 IELJ=1,NBELJ 258 : DO 230 IX=1,NBNN 259 : DO 230 JX=1,NBNN 260 : 230 RA(JX,IX) = 0.D0 261 : * 262 : NEAR=.FALSE. 263 : * 264 : * 265 : MMAT1=MMAT2 266 : SGAUSS=SGY 267 : * 268 : * ON CHERCHE LES COORDONNEES DES NOEUDS DE L'ELEMENT IEL, 269 : * DANS LE REPERE GLOBAL 270 : * 270 : CALL DOXE(XCOOR,IDIM,NBNN,IPT2.NUM,IELJ,XE) 271 : * 272 : * CALCUL DES COORDONNEES DES NOEUDS DANS LE REPERE LOCAL DE L' 273 : * ELEMENT COQUE 274 : * 274 : CALL COQLOC(NBNN,XE,COSD1,COSD2,XE1) 275 : * 275 : CALL ELGAUS(MINTE,MMAT1,SGAUSS,JFOIS,JFOI2) 276 : * 277 : IF(JFOIS.NE.0.AND.JFOIS.NE.NBPGAU)THEN 278 : * 279 : * LE JACOBIEN EST NEGATIF ,MAILLAGE INCORRECT 280 : INTERR(1)=IEL 280 : CALL ERREUR(195) 281 : GO TO 999 282 : ELSEIF(JFOI2.EQ.NBPGAU)THEN 283 : * 284 : * CAS OU LE JACOBIEN EST TRES PETIT 285 : * 286 : INTERR(1)=IEL 286 : CALL ERREUR (259) 287 : GO TO 999 288 : ENDIF 289 : * 290 : * CALCUL DE LA DISTANCE ENTRE LES DEUX ELEMENTS 291 : * 291 : CALL S2ELT3(NBNN,NUM,IEL,IPT2.NUM,IELJ,NEAR,SELF) 291 : CALL D2ELT3(NBNN,XE,MMATX.XE,DT3) 292 : NEAR=NEAR.OR.(DT3.LE.DREF) 293 : * 294 : * BOUCLE SUR LES POINTS DE GAUSS MAILLAGE 1 295 : * 296 : DO 22 IGAU=1,NBPGAU 297 : * 298 : MMAT1=MMATX 299 : SGAUSS=SGX 300 : * 301 : * ON CHERCHE LES VALEURS DE LA PERMEABILITE 302 : * 303 : MPTVAL=IVAMAT 304 : DO 30 IM=1,NMATR 305 : IF(IVAL(IM).NE.0)THEN 306 : MELVAL=IVAL(IM) 307 : IBMN=MIN(IEL,VELCHE(/2)) 308 : IGMN=MIN(IGAU,VELCHE(/1)) 309 : VALMAT(IM)=VELCHE(IGMN,IBMN) 310 : ELSE 311 : VALMAT(IM)=0 312 : ENDIF 313 : 30 CONTINUE 314 : PERM=VALMAT(1) 315 : XK=PERM*PM0S4P*DX(IGAU) 316 : * 317 : * BOUCLE SUR LES POINTS DE GAUSS MAILLAGE 2 318 : * 319 : DO 23 JGAU=1,NBPGAU 320 : * 321 : MMAT1=MMAT2 322 : SGAUSS=SGY 323 : * 324 : IF(SELF) THEN 325 : IF(JGAU.GT.1) GO TO 23 325 : CALL SELFT3(SGX.XGAUSS(1,IGAU),NBNN,XE,QQ) 326 : YK=XK*QQ 327 : ELSE IF(NEAR) THEN 328 : IF(JGAU.GT.1) GO TO 23 328 : CALL NEART3(SGX.XGAUSS(1,IGAU),NBNN,MMATX.XE,XE,QQ) 329 : YK=XK*QQ 330 : ELSE 331 : DIST=0. 332 : DO 120 I=1,IDIM 333 : DIST=DIST+(SGX.XGAUSS(I,IGAU)-XGAUSS(I,JGAU))**2 334 : 120 CONTINUE 335 : DIST=SQRT(DIST) 336 : YK=XK*DX(JGAU)/DIST 337 : ENDIF 338 : * 339 : * ON AJOUTE LE PRODUIT K*DJAC*TRANSPOSEE(GRADX)*GRADY 340 : * POUR LE POINT DE GAUSS CONSIDERE,A LA MATRICE RE 341 : * 342 : MMAT1=MMATX 342 : CALL WXWYSR(GRAD(1,1,IGAU),MMAT2.GRAD(1,1,JGAU), 343 : & YK,NBNN,NBNNJ,IDIM,RA) 344 : 23 CONTINUE 345 : * 346 : 22 CONTINUE 347 : MIJ=MIJ+1 348 : 349 : * 350 : * 351 : * realisation de l'assemblage 352 : * 353 : DO 240 IX=1,IPT2.NUM(/1) 354 : IA= IPT2.NUM(IX,IELJ) 355 : IB=ICPR2(IA) 356 : DO 240 JX=1,NUM(/1) 357 : IC=NUM(JX,IEL) 358 : ID=ICPR(IC) 359 : RE(IB,ID,1)=RA(IX,JX) + RE(IB,ID,1) 360 : 240 CONTINUE 361 : 362 : 363 : 364 : * 365 : 111 CONTINUE 366 : 110 CONTINUE 367 : * 368 : 10 CONTINUE 369 : * END DO 370 : * on symetrise la matrice 371 : DO 40 IX=1,RE(/1) 372 : DO 40 JX=1,IX 373 : XP =( RE(IX,JX,1) + RE(JX,IX,1)) / 2.D0 374 : RE(IX,JX,1)=XP 375 : RE(JX,IX,1)=XP 376 : 40 CONTINUE 377 : * 378 : * 379 : * DESACTIVATION DES SEGMENTS 380 : * 381 : 999 CONTINUE 382 : MMAT1=MMATX 383 : SEGSUP,MMAT1,MMAT2 384 : 99 CONTINUE 385 : SEGDES,MELEME*NOMOD 386 : SEGDES,MINTE*NOMOD 387 : * SEGDES,INFO*NOMOD 388 : * 389 : MPTVAL=IVAMAT 390 : DO 50 I=1,NMATR 391 : MELVAL=IVAL(I) 392 : SEGDES,MELVAL 393 : 50 CONTINUE 394 : SEGSUP MPTVAL 395 : NOMID=MOMATR 396 : SEGSUP NOMID,MAXT 397 : * 398 : RETURN 399 : END 400 : 401 : 402 : 403 : 404 : 405 : 406 : 407 : 408 : 409 : 410 : 411 : 412 : 413 : 414 : 415 : 416 : 417 : 418 : 419 :
© Cast3M 2003 - All rights reserved.
Disclaimer