1 : C DCOV SOURCE PV 13/04/12 21:15:27 7756 2 : SUBROUTINE DCOV 3 : C 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8(A-H,O-Z) 6 : C 7 : -INC CCOPTIO 8 : -INC SMELEME 9 : -INC SMCOORD 10 : -INC SMRIGID 11 : C 12 : DIMENSION ZLAM(3),ZINVL(3),XCOV(3,3) 13 : CHARACTER*8 MOTS(1),MOTLA(4) 14 : CHARACTER*4 MOTLO(2) 15 : CHARACTER*9 MOTD(1) 16 : C 17 : DATA MOTS/'SIGMA '/ 18 : DATA MOTLA/'LAMBDA ','LAMBDA1 ','LAMBDA2 ','LAMBDA3 '/ 19 : DATA MOTLO/'EXPO','GAUS'/ 20 : DATA MOTD/'DIRECTION'/ 21 : C 22 : C epsilon servant à différents tests 23 : C 24 : EPS = 1.D-12 25 : C 26 : IDIR=0 27 : C 28 : C Lecture du maillage 29 : C 29 : CALL LIROBJ('MAILLAGE',MELEME,1,IRET) 30 : IF (IERR.NE.0) RETURN 31 : C 32 : C Lecture de la loi utilisée 33 : C 33 : CALL LIRMOT(MOTLO,2,ILOI,1) 34 : IF (ILOI.EQ.0) RETURN 35 : C 36 : C Lecture du mot-clé 'SIGMA' 37 : C 37 : CALL LIRMOT(MOTS,1,IRET,1) 38 : IF (IRET.EQ.0) RETURN 39 : C 40 : C Lecture de la valeur de sigma (strictement supérieure à 0.D0) 41 : C 41 : CALL LIRREE(ZSIG,1,IRET) 42 : IF (IERR.NE.0) RETURN 43 : IF (ZSIG.LE.0.D0) THEN 44 : REAERR(1) = REAL(ZSIG) 45 : REAERR(2) = REAL(0.D0) 46 : MOTERR(1:8) = ' SIGMA ' 46 : CALL ERREUR(41) 47 : RETURN 48 : ENDIF 49 : C 50 : C Lecture du mot-clé 'LAMBDA' ou 'LAMBDA1' 51 : C 52 : IF (IDIM.EQ.1) THEN 52 : CALL LIRMOT(MOTLA,1,IANI,1) 53 : ELSE 53 : CALL LIRMOT(MOTLA,2,IANI,1) 54 : ENDIF 55 : IF (IANI.EQ.0) RETURN 56 : C 57 : C Lecture de la valeur de lambda (ou lambda1) 58 : C (strictement supérieure à 0.D0) 59 : C 59 : CALL LIRREE(ZLAM(1),1,IRET) 60 : IF (IERR.NE.0) RETURN 61 : IF (ZLAM(1).LE.0.D0) THEN 62 : REAERR(1) = REAL(ZLAM(1)) 63 : REAERR(2) = REAL(0.D0) 64 : IF (IANI.EQ.1) MOTERR(1:8) = ' LAMBDA ' 65 : IF (IANI.EQ.2) MOTERR(1:8) = ' LAMBDA1' 65 : CALL ERREUR(41) 66 : RETURN 67 : ENDIF 68 : C 69 : IF (IANI.EQ.2) THEN 70 : C 71 : C Lecture éventuelle du mot-clé 'LAMBDA2' (cas anisotrope) 72 : C 72 : CALL LIRMOT(MOTLA(3),1,IRET,1) 73 : IF (IRET.EQ.0) RETURN 74 : C 75 : C Lecture éventuelle de la valeur de lambda2 (cas anisotrope) 76 : C (strictement supérieure à 0.D0) 77 : C 77 : CALL LIRREE(ZLAM(2),1,IRET) 78 : IF (IERR.NE.0) RETURN 79 : IF (ZLAM(2).LE.0.D0) THEN 80 : REAERR(1) = REAL(ZLAM(2)) 81 : REAERR(2) = REAL(0.D0) 82 : MOTERR(1:8) = ' LAMBDA2' 82 : CALL ERREUR(41) 83 : RETURN 84 : ENDIF 85 : C 86 : IF (IDIM.EQ.3) THEN 87 : C 88 : C Lecture éventuelle du mot-clé 'LAMBDA3' (cas anisotrope en 3d) 89 : C 89 : CALL LIRMOT(MOTLA(4),1,IRET,1) 90 : IF (IRET.EQ.0) RETURN 91 : C 92 : C Lecture éventuelle de la valeur de lambda3 (cas anisotrope en 3d) 93 : C (strictement supérieure à 0.D0) 94 : C 94 : CALL LIRREE(ZLAM(3),1,IRET) 95 : IF (IERR.NE.0) RETURN 96 : IF (ZLAM(3).LE.0.D0) THEN 97 : REAERR(1) = REAL(ZLAM(3)) 98 : REAERR(2) = REAL(0.D0) 99 : MOTERR(1:8) = ' LAMBDA3' 99 : CALL ERREUR(41) 100 : RETURN 101 : ENDIF 102 : C 103 : ENDIF 104 : C 105 : C Lecture éventuelle et optionnelle du mot-clé 'DIRECTION' 106 : C (cas anisotrope, directions d'anisotropie non parallèles aux axes) 107 : C 107 : CALL LIRMOT(MOTD,1,IDIR,0) 108 : IF (IERR.NE.0) RETURN 109 : C 110 : IF (IDIR.NE.0) THEN 111 : C 112 : C Lecture éventuelle de la valeur de VEC1 113 : C (cas anisotrope, directions d'anisotropie non parallèles aux axes) 114 : C 114 : CALL LIROBJ('POINT',IPTV1,1,IRET) 115 : IF (IERR.NE.0) RETURN 116 : C 117 : IF (IDIM.EQ.3) THEN 118 : C 119 : C Lecture éventuelle de la valeur de VEC2 120 : C (cas anisotrope en 3d, 121 : C directions d'anisotropie non parallèles aux axes) 122 : C 122 : CALL LIROBJ('POINT',IPTV2,1,IRET) 123 : IF (IERR.NE.0) RETURN 124 : C 125 : ENDIF 126 : C 127 : ENDIF 128 : C 129 : ENDIF 130 : C 131 : C Dans le cas anisotrope, 132 : C directions d'anisotropie non parallèles aux axes, 133 : C on vérifie que le(s) vecteur(s) introduit(s) sous forme de point ne 134 : C sont pas de longueur nulle 135 : C si tel est le cas, on le(s) norme 136 : C 137 : IF (IDIR.NE.0) THEN 138 : C 139 : IREF1 = (IPTV1-1)*(IDIM+1) 140 : SCOV1 = 0.D0 141 : DO 1 I=1,IDIM 142 : XCOV(1,I) = XCOOR(IREF1+I) 143 : SCOV1 = SCOV1+(XCOV(1,I)*XCOV(1,I)) 144 : 1 CONTINUE 145 : C 146 : IF (SCOV1.LT.EPS) THEN 146 : CALL ERREUR(239) 147 : RETURN 148 : ELSE 149 : DO 2 I=1,IDIM 150 : XCOV(1,I) = XCOV(1,I)/SCOV1 151 : 2 CONTINUE 152 : ENDIF 153 : C 154 : IF (IDIM.EQ.3) THEN 155 : C 156 : IREF2 = (IPTV2-1)*(IDIM+1) 157 : SCOV2 = 0.D0 158 : DO 3 I=1,IDIM 159 : XCOV(2,I) = XCOOR(IREF2+I) 160 : SCOV2 = SCOV2+(XCOV(2,I)*XCOV(2,I)) 161 : 3 CONTINUE 162 : C 163 : IF (SCOV2.LT.EPS) THEN 163 : CALL ERREUR(239) 164 : RETURN 165 : ELSE 166 : DO 4 I=1,IDIM 167 : XCOV(2,I) = XCOV(2,I)/SCOV2 168 : 4 CONTINUE 169 : ENDIF 170 : C 171 : C Dans le cas anisotrope 3d, 172 : C directions d'anisotropie non parallèles aux axes, 173 : C on vérifie que les vecteurs introduits sous forme de point 174 : C sont bien orthogonaux 175 : C 176 : SCOV12 = 0.D0 177 : DO 5 I=1,IDIM 178 : SCOV12 = SCOV12+(XCOV(1,I)*XCOV(2,I)) 179 : 5 CONTINUE 180 : IF (SCOV12.GT.EPS) THEN 180 : CALL ERREUR(161) 181 : RETURN 182 : ENDIF 183 : C 184 : ENDIF 185 : C 186 : ENDIF 187 : C 188 : IF (IDIR.NE.0) THEN 189 : C 190 : C Dans le cas anisotrope 2d, 191 : C directions d'anisotropie non parallèles aux axes, 192 : C on complète le repère par rotation de +90 degrés 193 : C 194 : IF (IDIM.EQ.2) THEN 195 : XCOV(2,1) = -1.D0*XCOV(1,2) 196 : XCOV(2,2) = XCOV(1,1) 197 : ENDIF 198 : C 199 : C Dans le cas anisotrope 3d, 200 : C directions d'anisotropie non parallèles aux axes, 201 : C on complète le repère par produit vectoriel des deux premiers 202 : C vecteurs 203 : C 204 : IF (IDIM.EQ.3) THEN 205 : XCOV(3,1) = (XCOV(1,2)*XCOV(2,3))-(XCOV(1,3)*XCOV(2,2)) 206 : XCOV(3,2) = (XCOV(1,3)*XCOV(2,1))-(XCOV(1,1)*XCOV(2,3)) 207 : XCOV(3,3) = (XCOV(1,1)*XCOV(2,2))-(XCOV(1,2)*XCOV(2,1)) 208 : ENDIF 209 : C 210 : ENDIF 211 : C 212 : C Il vaut mieux faire * (1./lambda) que / lambda 213 : C 214 : ZINVL(1) = 1.D0/ZLAM(1) 215 : IF (IANI.EQ.2) THEN 216 : ZINVL(2) = 1.D0/ZLAM(2) 217 : IF (IDIM.EQ.3) ZINVL(3) = 1.D0/ZLAM(3) 218 : ENDIF 219 : C 220 : C Les éléments du MELEME sont transformés en POI1 si nécessaire 221 : C 222 : SEGACT MELEME 223 : IPT1=MELEME 224 : IPT2=MELEME 225 : IF(ITYPEL.NE.1.OR.(LISOUS(/1)).NE.0) THEN 225 : CALL CHANGE(IPT1,1) 226 : IF (IERR.NE.0) RETURN 227 : ENDIF 228 : C 229 : NBLI=IPT1.NUM(/2) 230 : C 231 : C On définit un super-élément correspondant aux noeuds 232 : C 233 : NBSOUS = 0 234 : NBREF = 0 235 : NBNN = NBLI 236 : NBELEM = 1 237 : SEGINI MELEME 238 : ITYPEL = 28 239 : DO 7 I=1,NBLI 240 : NUM(I,1) = IPT1.NUM(1,I) 241 : 7 CONTINUE 242 : SEGDES IPT2 243 : C 244 : C On met la matrice dans un objet rigidité 245 : C 246 : NRIGE = 8 247 : NRIGEL = 1 248 : SEGINI MRIGID 249 : ICHOLE = 0 250 : IMGEO1 = 0 251 : IMGEO2 = 0 252 : IFORIG = IFOUR 253 : ISUPEQ = 0 254 : NELRIG = 1 255 : COERIG(1) = 1.D0 256 : * SEGINI IMATRI 257 : NLIGRD = NBLI 258 : NLIGRP = NBLI 259 : segini xmatri 260 : SEGINI DESCR 261 : IRIGEL(1,1) = MELEME 262 : IRIGEL(2,1) = 0 263 : IRIGEL(3,1) = DESCR 264 : IRIGEL(4,1) = xMATRI 265 : IRIGEL(5,1) = 0 266 : IRIGEL(6,1) = 0 267 : DO 8 I=1,NBLI 268 : NOELEP(I) = I 269 : NOELED(I) = I 270 : LISINC(I) = 'SCAL' 271 : LISDUA(I) = 'SCAL' 272 : 8 CONTINUE 273 : SEGDES DESCR,MRIGID 274 : * SEGINI XMATRI 275 : * IMATTT(1) = XMATRI 276 : * DO 9 I=1,NBLI 277 : * DO 10 J=1,NBLI 278 : * RE(I,J,1) = 0.D0 279 : * 10 CONTINUE 280 : * 9 CONTINUE 281 : C 282 : IF (ILOI.EQ.1) THEN 283 : C 284 : C Loi exponentielle 285 : C 286 : IF (IANI.EQ.1) THEN 287 : C 288 : C Cas isotrope 289 : C 290 : DO 100 I=1,NBLI 291 : IREF=(NUM(I,1)-1)*(IDIM+1) 292 : DO 110 J=1,I 293 : JREF=(NUM(J,1)-1)*(IDIM+1) 294 : VM=0.D0 295 : DO 120 I1=1,IDIM 296 : VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2 297 : 120 CONTINUE 298 : VM=SQRT(VM) 299 : RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1)) 300 : 110 CONTINUE 301 : 100 CONTINUE 302 : C 303 : ELSEIF (IANI.EQ.2) THEN 304 : C 305 : C Cas anisotrope : 306 : C directions d'anisotropie parallèles aux axes (IDIR nul) 307 : C 308 : IF (IDIR.EQ.0) THEN 309 : C 310 : DO 200 I=1,NBLI 311 : IREF=(NUM(I,1)-1)*(IDIM+1) 312 : DO 210 J=1,I 313 : JREF=(NUM(J,1)-1)*(IDIM+1) 314 : VM=0.D0 315 : DO 220 I1=1,IDIM 316 : VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2 317 : 220 CONTINUE 318 : VM=SQRT(VM) 319 : RE(I,J,1)=ZSIG*ZSIG*EXP(-VM) 320 : 210 CONTINUE 321 : 200 CONTINUE 322 : C 323 : C Cas anisotrope : 324 : C directions d'anisotropie non parallèles aux axes (IDIR non nul) 325 : C 326 : ELSEIF (IDIR.NE.0) THEN 327 : C 328 : DO 300 I=1,NBLI 329 : IREF=(NUM(I,1)-1)*(IDIM+1) 330 : DO 310 J=1,I 331 : JREF=(NUM(J,1)-1)*(IDIM+1) 332 : VM=0.D0 333 : DO 320 I1=1,IDIM 334 : PSCAL=0.D0 335 : DO 330 J1=1,IDIM 336 : PSCAL = PSCAL + 337 : * (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1) 338 : 330 CONTINUE 339 : VM=VM+(PSCAL*ZINVL(I1))**2 340 : 320 CONTINUE 341 : VM=SQRT(VM) 342 : RE(I,J,1)=ZSIG*ZSIG*EXP(-VM) 343 : 310 CONTINUE 344 : 300 CONTINUE 345 : C 346 : ENDIF 347 : C 348 : ENDIF 349 : C 350 : ELSEIF (ILOI.EQ.2) THEN 351 : C 352 : C Loi gaussienne 353 : C 354 : IF (IANI.EQ.1) THEN 355 : C 356 : C Cas isotrope 357 : C 358 : DO 400 I=1,NBLI 359 : IREF=(NUM(I,1)-1)*(IDIM+1) 360 : DO 410 J=1,I 361 : JREF=(NUM(J,1)-1)*(IDIM+1) 362 : VM=0.D0 363 : DO 420 I1=1,IDIM 364 : VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2 365 : 420 CONTINUE 366 : RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1)) 367 : 410 CONTINUE 368 : 400 CONTINUE 369 : C 370 : ELSEIF (IANI.EQ.2) THEN 371 : C 372 : C Cas anisotrope : 373 : C directions d'anisotropie parallèles aux axes (IDIR nul) 374 : C 375 : IF (IDIR.EQ.0) THEN 376 : C 377 : DO 500 I=1,NBLI 378 : IREF=(NUM(I,1)-1)*(IDIM+1) 379 : DO 510 J=1,I 380 : JREF=(NUM(J,1)-1)*(IDIM+1) 381 : VM=0.D0 382 : DO 520 I1=1,IDIM 383 : VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2 384 : 520 CONTINUE 385 : RE(I,J,1)=ZSIG*ZSIG*EXP(-VM) 386 : 510 CONTINUE 387 : 500 CONTINUE 388 : C 389 : C Cas anisotrope : 390 : C directions d'anisotropie non parallèles aux axes (IDIR non nul) 391 : C 392 : ELSEIF (IDIR.NE.0) THEN 393 : C 394 : DO 600 I=1,NBLI 395 : IREF=(NUM(I,1)-1)*(IDIM+1) 396 : DO 610 J=1,I 397 : JREF=(NUM(J,1)-1)*(IDIM+1) 398 : VM=0.D0 399 : DO 620 I1=1,IDIM 400 : PSCAL=0.D0 401 : DO 630 J1=1,IDIM 402 : PSCAL = PSCAL + 403 : * (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1) 404 : 630 CONTINUE 405 : VM=VM+(PSCAL*ZINVL(I1))**2 406 : 620 CONTINUE 407 : RE(I,J,1)=ZSIG*ZSIG*EXP(-VM) 408 : 610 CONTINUE 409 : 600 CONTINUE 410 : C 411 : ENDIF 412 : C 413 : ENDIF 414 : C 415 : ENDIF 416 : C 417 : C Calcul de M telle que (RE =) COV = M Mt 418 : C M est stockée dans la partie inférieure de RE 419 : C 419 : CALL SYCHOM(RE,NBLI) 420 : C 421 : C Création de l'objet rigidité 422 : C 422 : CALL ECROBJ('RIGIDITE',MRIGID) 423 : C 424 : SEGDES XMATRI,MELEME 425 : C 426 : END 427 : 428 : 429 : 430 : 431 : 432 : 433 : 434 : 435 : 436 : 437 :
© Cast3M 2003 - All rights reserved.
Disclaimer