1 : C BLOQUE SOURCE PV 13/04/11 21:15:26 7755 2 : 3 : C----------------------------------------------------------------------- 4 : C Cet operateur impose les BLOCAGES 5 : C 6 : C Syntaxe : 7 : C 8 : C ENC1 = BLOQUER ( DEPL ) ( ROTA ) POI1 9 : C 10 : C ou ENC1 = BLOQUER RADIAL P1 (P2) MELEME 11 : C ORTHOR P1 (P2) MELEME 12 : C 13 : C ou ENC1 = BLOQUER (DEPL) (ROTA) DIRECTION V1 MELEME 14 : C 15 : C DIM = 1 ( UX UY UZ ) ou ( UR UZ ) | 16 : C DIM = 2 OU 3 ( UX UY UZ RX RY RZ ) | MELEME 17 : C AXISYM ( RX RZ RT UT ) | 18 : Clist t 19 : C POI1 = OBJET DE TYPE POINT (entree) 20 : C MELEME = OBJET DE TYPE MELEME (entree) 21 : C ENC1 = OBJET DE TYPE RIGIDITE (sortie) 22 : C 23 : C 1) On peut imposer des BLOQUAGES UNILATERAUX en specifiant les 24 : C mots-cles MINIMUM ou MAXIMUM. 25 : C 2) On peut imposer des BLOQUAGES de type FROTTEMENT en specifiant 26 : C le mot-cle FROTTEMENT. 27 : C----------------------------------------------------------------------- 28 : C Juillet 2003 : passage a un seul multiplicateur 29 : C----------------------------------------------------------------------- 30 : 31 : SUBROUTINE BLOQUE 32 : 33 : IMPLICIT INTEGER(I-N) 34 : IMPLICIT REAL*8 (A-H,O-Z) 35 : 36 : -INC CCOPTIO 37 : -INC CCGEOME 38 : -INC CCREEL 39 : -INC CCHAMP 40 : 41 : -INC SMCHPOI 42 : -INC SMCOORD 43 : -INC SMELEME 44 : -INC SMLMOTS 45 : POINTEUR LMDEPL.MLMOTS 46 : -INC SMMODEL 47 : -INC SMRIGID 48 : -INC SMTABLE 49 : -INC TMTRAV 50 : 51 : SEGMENT MSWBLO 52 : CHARACTER*4 MOTDDL(0) 53 : ENDSEGMENT 54 : 55 : segment lispoi 56 : INTEGER pilpoi(mpoin),pilmul(mpoin) 57 : endsegment 58 : 59 : DIMENSION XNOR(3),U1(3),U2(3) 60 : 61 : CHARACTER*4 CHADDL 62 : CHARACTER*8 MOTRIG 63 : CHARACTER*4 MOTPV(3), MOTBLO(5) 64 : CHARACTER*4 MODEPL(6),MODEDU(6), MOROTA(5),MORODU(5) 65 : CHARACTER*4 MODE1D(2),MOFO1D(2) 66 : 67 : DATA EPSI / 1.D-12 / 68 : DATA MOTRIG / 'RIGIDITE' / 69 : DATA MOTPV / 'MINI','MAXI','FROT' / 70 : DATA MOTBLO / 'DEPL','ROTA','RADI','ORTH','DIRE' / 71 : DATA MODEPL / 'UX ','UY ','UZ ','UR ','UZ ','UT ' / 72 : DATA MODEDU / 'FX ','FY ','FZ ','FR ','FZ ','FT ' / 73 : DATA MOROTA / 'RX ','RY ','RZ ','RT ','RS ' / 74 : DATA MORODU / 'MX ','MY ','MZ ','MT ','MS ' / 75 : C Tableaux MODE1D et MOFO1D sont utilises pour certains modes 1D 76 : DATA MODE1D / 'UX ','UZ ' / 77 : DATA MOFO1D / 'FX ','FZ ' / 78 : 79 : C Pour ne pas avoir de verrouillage sur MCOORD en // 80 : SEGDES,MCOORD 81 : SEGACT,MCOORD*MOD 82 : 83 : c lecture table liaisons statiques 83 : CALL LIRTAB('LIAISONS_STATIQUES',ipt,0,iretou) 84 : c traitement table liaisons statiques 85 : IF (iretou.NE.0) THEN 85 : CALL BLOQU2(ipt) 86 : RETURN 87 : ENDIF 88 : 89 : C Est-ce une condition unilaterale ? 90 : NILATE=0 90 : CALL LIRMOT(MOTPV,3,IPO,0) 91 : IF (IPO.EQ.1) NILATE=-1 92 : IF (IPO.EQ.2) NILATE=1 93 : IF (IPO.EQ.3) NILATE=2 94 : C Pas de frottement en 1D 95 : IF (IPO.EQ.3.AND.IDIM.EQ.1) THEN 96 : INTERR(1)=IDIM 97 : MOTERR(1:4)=MOTPV(3) 97 : CALL ERREUR(971) 98 : RETURN 99 : ENDIF 100 : 101 : idimp1=IDIM+1 102 : C Quelques initialisations selon le type de probleme 103 : C Cas IDIM = 1 : 104 : C ISPE1D = 1 si IDIM=1 et IFOUR=9 ou 10, car les noms de DDL primaux et 105 : C variables duales ne sont pas dans l'ordre "classique" (un traitement 106 : C specifique est alors necessaire). 107 : ISPE1D=0 108 : C Deformations planes ou contraintes planes ou defo. plane gene : 109 : IF (IFOUR.EQ.-1.OR.IFOUR.EQ.-2.OR.IFOUR.EQ.-3) THEN 110 : LDEPL=2 111 : IADEPL=0 112 : LROTA=1 113 : IAROTA=2 114 : C Axisymetrique : 115 : ELSE IF (IFOUR.EQ.0) THEN 116 : LDEPL=2 117 : IADEPL=3 118 : LROTA=1 119 : IAROTA=3 120 : C Fourier : 121 : ELSE IF (IFOUR.EQ.1) THEN 122 : LDEPL=3 123 : IADEPL=3 124 : LROTA=1 125 : IAROTA=3 126 : C Tridimensionnel : 127 : ELSE IF (IFOUR.EQ.2) THEN 128 : LDEPL=3 129 : IADEPL=0 130 : LROTA=3 131 : IAROTA=0 132 : C Massif 1D (IDIM=1) : 133 : ELSE IF (IFOUR.GE.3.AND.IFOUR.LE.15) THEN 134 : IF (IFOUR.LE.6) THEN 135 : LDEPL=1 136 : IADEPL=0 137 : ELSE IF (IFOUR.GE.7.AND.IFOUR.LE.10) THEN 138 : LDEPL=2 139 : IADEPL=0 140 : IF (IFOUR.EQ.9.OR.IFOUR.EQ.10) ISPE1D=1 141 : ELSE IF (IFOUR.EQ.11) THEN 142 : LDEPL=3 143 : IADEPL=0 144 : ELSE IF (IFOUR.EQ.15) THEN 145 : LDEPL=2 146 : IADEPL=3 147 : ELSE 148 : LDEPL=1 149 : IADEPL=3 150 : ENDIF 151 : LROTA=0 152 : IAROTA=0 153 : C Autres cas : 154 : ELSE 155 : LDEPL=0 156 : IADEPL=0 157 : LROTA=0 158 : IAROTA=0 159 : ENDIF 160 : 161 : C Initialisation de la liste des ddls MOTDDL (segment MSWBLO) 162 : SEGINI,MSWBLO 163 : 164 : C----------------------------------------------------------------------- 165 : C Lecture eventuelle des MOTS autres que des DDLS 166 : C----------------------------------------------------------------------- 167 : C Lecture eventuelle de 'DEPL' et/ou 'ROTA' 168 : C -------------------- 169 : IDEPL=0 170 : IROTA=0 170 : 481 CALL LIRMOT(MOTBLO,2,IMOT,0) 171 : IF (IMOT.EQ.1) IDEPL=1 172 : IF (IMOT.EQ.2) IROTA=1 173 : IF (IMOT.NE.0) GOTO 481 174 : C Lecture eventuelle de 'RADI','ORTH' ou 'DIRE' 175 : C -------------------- 176 : IRADIA=0 177 : IDIREC=0 178 : IDICHP=0 178 : 4480 CALL LIRMOT(MOTBLO(3),3,IMOT,0) 179 : IF (IMOT.EQ.0) THEN 180 : IF (IDEPL.EQ.1) GOTO 44801 181 : IF (IROTA.EQ.1) GOTO 44802 182 : IBDDL=MOTDDL(/2) 183 : IF (IBDDL.NE.0) GOTO 449 184 : IF (IBDDL.EQ.0) GOTO 445 185 : ENDIF 186 : C En DIMENSION 1, les mots-cles 'RADI,'ORTH' et 'DIRE' sont interdits. 187 : IF (IDIM.EQ.1) THEN 188 : INTERR(1)=IDIM 189 : MOTERR(1:4)=MOTBLO(2+IMOT) 189 : CALL ERREUR(971) 190 : GOTO 1000 191 : ENDIF 192 : GO TO (44803,44803,44804),IMOT 193 : C Traitement des mots-cles : Mise a jour de MOTDDL 194 : C ---------------------------- 195 : C On a trouve le mot DEPLAcement 196 : 44801 IDEPL=0 197 : C Cas particulier pour certains modes de IDIM=1 198 : IF (ISPE1D.EQ.1) THEN 199 : DO i=1,LDEPL 200 : MOTDDL(**)=MODE1D(IADEPL+i) 201 : MOTDDL(**)=MOFO1D(IADEPL+i) 202 : ENDDO 203 : C Cas general 204 : ELSE 205 : DO i=1,LDEPL 206 : MOTDDL(**)=MODEPL(IADEPL+i) 207 : MOTDDL(**)=MODEDU(IADEPL+i) 208 : ENDDO 209 : ENDIF 210 : GOTO 4480 211 : C On a trouve le mot ROTAtion 212 : 44802 IROTA=0 213 : DO i=1,LROTA 214 : MOTDDL(**)=MOROTA(IAROTA+i) 215 : MOTDDL(**)=MORODU(IAROTA+i) 216 : ENDDO 217 : GOTO 4480 218 : C On a trouve le mot RADial ou le mot ORTHoradial 219 : 44803 IRADIA=IMOT 219 : CALL LIROBJ('POINT',KPOINT,1,IRETOU) 220 : IF (IRETOU.EQ.0) GOTO 1000 221 : j=(KPOINT-1)*idimp1 222 : DO i=1,IDIM 223 : U1(i)=XCOOR(j+i) 224 : ENDDO 225 : IF (IDIM.EQ.3) THEN 226 : C Lecture du 2nd point de l'axe (en 3D) 226 : CALL LIROBJ('POINT',KPOINT,1,IRETOU) 227 : IF (IRETOU.EQ.0) GOTO 1000 228 : j=(KPOINT-1)*idimp1 229 : YL=0.D0 230 : DO i=1,IDIM 231 : U2(i)=XCOOR(j+i)-U1(i) 232 : YL=YL+U2(i)*U2(i) 233 : ENDDO 234 : C Calcul du vecteur directeur unitaire de l'axe (U2) 235 : IF (YL.LT.EPSI) THEN 235 : CALL ERREUR(237) 236 : GOTO 1000 237 : ENDIF 238 : YL=1.D0/SQRT(YL) 239 : DO i=1,IDIM 240 : U2(i)=U2(i)*YL 241 : ENDDO 242 : ENDIF 243 : GOTO 449 244 : C On a trouve le mot DIREction 245 : 44804 IDIREC=1 245 : CALL LIROBJ('POINT',KPOINT,0,IRETOU) 246 : IF (IRETOU.EQ.0) THEN 247 : * WRITE(IOIMP,*) 'Lecture dun chpoint' 247 : CALL LIROBJ('CHPOINT',MCHPOI,1,IRETOU) 248 : IF (IRETOU.EQ.0) GOTO 1000 249 : IDICHP=1 250 : ELSE 251 : j=(KPOINT-1)*idimp1 252 : YL=0.D0 253 : DO i=1,IDIM 254 : XNOR(i)=XCOOR(j+i) 255 : YL=YL+XNOR(i)*XNOR(i) 256 : ENDDO 257 : IF (YL.LT.EPSI) THEN 257 : CALL ERREUR(239) 258 : GOTO 1000 259 : ENDIF 260 : YL=1.D0/SQRT(YL) 261 : DO i=1,IDIM 262 : XNOR(i)=XNOR(i)*YL 263 : ENDDO 264 : ENDIF 265 : GOTO 449 266 : 267 : C Lecture eventuelle de DDLs : 268 : C ------------------------------ 269 : 445 CONTINUE 270 : C La liste des grandeurs PRIMALEs/DDLs est un LISTMOTS. 271 : C La liste des grandeurs DUALEs est un LISTMOTS obligatoire. 271 : CALL LIROBJ('LISTMOTS',MLMOT1,0,IRETOU) 272 : IF (IERR.NE.0) GOTO 1000 273 : IF (IRETOU.EQ.0) GOTO 446 273 : CALL LIROBJ('LISTMOTS',MLMOT2,1,IRETOU) 274 : IF (IERR.NE.0) GOTO 1000 275 : SEGACT,MLMOT1,MLMOT2 276 : IMOT1 = MLMOT1.MOTS(/2) 277 : IF (IMOT1.LE.0) THEN 277 : CALL ERREUR(643) 278 : SEGDES,MLMOT1,MLMOT2 279 : GOTO 1000 280 : ELSE IF (IMOT1.NE.MLMOT2.MOTS(/2)) THEN 280 : CALL ERREUR(854) 281 : SEGDES,MLMOT1,MLMOT2 282 : GOTO 1000 283 : ENDIF 284 : DO IMOT = 1, IMOT1 285 : MOTDDL(**)=MLMOT1.MOTS(IMOT) 286 : MOTDDL(**)=MLMOT2.MOTS(IMOT) 287 : ENDDO 288 : SEGDES,MLMOT1,MLMOT2 289 : GOTO 449 290 : C La liste des ddls autorises NOMDD est dans BDATA.ESO 291 : C On doit lire au moins 1 ddl (car sinon MOTDDL est vide !) 292 : 446 CONTINUE 292 : CALL LIRMOT(NOMDD,LNOMDD,IMOT,0) 293 : IF (IERR.NE.0) GOTO 1000 294 : IF (IMOT.NE.0) THEN 295 : MOTDDL(**)=NOMDD(IMOT) 296 : MOTDDL(**)=NOMDU(IMOT) 297 : ELSE 297 : CALL LIRCHA(CHADDL,0,IMOT) 298 : IF (IERR.NE.0) GOTO 1000 299 : IF (IMOT.EQ.0) GOTO 447 300 : MOTDDL(**)=CHADDL 301 : MOTDDL(**)=' ' 302 : ENDIF 303 : GOTO 446 304 : 447 CONTINUE 305 : IBDDL=MOTDDL(/2) 306 : IPMODL = 0 306 : CALL LIROBJ('MMODEL ',IPMODL,0,IRETOU) 307 : IF (IERR.NE.0) GOTO 1000 308 : IF (IPMODL.EQ.0) THEN 309 : DO i = 2, IBDDL, 2 310 : IF (MOTDDL(i).EQ.' ') THEN 311 : MOTERR(1:4)=MOTDDL(i-1) 311 : CALL ERREUR(108) 312 : ENDIF 313 : ENDDO 314 : ELSE 314 : CALL NOVARD(IPMODL,'DEPL') 314 : CALL LIROBJ('LISTMOTS',MLMOT1,1,IRETOU) 315 : IF (IERR.NE.0) GOTO 1000 315 : CALL NOVARD(IPMODL,'FORC') 315 : CALL LIROBJ('LISTMOTS',MLMOT2,1,IRETOU) 316 : IF (IERR.NE.0) GOTO 1000 317 : SEGACT,MLMOT1,MLMOT2 318 : IMOT1 = MLMOT1.MOTS(/2) 319 : DO i = 2, IBDDL, 2 320 : IF (MOTDDL(i).EQ.' ') THEN 320 : CALL PLACE(MLMOT1.MOTS(1),IMOT1,iplac,MOTDDL(i-1)) 321 : IF (iplac.EQ.0) THEN 322 : MOTERR(1:4)=MOTDDL(i-1) 322 : CALL ERREUR(197) 323 : ELSE 324 : MOTDDL(i)=MLMOT2.MOTS(iplac) 325 : ENDIF 326 : ENDIF 327 : ENDDO 328 : SEGSUP,MLMOT1,MLMOT2 329 : ENDIF 330 : IF (IERR.NE.0) GOTO 1000 331 : 332 : 449 IBDDL=MOTDDL(/2) 333 : C Verification que le nombre de DDLs a bloquer n'est pas nul 334 : C IF (IBDDL.EQ.0) GOTO 1000 335 : C----------------------------------------------------------------------- 336 : C Fin de la lecture des mots (DEPL,ROTA...) ou des DDLs 337 : C----------------------------------------------------------------------- 338 : 339 : C Recherche du maillage MELEME de type POINT : 340 : C ---------------------------------------------- 341 : KOBJET=0 342 : MELEME=0 343 : C On cherche d'abord si on a un POINT que l'on transformera en POI1 344 : C sinon on cherche un maillage que l'on transforme en POI1 si besoin 344 : CALL LIROBJ('POINT',KPOINT,0,IRETOU) 345 : IF (IERR.NE.0) GOTO 1000 346 : IF (IRETOU.EQ.0) THEN 346 : CALL LIROBJ('MAILLAGE',KOBJET,1,IRETOU) 347 : IF (IERR.NE.0) GOTO 1000 348 : MELEME=KOBJET 349 : SEGACT,MELEME 349 : IF (ITYPEL.NE.1) CALL CHANGE(MELEME,1) 350 : NBPOIN=NUM(/2) 351 : ELSE 352 : MELEME=KPOINT 352 : CALL CRELEM(MELEME) 353 : NBPOIN=1 354 : C** Ce bloc est ancien mais conserve. 355 : C** Cas du POINT SUPPORT des DDLS de DEFO PLANE GENE 356 : C** On rajoute les DDLS globaux si 'DEPL' ou 'ROTA' 357 : C** IF (KPOINT.EQ.NSDPGE) THEN 358 : C** IF (IDGDEP.EQ.1) THEN 359 : C** MOTDDL(**)=MODEPL(5) 360 : C** MOTDDL(**)=MODEDU(5) 361 : C** ENDIF 362 : C** IF (IDGROT.EQ.1) THEN 363 : C** MOTDDL(**)=MOROTA(1) 364 : C** MOTDDL(**)=MORODU(1) 365 : C** MOTDDL(**)=MOROTA(2) 366 : C** MOTDDL(**)=MORODU(2) 367 : C** ENDIF 368 : C** IBDDL=MOTDDL(/2) 369 : C** ENDIF 370 : ENDIF 371 : * 372 : IF (IDICHP.EQ.1) THEN 373 : * Cas où la direction est donnée par un chpoint 374 : * On construit un segment MTRAV contenant les directions normées 375 : * WRITE(IOIMP,*) 'Transfo du chpoint en MTRA' 376 : * 1) Recopie des composantes et valeurs pertinentes du chpoint 377 : * dans le TMTRAV 378 : * Composantes recherchées (LMDEPL.MLMOTS) 379 : * sur les points du maillage MELEME 380 : JGN=4 381 : JGM=IDIM 382 : SEGINI,LMDEPL 383 : DO I=1,IDIM 384 : LMDEPL.MOTS(I)=MODEPL(IADEPL+I) 385 : ENDDO 385 : CALL CP2TR2(LMDEPL,MELEME,MCHPOI,MTRAV,IRET) 386 : SEGSUP LMDEPL 387 : SEGACT MELEME 388 : IF (IRET.NE.0) THEN 388 : CALL ERREUR(5) 389 : GOTO 1000 390 : ENDIF 391 : * 2) Calcul des normes 392 : * WRITE(IOIMP,*) 'Calcul des normes ' 393 : SEGACT MTRAV*MOD 394 : DO IBPOIN=1,NBPOIN 395 : YL=0.D0 396 : DO I=1,IDIM 397 : XNOR(I)=BB(I,IBPOIN) 398 : YL=YL+XNOR(I)*XNOR(I) 399 : ENDDO 400 : IF (YL.LT.XPETIT) THEN 400 : CALL ERREUR(239) 401 : GOTO 1000 402 : ENDIF 403 : YL=1.D0/SQRT(YL) 404 : DO I=1,IDIM 405 : BB(I,IBPOIN)=XNOR(I)*YL 406 : ENDDO 407 : ENDDO 408 : * WRITE(IOIMP,*) 'Fin du calcul des normes ' 409 : * SEGPRT,MTRAV 410 : ENDIF 411 : C Determination du nombre de multiplicateurs NNMAT par noeud de MELEME 412 : C NNMAT correspond au nombre de DDLs a bloquer par noeud = nombre de 413 : C multiplicateurs a creer par noeud (1 multiplicateur) = NRIGEL 414 : C Dans les cas RADIal et DIREction, on a une seule matrice par noeud. 415 : C Dans les autres cas, autant de matrices que MOTDDL(/2)/2. 416 : NNMAT=1 417 : IF (IDIREC+IRADIA.EQ.0) NNMAT=IBDDL/2 418 : 419 : C Initialisation de l'objet RIGIDITE associe aux BLOCAGES 420 : NRIGE=8 421 : NRIGEL=NNMAT 422 : SEGINI,MRIGID 423 : MTYMAT=MOTRIG 424 : C* IFORIG=IFOMOD 425 : IFORIG=IFOUR 426 : ICHOLE=0 427 : IMGEO1=0 428 : IMGEO2=0 429 : KRIGI=MRIGID 430 : 431 : C A chaque multplicateur est associe un nouveau noeud 432 : C NBPOIN : nombre de points du maillage MELEME a bloquer 433 : NBNO=XCOOR(/1)/idimp1 434 : NBNOI=NBNO 435 : NBPTS=NBNO+NNMAT*NBPOIN 436 : SEGADJ,MCOORD 437 : 438 : C Boucle sur le nombre de DDLs a bloquer 439 : DO IAA=1,NNMAT 440 : C Pour chaque noeud du maillage, creation d'un noeud associe (place 441 : C initialement a l'origine) au IAA-eme DDL 442 : DO i=1,NBPOIN 443 : DO j=1,idimp1 444 : XCOOR(NBNOI*idimp1+j)=0.D0 445 : ENDDO 446 : NBNOI=NBNOI+1 447 : ENDDO 448 : C Creation du maillage MELEME de MULTiplicateurs associe aux BLOCAGES 449 : NBSOUS=0 450 : NBREF=0 451 : NBNN=2 452 : NBELEM=NBPOIN 453 : SEGINI,IPT1 454 : KIPT1=IPT1 455 : IPT1.ITYPEL=22 456 : DO i=1,NBPOIN 457 : j=(IAA-1)*NBPOIN+i 458 : IPT1.NUM(1,i)=NBNO+j 459 : IPT1.NUM(2,i)=NUM(1,i) 460 : IPT1.ICOLOR(i)=IDCOUL 461 : C Correction : Les pts mult sont a la meme position que les noeuds 462 : IREF3=(NUM(1,i)-1)*idimp1 463 : IREF1=(NBNO+j-1)*idimp1 464 : DO j=1,IDIM 465 : XCOOR(IREF1+j)=XCOOR(IREF3+j) 466 : ENDDO 467 : ENDDO 468 : C Creation des RAIDEURS associees au IAA-eme multplicateur (DDL) 469 : IRIGEL(1,IAA)=KIPT1 470 : IRIGEL(2,IAA)=0 471 : IRIGEL(5,IAA)=NIFOUR 472 : IRIGEL(6,IAA)=NILATE 473 : NLIGRP=2 474 : IF (IDIREC+IRADIA.NE.0) NLIGRP=1+LDEPL 475 : NLIGRD=NLIGRP 476 : C Remplissage du tableau des DESCripteurs de RIG 477 : SEGINI,DESCR 478 : IRIGEL(3,IAA)=DESCR 479 : NOELEP(1)=1 480 : NOELEP(2)=2 481 : NOELED(1)=1 482 : NOELED(2)=2 483 : IF (IDIREC+IRADIA.EQ.0) THEN 484 : j=2*(IAA-1) 485 : LISINC(1)='LX ' 486 : LISINC(2)=MOTDDL(j+1) 487 : LISDUA(1)='FLX ' 488 : LISDUA(2)=MOTDDL(j+2) 489 : ELSE 490 : DO i=1,LDEPL 491 : NOELEP(1+i)=2 492 : NOELED(1+i)=2 493 : IF (IROTA.NE.1) THEN 494 : LISINC(1+i)=MODEPL(IADEPL+i) 495 : LISDUA(1+i)=MODEDU(IADEPL+i) 496 : ELSE 497 : LISINC(1+i)=MOROTA(IADEPL+i) 498 : LISDUA(1+i)=MORODU(IADEPL+i) 499 : ENDIF 500 : ENDDO 501 : LISINC(1)='LX ' 502 : LISDUA(1)='FLX ' 503 : ENDIF 504 : SEGDES DESCR 505 : 506 : C** NLIGRP=3 507 : C** IF (IDIREC+IRADIA.NE.0) NLIGRP=IDIM+2 508 : C** NLIGRD=NLIGRP 509 : NELRIG=NBPOIN 510 : SEGINI,XMATRI 511 : IRIGEL(4,IAA)=XMATRI 512 : COERIG(IAA)=1.D0 513 : C Remplissage de la matrice de rigidite RE : 514 : C Il faut distinguer les cas IRADIA, IDIREC et les autres 515 : C IRADIA : Il faut calculer la DIREction puis identique a IDIREC 516 : C IDIREC : La DIRECTION est stockee dans le vecteur norme XNOR 517 : C AUTRES : La matrice est predefinie dans RIG 518 : C Option RADIAL : Calcul prealable de la direction pour chaque noeud 519 : IF (IRADIA.NE.0) THEN 520 : DO IB=1,NBPOIN 521 : j=(NUM(1,IB)-1)*idimp1 522 : DO i=1,IDIM 523 : XNOR(i)=XCOOR(j+i)-U1(i) 524 : ENDDO 525 : IF (IDIM.EQ.2) THEN 526 : YL=XNOR(1)*XNOR(1)+XNOR(2)*XNOR(2) 527 : IF (YL.LT.EPSI) THEN 527 : CALL ERREUR(238) 528 : RETURN 529 : ENDIF 530 : YL=1.D0/SQRT(YL) 531 : IF (IRADIA.EQ.1) THEN 532 : XNOR(1)=XNOR(1)*YL 533 : XNOR(2)=XNOR(2)*YL 534 : ELSE IF (IRADIA.EQ.2) THEN 535 : XX=XNOR(1) 536 : XNOR(1)=-XNOR(2)*YL 537 : XNOR(2)=XX*YL 538 : ENDIF 539 : ELSE 540 : YL=XNOR(1)*U2(1)+XNOR(2)*U2(2)+XNOR(3)*U2(3) 541 : XL=0.D0 542 : DO i=1,3 543 : XNOR(i)=XNOR(i)-YL*U2(i) 544 : XL=XL+XNOR(i)*XNOR(i) 545 : ENDDO 546 : IF (XL.LT.EPSI) THEN 546 : CALL ERREUR(238) 547 : RETURN 548 : ENDIF 549 : IF (IRADIA.EQ.1) THEN 550 : XL=1.D0/SQRT(XL) 551 : XNOR(1)=XNOR(1)*XL 552 : XNOR(2)=XNOR(2)*XL 553 : XNOR(3)=XNOR(3)*XL 554 : ELSE IF (IRADIA.EQ.2) THEN 555 : XX=XNOR(1) 556 : YY=XNOR(2) 557 : ZZ=XNOR(3) 558 : XNOR(1)=YY*U2(3)-ZZ*U2(2) 559 : XNOR(2)=ZZ*U2(1)-XX*U2(3) 560 : XNOR(3)=XX*U2(2)-YY*U2(1) 561 : ENDIF 562 : ENDIF 563 : C XNOR contient la direction normee 564 : * SEGINI,XMATRI 565 : * IMATTT(IB)=XMATRI 566 : RE(1,1,IB)=0.D0 567 : RE(2,1,IB)=XNOR(1) 568 : RE(3,1,IB)=XNOR(2) 569 : RE(1,2,IB)=RE(2,1,IB) 570 : RE(1,3,IB)=RE(3,1,IB) 571 : IF (IDIM.EQ.3) THEN 572 : RE(4,1,IB)=XNOR(3) 573 : RE(1,4,IB)=RE(4,1,IB) 574 : ENDIF 575 : * SEGDES,XMATRI 576 : ENDDO 577 : SEGDES,XMATRI 578 : C Option DIRECTION 579 : ELSE IF (IDIREC.EQ.1) THEN 580 : IF (IDICHP.EQ.1) THEN 581 : DO I=1,NBPOIN 582 : * SEGINI,XMATRI 583 : RE(1,1,i)=0.D0 584 : RE(2,1,i)=BB(1,I) 585 : RE(3,1,i)=BB(2,I) 586 : RE(1,2,i)=RE(2,1,i) 587 : RE(1,3,i)=RE(3,1,i) 588 : RE(2,3,i)=RE(3,2,i) 589 : IF (IDIM.EQ.3) THEN 590 : RE(4,1,i)=BB(3,I) 591 : RE(1,4,i)=RE(4,1,i) 592 : ENDIF 593 : ** SEGDES,XMATRI 594 : * IMATTT(I)=XMATRI 595 : ENDDO 596 : * SEGDES,XMATRI 597 : ELSE 598 : * SEGINI,XMATRI 599 : DO i=1,NBPOIN 600 : RE(1,1,i)=0.D0 601 : RE(2,1,i)=XNOR(1) 602 : RE(3,1,i)=XNOR(2) 603 : RE(1,2,i)=RE(2,1,i) 604 : RE(1,3,i)=RE(3,1,i) 605 : RE(2,3,i)=RE(3,2,i) 606 : IF (IDIM.EQ.3) THEN 607 : RE(4,1,i)=XNOR(3) 608 : RE(1,4,i)=RE(4,1,i) 609 : ENDIF 610 : * DO i=1,NBPOIN 611 : * IMATTT(i)=XMATRI 612 : ENDDO 613 : * SEGDES,XMATRI 614 : SEGDES,xMATRI 615 : ENDIF 616 : C Autres options : 617 : ELSE 618 : * SEGINI,XMATRI 619 : * IXMATR=XMATRI 620 : DO i=1,NBPOIN 621 : RE(1,1,i)=0.D0 622 : RE(2,1,i)=1.D0 623 : RE(2,2,i)=0.D0 624 : RE(1,2,i)=RE(2,1,i) 625 : SEGDES,XMATRI 626 : * DO i=1,NBPOIN 627 : * IMATTT(i)=IXMATR 628 : ENDDO 629 : SEGDES,xMATRI 630 : SEGDES,IPT1 631 : ENDIF 632 : ENDDO 633 : C Fin de la boucle sur les IAA DDLs a bloquer 634 : IF (IDICHP.EQ.1) THEN 635 : * WRITE(IOIMP,*) 'Destruction du segment de travail' 636 : SEGSUP,MTRAV 637 : ENDIF 638 : C Fin normale du traitement 639 : SEGDES,MRIGID 639 : CALL ECROBJ('RIGIDITE',KRIGI) 640 : IF (MELEME.GT.0) THEN 641 : IF (KOBJET.EQ.0) THEN 642 : SEGSUP,MELEME 643 : ELSE 644 : IF (MELEME.EQ.KOBJET) THEN 645 : SEGDES,MELEME 646 : ELSE 647 : ** SEGSUP,MELEME 648 : SEGDES,MELEME 649 : ENDIF 650 : ENDIF 651 : ENDIF 652 : C Il peut rester des maillages actifs 653 : MRIGID=KRIGI 654 : SEGACT,MRIGID 655 : DO i=1,IRIGEL(/2) 656 : MELEME=IRIGEL(1,i) 657 : IF (MELEME.GT.0) SEGDES,MELEME 658 : ENDDO 659 : SEGDES,MRIGID 660 : 661 : C Sortie du sousprogramme (normale ou en cas d'erreur) 662 : 1000 SEGSUP,MSWBLO 663 : SEGDES,MCOORD 664 : SEGACT,MCOORD 665 : 666 : RETURN 667 : END 668 : 669 : 670 : 671 :
© Cast3M 2003 - All rights reserved.
Disclaimer