Download rela2.eso

Back to the list

   1 : C RELA2     SOURCE    BP208322  13/12/20    21:15:11     7890           
   2 :       SUBROUTINE RELA2(LEMOT)
   3 : ************************************************************************
   4 : *   cet operateur impose des mouvements de corps rigide 
   5 : *   ou l accrochage de deux maillages 
   6 : *
   7 : *   syntaxe : 
   8 : *          rig1 = rela   |   cori   depl    (rota)      |
   9 : *                        |                              |  geo1 (geo2)
  10 : *                        | ense ux uy uz ut rx ry rz rt |
  11 : *
  12 : *          cori  = mot cle
  13 : *          ense  = mot cle
  14 : *          depl  = pour imposer que la distance entre les points
  15 : *                  reste constante
  16 : *          rota  = pour imposer que tous les points ont les memes
  17 : *                  rotations
  18 : *          ux    = pour imposer que la valeur sur la composante ux
  19 : *                  est la meme pour tous les points
  20 : *          geo1  = objet de type meleme ou point
  21 : *          geo2  = objet de type meleme
  22 : *          rig1  = objet de type rigidite
  23 : *
  24 : *   accrochage syntaxe : 
  25 : *          rig1 = RELA  geo1 ACCR  geo2  ; ( geo2 massif ) 
  26 : *
  27 : ************************************************************************
  28 : 
  29 :       IMPLICIT INTEGER(I-N)
  30 :       IMPLICIT REAL*8 (A-H,O-Z)
  31 : -INC CCOPTIO
  32 : -INC CCGEOME
  33 : -INC SMELEME
  34 : -INC SMCOORD
  35 : -INC SMRIGID
  36 : -INC CCHAMP
  37 : c
  38 :       SEGMENT MOTDDV
  39 :        CHARACTER*4 MOTDDL(0)
  40 :       ENDSEGMENT
  41 :       SEGMENT NOMINV
  42 :        CHARACTER*4 NOMINC(0)
  43 :       ENDSEGMENT
  44 :       SEGMENT NUMAIT(NBMAIT)
  45 :       SEGMENT NUESCL(NBESCL)
  46 : c
  47 :       CHARACTER*4 ITCORI(3),LEMOT
  48 :       CHARACTER*4 MOACCR(2)
  49 :       CHARACTER*4 MODEPL(6),MODEDU(6)
  50 :       CHARACTER*8 CTYP
  51 :       DIMENSION COEF(3)
  52 : c
  53 :       DATA ITCORI(1)/'DEPL'/
  54 :       DATA ITCORI(2)/'ROTA'/
  55 :       DATA ITCORI(3)/'NOVE'/
  56 :       DATA MODEPL(1)/'UX  '/
  57 :       DATA MODEPL(2)/'UY  '/
  58 :       DATA MODEPL(3)/'UZ  '/
  59 :       DATA MODEPL(4)/'UR  '/
  60 :       DATA MODEPL(5)/'UZ  '/
  61 :       DATA MODEPL(6)/'UT  '/
  62 :       DATA MODEDU(1)/'FX  '/
  63 :       DATA MODEDU(2)/'FY  '/
  64 :       DATA MODEDU(3)/'FZ  '/
  65 :       DATA MODEDU(4)/'FR  '/
  66 :       DATA MODEDU(5)/'FZ  '/
  67 :       DATA MODEDU(6)/'FT  '/
  68 :       DATA MOACCR(1)/'FORT'/
  69 :       DATA MOACCR(2)/'FAIB'/
  70 : c
  71 :       DATA ZERO,UNO   /0.D0,1.D0/
  72 : 
  73 : 
  74 : * on teste la compatibilité des données ******************************
  75 : * en fourier le mouvement de corps rigide n'a pas toujours de sens
  76 :       IF ( LEMOT.NE.'ENSE') THEN
  77 : * cette remarque ne vaut que pour l'option cori
  78 :        IF (((IDIM.EQ.3).AND.(IFOMOD.NE.2)).OR.
  79 :      +   ((IDIM.EQ.2).AND. (IFOMOD.EQ.2)).OR.
  80 :      +   ((IFOUR.EQ.1).AND.( ABS(NIFOUR).GT.1))) THEN
  80 :          CALL ERREUR(21)
  81 :          RETURN
  82 :        ENDIF
  83 :       ENDIF
  84 : c
  85 : c on cherche a lire les directions ***********************************
  86 :       SEGINI MOTDDV,NOMINV
  87 :       IDEPL=0
  88 :       IROTA=0
  89 :       IENSE=0
  90 :       JVER=1
  91 : 
  92 : c
  93 :       IF(IFOUR.NE.-2.AND.IFOUR.NE.-1.AND.IFOUR.NE.-3)GOTO 101
  94 : c
  95 : c     deformations planes ou contraintes planes ou def. planes generalisées
  96 :       LDEPL=2
  97 :       IADEPL=0
  98 :       IAROTA=2
  99 :       GOTO 107
 100 :       
 101 :   101 CONTINUE
 102 :       IF(IFOUR.NE.1) GOTO 103
 103 : c
 104 : c     fourier
 105 :       LDEPL=3
 106 :       IADEPL=3
 107 :       IAROTA=2
 108 :       GOTO 107
 109 :       
 110 :   103 CONTINUE
 111 :       IF(IFOUR.NE.2) GOTO 107
 112 : c
 113 : c     tridim
 114 :       LDEPL=3
 115 :       IADEPL=0
 116 :       IAROTA=0
 117 : 
 118 :   107 CONTINUE
 119 :   
 120 : C option ACCRO ou GLIS *********************************************
 121 :       IF(LEMOT.EQ.'ACCR'.OR.LEMOT.EQ.'GLIS') THEN
 122 :       igli=0
 123 :       iaccr=0
 124 :       if (LEMOT.EQ.'GLIS' ) then
 125 :         igli=1
 126 :       else
 127 : *       ACCRO : lecture facultative du mot clé FORT ou FAIBLE 
 127 :         CALL LIRMOT(MOACCR,2,iaccr,0)
 128 :       endif
 129 :       if (iaccr.eq.0) then
 130 : c       RELA GLIS ou RELA ACCRO => ACCRO
 130 :         CALL ACCRO (igli)
 131 :       elseif (iaccr.eq.1) then
 132 : c       RELA ACCRO 'FORT' => ACCRO3
 132 :         CALL ACCRO3
 133 :       elseif (iaccr.eq.2) then
 134 : c       RELA ACCRO 'FAIBL' => ACCRO2     
 134 :         CALL ACCRO2    
 135 :       else
 136 : c       Option indisponible
 136 :         call erreur(19)
 137 :       endif
 138 :       RETURN
 139 :       ENDIF    
 140 :       
 141 : c option ENSE ******************************************************
 142 :       IF(LEMOT.EQ.'ENSE') THEN
 143 :       
 144 :       IENSE = 1
 145 :   445 CONTINUE
 145 :       CALL LIRMOT(NOMDD,LNOMDD,IMOT,0)
 146 :       IF(IMOT.EQ.0) GO TO 446
 147 :       MOTDDL(**)=NOMDD(IMOT)
 148 :       MOTDDL(**)=NOMDU(IMOT)
 149 :       GO TO 445
 150 :   446 CONTINUE
 151 :       IBDDL=MOTDDL(/2)
 152 :       
 153 : c option CORI ******************************************************
 154 :       ELSE IF(LEMOT.EQ.'CORI') THEN
 155 : 
 156 : cccccccc   a t on une directive depl
 157 : c
 158 :       IOBL=1
 159 :  4480 CONTINUE
 159 :       CALL LIRMOT(ITCORI,3,IMOT,IOBL)
 160 :       IF(IERR.NE.0) RETURN
 161 :       IF(IMOT.EQ.0) GOTO 448
 162 :       IOBL=0
 163 :       GO TO (44801,44802,44803),IMOT
 164 : 44801 continue
 165 : c
 166 : cccccccccccc on a trouve le mot deplacement
 167 : c
 168 :       IDEPL=1
 169 :       DO 4481 IA=1,LDEPL
 170 :       NOMINC(**)=MODEPL(IADEPL+IA)
 171 :       NOMINC(**)=MODEDU(IADEPL+IA)
 172 :  4481 CONTINUE
 173 :       GO TO 4480
 174 : 44802 continue
 175 : c
 176 : cccccccc      on a trouve le mot rotation
 177 : c
 178 :       IROTA=1
 179 :       GO TO 4480
 180 : 44803 continue
 181 :       JVER=0
 182 :       GO TO 4480
 183 :   448 CONTINUE
 184 :   
 185 :       ENDIF
 186 : c endif des options ENSE et CORI ***********************************
 187 : 
 188 : c
 189 : c on cherche a lire le(s) maillage(s) ******************************
 190 : c
 191 : * cas de mot clé rota ou fourier 0 ou axis 0
 192 :       IF ((IDEPL.EQ.1).AND.((IROTA.EQ.1).OR.
 193 :      +   (NIFOUR.EQ.0 .AND. IFOMOD.EQ.1).OR.(IFOMOD.EQ.0))) THEN
 193 :       CALL QUETYP(CTYP,1,IRETOU)
 194 :       IF (IRETOU.EQ.0) GOTO 555
 195 :         IF(CTYP.EQ.'MAILLAGE') THEN
 195 :          CALL LIROBJ('MAILLAGE',KOBJET,1,IRETOU)
 196 :          IF (IRETOU.EQ.0) GOTO 555
 197 :          LP=0
 198 : * on regarde si on a donné juste un point
 199 :         ELSE IF (CTYP.EQ.'POINT   ') THEN
 199 :          CALL LIROBJ('POINT   ',KOBJET,1,IRETOU)
 200 :          IF (IRETOU.EQ.0) GOTO 555
 201 :          LP=1
 202 :         ELSE
 203 :          GOTO 555
 204 :         ENDIF
 204 :        CALL LIROBJ('MAILLAGE',KOBJE2,0,IRETOU)
 205 :        IF((IRETOU.EQ.0).AND.(LP.EQ.1)) THEN
 205 :          CALL ERREUR(26)
 206 :          GOTO 555
 207 :        ENDIF
 207 :        CALL RELA3 (KOBJET,KOBJE2,LP,IROTA)
 208 :        GOTO 555
 209 :       ELSE
 210 : * cas sans mot clé rota (sauf  fourier 0 et axis 0)
 210 :         CALL LIROBJ('MAILLAGE',KOBJET,1,IRETOU)
 211 :         IF(IRETOU.EQ .0) GOTO 555
 212 :       ENDIF
 213 : 
 214 :       MELEME=KOBJET
 215 :       SEGACT MELEME
 215 :       IF(ITYPEL .NE.1) CALL CHANGE(MELEME,1)
 216 :       NBPOIN=NUM(/2)
 217 : c
 218 : c  dans le cas depl on ecrit que la distance entre les points reste
 219 : c  constante ,1 matrice avec
 220 : c  en 2d , 3 noeuds maitres et 2(n-3)+3 relations
 221 : c  en 3d , 4 noeuds maitres et 3(n-4)+6 relations
 222 : c   on prend comme noeuds maitres possibles
 223 : c       les 3 premiers noeuds non alignes en 2d
 224 : c       les 4 premiers noeuds non coplanaires en 3d
 225 : c
 226 :       IF (IDEPL.EQ.1) THEN
 227 :         IF (IDIM.EQ.2) NBMAIT=3
 228 :         IF (IDIM.EQ.3) NBMAIT=4
 229 : c   erreur nb de pts du maillage inferieur au nb de noeuds maitres
 230 :         IF (NBPOIN.LT.NBMAIT.AND.JVER.EQ.1) THEN
 230 :           CALL ERREUR(364)
 231 :           GOTO 998
 232 :         ENDIF
 233 :         IF(JVER.EQ.1) THEN
 234 :           NBESCL=NBPOIN-NBMAIT
 235 :           SEGINI NUMAIT,NUESCL
 236 : cc   cas 2d determination des 3 noeuds maitres
 237 :           IF (IDIM.EQ.2) THEN
 238 :             IP1=NUM(1,1)
 239 :             IP2=NUM(1,2)
 240 :             IREFP1=(IP1-1)*(IDIM+1)
 241 :             IREFP2=(IP2-1)*(IDIM+1)
 242 :             X12=XCOOR(IREFP2+1)-XCOOR(IREFP1+1)
 243 :             Y12=XCOOR(IREFP2+2)-XCOOR(IREFP1+2)
 244 :             P12N=SQRT(X12**2+Y12**2)
 245 :             DO 300 I=3,NBPOIN
 246 :               IP3=NUM(1,I)
 247 :               IREFP3=(IP3-1)*(IDIM+1)
 248 :               X13=XCOOR(IREFP3+1)-XCOOR(IREFP1+1)
 249 :               Y13=XCOOR(IREFP3+2)-XCOOR(IREFP1+2)
 250 :               P13N=SQRT(X13**2+Y13**2)
 251 :               TEST=X12*X13+Y12*Y13
 252 :               IF(ABS(TEST).GT.0.9997D0*(P12N*P13N)) GOTO 300
 253 :               GOTO 399
 254 : 300         continue
 255 : 310         continue
 256 : cccc   erreur  tous les points sont alignes ou quasiment
 256 :             CALL ERREUR(362)
 257 :             GOTO 999
 258 : cc   cas 3d determination des 4 noeuds maitres
 259 :           ELSE IF (IDIM.EQ.3) THEN
 260 :             IP1=NUM(1,1)
 261 :             IP2=NUM(1,2)
 262 :             IREFP1=(IP1-1)*(IDIM+1)
 263 :             IREFP2=(IP2-1)*(IDIM+1)
 264 :             X12=XCOOR(IREFP2+1)-XCOOR(IREFP1+1)
 265 :             Y12=XCOOR(IREFP2+2)-XCOOR(IREFP1+2)
 266 :             Z12=XCOOR(IREFP2+3)-XCOOR(IREFP1+3)
 267 :             P12N=SQRT(X12**2+Y12**2+Z12**2)
 268 :             DO 350 I=3,NBPOIN
 269 :               IP3=NUM(1,I)
 270 :               IREFP3=(IP3-1)*(IDIM+1)
 271 :               X13=XCOOR(IREFP3+1)-XCOOR(IREFP1+1)
 272 :               Y13=XCOOR(IREFP3+2)-XCOOR(IREFP1+2)
 273 :               Z13=XCOOR(IREFP3+3)-XCOOR(IREFP1+3)
 274 :               P13N=SQRT(X13**2+Y13**2+Z13**2)
 275 :               TEST=X12*X13+Y12*Y13+Z12*Z13
 276 :               IF(ABS(TEST).GT.0.9997D0*(P12N*P13N)) GOTO 350
 277 :               GOTO 360
 278 : 350         continue
 279 :             GOTO 380
 280 : 360         continue
 281 :             XN=Y12*Z13-Z12*Y13
 282 :             YN=Z12*X13-X12*Z13
 283 :             ZN=X12*Y13-Y12*X13
 284 :             PN=XN**2+YN**2+ZN**2
 285 :             DO 370 I=3,NBPOIN
 286 :               IP4=NUM(1,I)
 287 :               IREFP4=(IP4-1)*(IDIM+1)
 288 :               X14=XCOOR(IREFP4+1)-XCOOR(IREFP1+1)
 289 :               Y14=XCOOR(IREFP4+2)-XCOOR(IREFP1+2)
 290 :               Z14=XCOOR(IREFP4+3)-XCOOR(IREFP1+3)
 291 :               P14N=SQRT(X14**2+Y14**2+Z14**2)
 292 :               TEST=XN*X14+YN*Y14+ZN*Z14
 293 :               IF(ABS(TEST).LT.0.025D0*(P12N*P13N)) GOTO 370
 294 :               GOTO 399
 295 : 370         continue
 296 : 380         continue
 297 : cccc   erreur  tous les points sont coplanaires ou quasiment
 297 :             CALL ERREUR(363)
 298 :             GOTO 999
 299 :           ENDIF
 300 : 399       continue
 301 :           NUMAIT(1)=IP1
 302 :           NUMAIT(2)=IP2
 303 :           NUMAIT(3)=IP3
 304 :           IF (IDIM.EQ.3) NUMAIT(4)=IP4
 305 :           I1=0
 306 :           DO 397 I=3,NBPOIN
 307 :             DO 398 J=1,NBMAIT
 308 :               IF(NUM(1,I).EQ.NUMAIT(J)) GOTO 397
 309 : 398         continue
 310 :             I1=I1+1
 311 :             NUESCL(I1)=NUM(1,I)
 312 : 397       continue
 313 :         ELSEIF(JVER.EQ.0) THEN
 314 :           NBESCL=NBPOIN-NBMAIT
 315 :           IF(NBESCL.GT.0) THEN
 316 :             SEGINI NUMAIT,NUESCL
 317 :           ELSE
 318 :             NBMAIT=NBPOIN
 319 :             NBESCL=0
 320 :             SEGINI NUMAIT
 321 :           ENDIF
 322 :           DO 396 I=1,NBMAIT
 323 :               NUMAIT(I)=NUM(1,I)
 324 :   396     CONTINUE
 325 :           IF(NBESCL.GT.0) THEN
 326 :              DO 395 I=1,NBESCL
 327 :                 NUESCL(I)=NUM(1,I+NBMAIT)
 328 :   395        CONTINUE
 329 :           ENDIF
 330 :         ENDIF
 331 :       ENDIF
 332 : c
 333 : c  dans le cas ense ,on ecrit que la valeur sur un ddl
 334 : c  donne est la meme pour tous les points ,ibddl/2 matrices avec
 335 : c  (n-1) relations
 336 : c
 337 :       IF(IENSE.EQ.1) NNMAT=IBDDL/2
 338 :       IF(IDEPL.EQ.1) THEN
 339 :        NNMAT=1
 340 : * en fourier 1 on a en plus que ur=-ut a tout point et une relation d'analogie
 341 : * sur uz entre deux noeuds (p.ex. 2 noeuds maitres)
 342 :        IF(ABS(NIFOUR).EQ.1) NNMAT=3
 343 :       ENDIF
 344 :       IF(NNMAT.EQ.0) GO TO 555
 345 :       NRIGE=7
 346 :       NRIGEL=NNMAT
 347 :       SEGINI MRIGID
 348 :       ICHOLE=0
 349 :       IMGEO1=0
 350 :       IMGEO2=0
 351 :       ISUPEQ=0
 352 :       IFORIG=IFOMOD
 353 :       MTYMAT='RIGIDITE'
 354 :       KRIGI=MRIGID
 355 : c
 356 : cccc on rajoute les noeuds associes aux multiplicateurs au maillage
 357 : c
 358 : cc      segact mcoord
 359 : ccccccccc   nombre de noeuds
 360 :       NBNO=XCOOR(/1)/(IDIM+1)
 361 :       NBNOI=NBNO
 362 :       IF(IENSE.EQ.1) NBPTS=NBNO+NNMAT*(NBPOIN-1)
 363 :       IF(IDEPL.EQ.1) THEN
 364 :          IF(IDIM.EQ.2) NBPTS=NBNO+(2*(NBPOIN-3)+3)
 365 :          IF(IDIM.EQ.3) NBPTS=NBNO+(3*(NBPOIN-4)+6)
 366 :          IF(NBPOIN.EQ.2) NBPTS=NBNO +1
 367 :          IF(ABS(NIFOUR).EQ.1)NBPTS = NBPTS + NBPOIN + 1
 368 :       ENDIF
 369 :       SEGADJ MCOORD
 370 : c
 371 : cccc  boucle sur les rigidites elementaires
 372 : c   si idepl=1 et irota=1 la 1ere rigidite elementaire est associee a
 373 : c   l option depl et les suivantes a l option rota
 374 : c
 375 : 
 376 :       IF (IENSE.EQ.1) NBRELA=NBPOIN-1
 377 :       I7=0
 378 :       DO 7 IAA=1,NNMAT
 379 :       IF (IDEPL.EQ.1) THEN
 380 :         IF (IAA.EQ.1) THEN
 381 :           IF(IDIM.EQ.2) NBRELA=2*(NBPOIN-3)+3
 382 :           IF(IDIM.EQ.3) NBRELA=3*(NBPOIN-4)+6
 383 :           IF(NBPOIN.EQ.2) NBRELA=1
 384 :         ELSE IF (IAA.eq.2) THEN
 385 : * fourier 1 ur=-ut
 386 :           NBRELA=NBPOIN
 387 : * fourier 1 relation d'analogie sur uz
 388 :         ELSE IF (IAA.eq.3) THEN
 389 :           nbrela = 1
 390 :         ENDIF
 391 :       ENDIF
 392 : *      print*,'idepl=',idepl
 393 : *      print*,'iaa=',iaa
 394 : cccccccc on cree le noeud nbno+1 et nbno+2
 395 : cccccccc on les met a l origine
 396 :       DO 5 IA=1,NBRELA
 397 :       XCOOR(NBNOI*(IDIM+1)+1)=0.D0
 398 :       XCOOR(NBNOI*(IDIM+1)+2)=0.D0
 399 :       IF (IDIM.EQ.3) XCOOR(NBNOI*(IDIM+1)+3)=0.D0
 400 :       XCOOR(NBNOI*(IDIM+1)+(IDIM+1))=0.D0
 401 :       NBNOI=NBNOI+1
 402 :     5 CONTINUE
 403 : c
 404 : c     on initialise le segment meleme associe aux blocages
 405 : c
 406 :       NBSOUS=0
 407 :       NBREF=0
 408 :       NBNN=3
 409 :       NBELEM=NBRELA
 410 :       SEGINI IPT1
 411 :       IRIGEL(1,IAA)=IPT1
 412 :       IPT1.ITYPEL=22
 413 :       DO 400 I4=1,NBRELA
 414 :       IPT1.ICOLOR(I4)=IDCOUL
 415 : c punti associati ai moltiplicatori
 416 :       IPT1.NUM(1,I4)=NBNO+I7+1
 417 :       I7=I7+1
 418 : 400   continue
 419 : c caricamento nodi pseudo-elementi
 420 : * option ensemble
 421 :       IF (IENSE.EQ.1) THEN
 422 :         DO 420 I9=1,NBRELA
 423 :         IPT1.NUM(2,I9)=NUM(1,I9)
 424 :         IPT1.NUM(3,I9)=NUM(1,I9+1)
 425 : 420     continue
 426 : * fourier 1 ur = -ut
 427 :       ELSE IF (IDEPL.EQ.1.AND.IAA.EQ.2) THEN
 428 :         DO 421 I9=1,NBRELA
 429 :         IPT1.NUM(2,I9)=NUM(1,I9)
 430 :         IPT1.NUM(3,I9)=NUM(1,I9)
 431 : 421     continue
 432 :       ELSE IF (IDEPL.EQ.1.AND.IAA.EQ.1) THEN
 433 :         I9=0
 434 :         IF(NBPOIN.GT.NBMAIT) THEN
 435 : c pseudo-elements pour les relations des autres noeuds avec
 436 : c         les noeuds maitres
 437 :           DO 430 IC=1, NBESCL
 438 :           ICOREC = 0
 439 :           DO 431 IB=1,NBMAIT-1
 440 :             I9=I9+1
 441 :             IPT1.NUM(3,I9)=NUESCL(IC)
 442 :             IPT1.NUM(2,I9)=NUMAIT(IB)
 443 : * on teste si le noeud esclave est colineaire avec les 2 maitres
 444 :          IF (IB.EQ.2) THEN
 445 :             IREFES = (NUESCL(IC)-1)*(IDIM+1)
 446 :             XX13=XCOOR(IREFES+1)-XCOOR(IREFP1+1)
 447 :             YY13=XCOOR(IREFES+2)-XCOOR(IREFP1+2)
 448 :             PP13N=SQRT(XX13**2+YY13**2)
 449 :             TTEST=X12*XX13+Y12*YY13
 450 :              IF (IDIM.EQ.3) THEN
 451 :                ZZ13 = XCOOR(IREFES+3)-XCOOR(IREFP1+3)
 452 :                PP13N=PP13N+ZZ13**2
 453 :                TTEST=TTEST+Z12*ZZ13
 454 :               ENDIF
 455 :              IF(ABS(TTEST).GT.0.9997D0*(P12N*PP13N)) THEN
 456 :               IPT1.NUM(2,I9)=NUMAIT(IB+1)
 457 :               ICOREC = 1
 458 :              ENDIF
 459 :           ELSE IF (IB.EQ.3) THEN
 460 :            IF (ICOREC.EQ.1) THEN
 461 :                 IPT1.NUM(2,I9)=NUMAIT(IB+1)
 462 : * on teste si le noeud esclave est coplanaire avec les 3 maitres
 463 :            ELSE
 464 :             TTEST=XN*XX13+YN*YY13+ZN*ZZ13
 465 :        IF(ABS(TTEST).LT.0.025D0*(P12N*P13N))IPT1.NUM(2,I9)=NUMAIT(IB+1)
 466 :            ENDIF
 467 :           ENDIF
 468 : 431       continue
 469 : 430       continue
 470 :         ENDIF
 471 : c pseudo-elements pour les relations entre les noeuds maitres
 472 :         DO 432 IB=2,NBMAIT
 473 :           I9=I9+1
 474 :           IPT1.NUM(2,I9)=NUMAIT(1)
 475 :           IPT1.NUM(3,I9)=NUMAIT(IB)
 476 : 432     continue
 477 :         IF(NBMAIT.GE.3) THEN
 478 :           DO 433 IB=3,NBMAIT
 479 :              I9=I9+1
 480 :              IPT1.NUM(2,I9)=NUMAIT(2)
 481 :              IPT1.NUM(3,I9)=NUMAIT(IB)
 482 : 433       continue
 483 :         ENDIF
 484 :         IF (NBMAIT.GE.4) THEN
 485 :           DO 434 IB=4,NBMAIT
 486 :              I9=I9+1
 487 :              IPT1.NUM(2,I9)=NUMAIT(3)
 488 :              IPT1.NUM(3,I9)=NUMAIT(IB)
 489 : 434       continue
 490 :         ENDIF
 491 :       ENDIF
 492 : * fourier 1, relation d'analogie sur uz
 493 : * on cherche les 2 noeuds maitres avec xmax et xmin
 494 : * pour etre sur de ne pas avoir choisi 2 noeuds ayan la 
 495 : * meme coordonnee x
 496 :        if ((idepl.eq.1).and.(iaa.eq.3)) then
 497 :           xx = xcoor(irefp1+1)
 498 :           x1 = xx
 499 :           ipp1 = ip1
 500 :          if (xcoor(irefp2+1).gt.xx) then
 501 :            x1 =xcoor(irefp2+1)
 502 :            ipp1 = ip2
 503 :          else if (xcoor(irefp3+1).gt.xx) then 
 504 :            x1 =xcoor(irefp3+1)
 505 :            ipp1 = ip3
 506 :          endif
 507 :           x2 = xx
 508 :           ipp2 =ip1
 509 :          if (xcoor(irefp2+1).lt.xx) then
 510 :            x2 =xcoor(irefp2+1)
 511 :            ipp2 = ip2
 512 :          else if (xcoor(irefp3+1).lt.xx) then 
 513 :            x2 =xcoor(irefp3+1)
 514 :            ipp2 = ip3
 515 :          endif
 516 :         ipt1.num(2,1)=ipp2
 517 :         ipt1.num(3,1)=ipp1
 518 :        endif
 519 : *  correction on place les pts mult au barycentre des noeuds associes
 520 : *   a chaque relation
 521 :         DO 425 IA=1,NBRELA
 522 :         IREF1=(IPT1.NUM(1,IA)-1)*(IDIM+1)
 523 :         IREF2=(IPT1.NUM(2,IA)-1)*(IDIM+1)
 524 :         IREF3=(IPT1.NUM(3,IA)-1)*(IDIM+1)
 525 :         XCOOR(IREF1+1)=(XCOOR(IREF2+1)+XCOOR(IREF3+1))/2
 526 :         XCOOR(IREF1+2)=(XCOOR(IREF2+2)+XCOOR(IREF3+2))/2
 527 :         IF(IDIM.EQ.3) XCOOR(IREF1+3)=(XCOOR(IREF2+3)+XCOOR(IREF3+3))/2
 528 :  425    CONTINUE
 529 : cccccccc   on vient de crer les segment meleme resultat
 530 : cccccccc   on va creer les raideurs
 531 :       IRIGEL(2,IAA)=0
 532 :       IRIGEL(5,IAA)=NIFOUR
 533 :       IRIGEL(6,IAA)=0
 534 :       IRIGEL(7,IAA)=0
 535 :       IF(IENSE.EQ.1.OR.IAA.GT.1) NLIGRE=3
 536 :       IF(IDEPL.EQ.1.AND.IAA.EQ.1) NLIGRE=2*IDIM+1
 537 :       NLIGRP=NLIGRE
 538 :       NLIGRD=NLIGRE
 539 :       SEGINI DESCR
 540 :       IRIGEL(3,IAA)=DESCR
 541 : c         on remplit le tableau des descripteurs de rig
 542 : 
 543 :         LISINC(1)='LX'
 544 :         LISDUA(1)='FLX'
 545 :         NOELEP(1)=1
 546 :         NOELED(1)=1
 547 :       IF(IENSE.EQ.1) THEN
 548 :         JAA=2*(IAA-1)
 549 :         LISINC(2)=MOTDDL(JAA+1)
 550 :         LISINC(3)=LISINC(2)
 551 :         LISDUA(2)=MOTDDL(JAA+2)
 552 :         LISDUA(3)=LISDUA(2)
 553 :         NOELEP(2)=2
 554 :         NOELEP(3)=3
 555 :         NOELED(2)=2
 556 :         NOELED(3)=3
 557 :       ELSE IF (IDEPL.EQ.1.AND.IAA.EQ.1) THEN
 558 :         DO 250 I=1,IDIM
 559 :         JAA=2*(I-1)
 560 :         J=1+2*(I-1)
 561 :         LISINC(J+1)=NOMINC(JAA+1)
 562 :         LISINC(J+2)=LISINC(J+1)
 563 :         LISDUA(J+1)=NOMINC(JAA+2)
 564 :         LISDUA(J+2)=LISDUA(J+1)
 565 :         NOELEP(J+1)=2
 566 :         NOELEP(J+2)=3
 567 :         NOELED(J+1)=2
 568 :         NOELED(J+2)=3
 569 : 250   continue
 570 : 
 571 : * fourier 1 ur = -ut
 572 :       ELSE IF (IDEPL.EQ.1.AND.IAA.EQ.2) THEN
 573 :         LISINC(2)=MODEPL(4)
 574 :         LISINC(3)=MODEPL(6)
 575 :         LISDUA(2)=MODEDU(4)
 576 :         LISDUA(3)=MODEDU(6)
 577 :         NOELEP(2)=2
 578 :         NOELEP(3)=3
 579 :         NOELED(2)=2
 580 :         NOELED(3)=3
 581 :       ELSE IF (IDEPL.EQ.1.AND.IAA.EQ.3) THEN
 582 :         LISINC(2)=MODEPL(3)
 583 :         LISINC(3)=MODEPL(3)
 584 :         LISDUA(2)=MODEDU(3)
 585 :         LISDUA(3)=MODEDU(3)
 586 :         NOELEP(2)=2
 587 :         NOELEP(3)=3
 588 :         NOELED(2)=2
 589 :         NOELED(3)=3
 590 :       ENDIF
 591 :       SEGDES DESCR
 592 : c
 593 :       NELRIG=NBRELA
 594 :       SEGINI xMATRI
 595 :       IRIGEL(4,IAA)=xMATRI
 596 :       COERIG(IAA)=1.D0
 597 : c
 598 : c  remplissage de  RE
 599 : 
 600 :       IF (IENSE.EQ.1.OR.IAA.GT.1) THEN
 601 :            do iou=1,nelrig
 602 : *      SEGINI XMATRI
 603 : *      IXMATR=XMATRI
 604 :             RE(1,1,iou)=  0.D0
 605 :             RE(2,1,iou)= 1.D0
 606 :          IF (IENSE.EQ.1) THEN
 607 :             RE(3,1,iou)= -1.D0
 608 : * ur = - ut
 609 :          ELSE IF (IAA.EQ.2) THEN
 610 :             RE(3,1,iou)= 1.D0
 611 : * analogie sur uz en fourier 1
 612 :          ELSE
 613 :             RE(3,1,iou)= -x2/x1
 614 :          ENDIF
 615 :               DO 777 I=1,1
 616 :                 DO 666 J=2,NLIGRE
 617 :                   RE(I,J,iou)=RE(J,I,iou)
 618 :  666             CONTINUE
 619 :  777           CONTINUE
 620 :        enddo
 621 : *      SEGDES XMATRI
 622 : *      DO 3 IA=1,NELRIG
 623 : *      IMATTT(IA)=IXMATR
 624 : *    3 CONTINUE
 625 :       ELSE IF (IDEPL.EQ.1.AND.IAA.EQ.1) THEN
 626 :       DO 740 IA=1,NELRIG
 627 :       DO 741 I=1,IDIM
 628 :         IREF2=(IPT1.NUM(2,IA)-1)*(IDIM+1)
 629 :         IREF3=(IPT1.NUM(3,IA)-1)*(IDIM+1)
 630 :       COEF(I)=XCOOR(IREF2+I)-XCOOR(IREF3+I)
 631 : 741   continue
 632 : *      SEGINI XMATRI
 633 : *      IMATTT(IA)=XMATRI
 634 :       RE(1,1,ia)= 0.D0
 635 : 
 636 :       DO 760 I1=2,NLIGRE
 637 :       I4=INT(I1/2)
 638 :       COEF1=COEF(I4)
 639 :       IF((2*I4).LT.I1) COEF1=-COEF(I4)
 640 :       DO 760 I2=1,1
 641 :       RE(I1,I2,ia)=COEF1
 642 :       RE(I2,I1,ia)=COEF1
 643 : 760   continue
 644 : *      SEGDES XMATRI
 645 : 740   continue
 646 :       ENDIF
 647 :       SEGDES xMATRI
 648 : c
 649 :       SEGDES IPT1
 650 :     7 CONTINUE
 651 : 
 652 : c
 653 :       SEGDES MRIGID
 653 :       CALL ECROBJ('RIGIDITE',KRIGI)
 654 : 999   continue
 655 :       IF (IDEPL.EQ.1) THEN
 656 :           SEGSUP NUMAIT
 657 :           IF(NBESCL.GT.0) SEGSUP NUESCL
 658 :       ENDIF
 659 : 998   continue
 660 :        SEGDES MELEME
 661 :   555 CONTINUE
 662 :       SEGSUP MOTDDV,NOMINV
 663 :       RETURN
 664 :       END
 665 :  
 666 :  
 667 :  
 668 :  
 669 :  

© Cast3M 2003 - All rights reserved.
Disclaimer