1 : C ACCRO SOURCE BP208322 12/12/17 21:15:06 7615 2 : SUBROUTINE ACCRO(igliss) 3 : C======================================================================== 4 : C ACCROCHAGE D UN MAILLAGE SUR UN MASSIF PAR L INTERMEDIAIRE 5 : C DES FONCTIONS DE FORME 6 : C SUBROUTINE APPELLEE PAR RELA2 7 : C JMB OCTOBRE 1996 8 : C======================================================================== 9 : 10 : IMPLICIT INTEGER(I-N) 11 : IMPLICIT REAL*8 (A-H,O-Z) 12 : -INC CCREEL 13 : -INC CCOPTIO 14 : -INC SMCOORD 15 : -INC SMELEME 16 : -INC SMRIGID 17 : -INC SMLMOTS 18 : -INC CCHAMP 19 : -INC CCGEOME 20 : -INC SMMODEL 21 : external shape 22 : DATA EPSI/1.D-9/ 23 : DIMENSION ICOR(100),IMEL(100),imtt(100) 24 : SEGMENT IEXCLU 25 : INTEGER IDEJVU(NBPTI) 26 : ENDSEGMENT 27 : SEGMENT ICOUNT 28 : INTEGER IEINT(2,NODES),MAXPZ(1) 29 : REAL*8 QQQ(3,NODES) 30 : ENDSEGMENT 31 : SEGMENT ICOSOM 32 : INTEGER IPTT(NBSZ) 33 : ENDSEGMENT 34 : SEGMENT MTRAV 35 : REAL*8 CORDEL(3,NBNN),QSI(3),XPO(3),SHPP(6,NBNN),COE(IDIM,30) 36 : REAL*8 CC(100),XEL(3,NBNN),BPSS(3,3),XPT(3),CORED(3,NBNN) 37 : ENDSEGMENT 38 : SEGMENT IOUBLI 39 : INTEGER IPTOU(NBOUB,2) 40 : ENDSEGMENT 41 : segment itange 42 : real*8 xperp1(idim,nbpti),xperp2(idim,nbpti) 43 : integer iextr(nbpti) 44 : endsegment 45 : segment inomil(0) 46 : segment icpr(xcoor(/1)/(idim+1)) 47 : segment imul(idim,nbpt) 48 : segment icomu(nbpt) 49 : logical ltelq 50 : 51 : 52 : C LECTURE des MAILLAGE et MMODEL *********************************** 53 : 54 : IFLAG= 0 55 : IFROCA=0 56 : MELEME=0 56 : CALL LIROBJ('MAILLAGE',IPP2,1,IRETOU) 57 : IF(IERR.NE.0) RETURN 57 : CALL LIROBJ('MAILLAGE',MELEME,0,IRETOU) 58 : MELEMI=MELEME 59 : IF(IRETOU.EQ.0) then 60 : melemi=ipp2 61 : meleme=ipp2 61 : CALL LIROBJ('MMODEL',IPMODL,1,IFROCA) 62 : IF (IERR.NE.0) RETURN 63 : endif 64 : c ! Erwan 65 : 66 : 67 : C**** CAS OU UN MODELE A ETE LU ************************************** 68 : if (IFROCA.ne.0) then 69 : mmode1 = ipmodl 70 : segact mmode1 71 : nbpt=0 72 : segini icpr 73 : C creation des l'ensemble des points a accrocher et des multiplicateurs 74 : C correspondant. imul(i,j) est le ieme multiplicateur a accrocher au j eme 75 : C point local ( icpr(ia)=j) 76 : C --- puis Assemblage du maillage global (tous les câbles) 77 : do 500 isouf = 1,mmode1.kmodel(/1) 78 : imode1 = mmode1.kmodel(isouf) 79 : segact imode1 80 : if(imode1.tymode(1).eq.'MMODEL') then 81 : mmode2 = imode1.ivamod(1) 82 : segact mmode2 83 : imode2 = mmode2.kmodel(1) 84 : else 84 : call erreur(975) 85 : return 86 : endif 87 : segact imode2 88 : ipt7 = imode2.imamod 89 : if(isouf.eq.1) then 90 : ipt5=ipt7 91 : else 92 : ltelq=.false. 92 : call fuse(ipp4,ipt7,ipt5,ltelq) 93 : ipt5=ipp4 94 : endif 95 : ipt6=imode1.imamod 96 : segact ipt6 97 : if(ipt6.itypel.ne.22) then 97 : call erreur (19) 98 : return 99 : endif 100 : segdes imode2,mmode2 101 : do 501 iel=1,ipt6.num(/2) 102 : ia=ipt6.num(2,iel) 103 : if(icpr(ia).eq.0 ) then 104 : nbpt=nbpt+1 105 : icpr(ia)=nbpt 106 : endif 107 : 501 continue 108 : 500 continue 109 : segini imul,icomu 110 : nbpt=0 111 : do 502 isouf = 1,mmode1.kmodel(/1) 112 : imode1 = mmode1.kmodel(isouf) 113 : ipt6=imode1.imamod 114 : do 503 iel=1,ipt6.num(/2) 115 : ia=ipt6.num(2,iel) 116 : ib=icpr(ia) 117 : icomu(ib)=icomu(ib)+1 118 : imul(icomu(ib),ib)=ipt6.num(1,iel) 119 : 503 continue 120 : segdes imode1,ipt6 121 : 502 continue 122 : segdes mmode1 123 : ipp2=ipt5 124 : * fin du if(ifroca..; 125 : endif 126 : C**** FIN DU PRE-TRAITEMENT DANS LE CAS D'UN MODELE ****************** 127 : 128 : 129 : C**** LECTURE DES NOMS D'INCONNUE ET DU CRITERE ********************* 130 : 130 : CALL LIROBJ('LISTMOTS',MLMOTX,0,IRETOU) 131 : CCRIT =1.D-5 131 : CALL LIRREE(CCRIT,0,IRETOU) 132 : CCRIT = -ABS(CCRIT) 133 : 134 : IF(MLMOTX.EQ.0) THEN 135 : C options par defaut 136 : JGN=4 137 : JGM=IDIM 138 : NINC= IDIM 139 : SEGINI MLMOTS 140 : IF(IFOUR.EQ.0) THEN 141 : C axisymetrique 142 : MOTS(1)='UR ' 143 : MOTS(2)='UZ ' 144 : ELSE 145 : MOTS(1)='UX ' 146 : MOTS(2)='UY ' 147 : IF (IDIM.EQ.3) MOTS(3)='UZ ' 148 : ENDIF 149 : ELSE 150 : C on a lu les noms d'inconnues 151 : if(igliss.eq.1) then 151 : call erreur (19) 152 : return 153 : endif 154 : MLMOTS = MLMOTX 155 : SEGACT MLMOTS 156 : NINC = MOTS(/2) 157 : IF(NINC.GT.100) THEN 158 : SEGDES MLMOTS 158 : CALL ERREUR(19) 159 : RETURN 160 : ENDIF 161 : ENDIF 162 : C IF (MLMOTX.EQ.0) THEN 163 : C MOTERR(1:8)='LISTMOTS' 163 : C CALL ERREUR(37) 164 : C RETURN 165 : C ENDIF 166 : 167 : 168 : C**** VERIFICATION DU MAILLAGE MELEME ******************************** 169 : 170 : C LE SECOND MAILLAGE DOIT ETRE MASSIF 171 : C POUR LE MOMENT TETRAEDRES ET PYRAMIDES SONT EXCLUS 172 : SEGACT MELEME 173 : IPT1 = MELEME 174 : DO 19 IOB=1,MAX(1,LISOUS(/1)) 175 : IF (LISOUS(/1).NE.0) THEN 176 : IPT1=LISOUS(IOB) 177 : SEGACT IPT1 178 : ENDIF 179 : ITY = IPT1.ITYPEL 180 : IF ((idim.eq.2.and.ITY.LT.4).or. 181 : & (idim.eq.3.and.ity.lt.14)) THEN 181 : CALL ERREUR(16) 182 : SEGDES IPT1*NOMOD 183 : RETURN 184 : ENDIF 185 : SEGDES IPT1*NOMOD 186 : 19 CONTINUE 187 : 188 : 189 : C**** TRAVAIL ******************************************************** 190 : 191 : SEGACT MCOORD 192 : NBPTI = XCOOR(/1)/(IDIM+1) 193 : C on veut savoir si il y a des point vus deux fois (interfaces sous zones) 194 : SEGINI IEXCLU 195 : IEXX = IEXCLU 196 : 197 : IPT2 = IPP2 198 : SEGACT IPT2 199 : C ON VA CHANGER LE MAILLAGE LINEAIRE EN ELEMENTS POI1 200 : C SI IL NE N L EST PAS DEJA 201 : IF((IPT2.ITYPEL).NE.1) THEN 201 : CALL CHANGE(IPP2,1) 202 : ENDIF 203 : SEGDES IPT2 204 : 205 : C-------------------------------------- 206 : *---- cas de l'accrochage glissant ---- 207 : IF ( igliss.eq.1) THEN 208 : 209 : * on va attacher à chaque point un vecteur tangent au cable 210 : * si on est arrive avec un modele, on prend le maillage du modele 211 : * en reference 212 : segini itange 213 : idim1=idim+1 214 : segact ipt2 215 : meleme = ipt2 216 : c boucle sur les maillages elementaires 217 : do 3657 i=1,max(1,ipt2.lisous(/1)) 218 : if(ipt2.lisous(/1).ne.0) then 219 : meleme= ipt2.lisous(i) 220 : segact meleme 221 : endif 222 : if(itypel.ne.2) then 222 : call erreur (946) 223 : return 224 : endif 225 : c boucle sur les elements 226 : do 3658 j=1,num(/2) 227 : ia=num(1,j) 228 : ib=num(2,j) 229 : if(iextr(ia).gt.2.Or.iextr(ib).gt.2) then 229 : call erreur (947) 230 : return 231 : endif 232 : xab=xcoor((ib-1)*idim1 + 1)- xcoor((ia-1)*idim1 + 1) 233 : yab=xcoor((ib-1)*idim1 + 2)- xcoor((ia-1)*idim1 + 2) 234 : zab=0. 235 : if(idim.eq.3) then 236 : zab = xcoor((ib-1)*idim1 + 3)- xcoor((ia-1)*idim1 + 3) 237 : endif 238 : xlo=sqrt(xab * xab + yab * yab + zab * zab) 239 : if(xlo.eq.0.) then 239 : call erreur(945) 240 : return 241 : endif 242 : xab=xab/xlo 243 : yab=yab/xlo 244 : zab=zab/xlo 245 : if(iextr(ia).eq.0) then 246 : xperp1(1,ia)=xab 247 : xperp1(2,ia)=yab 248 : if(idim.eq.3) xperp1(3,ia)=zab 249 : else 250 : xxx=1.d0 251 : if( iextr(ia).eq.1) xxx = -1.d0 252 : xperp1(1,ia)=xab*xxx+ xperp1(1,ia) 253 : xperp1(2,ia)=yab*xxx+ xperp1(2,ia) 254 : if(idim.eq.3) xperp1(3,ia)=zab*xxx+ xperp1(3,ia) 255 : xlo=sqrt (xperp1(1,ia)**2+xperp1(2,ia)**2) 256 : if(idim.eq.3)xlo=sqrt (xperp1(1,ia)**2+xperp1(2,ia)**2 257 : $ + xperp1(3,ia)**2) 258 : xperp1(1,ia)=xperp1(1,ia)/xlo 259 : xperp1(2,ia)=xperp1(2,ia)/xlo 260 : if(idim.eq.3)xperp1(3,ia)=xperp1(3,ia)/xlo 261 : endif 262 : if(iextr(ib).eq.0) then 263 : xperp1(1,ib)=xab 264 : xperp1(2,ib)=yab 265 : if(idim.eq.3) xperp1(3,ib)=zab 266 : else 267 : xxx=1.d0 268 : if( iextr(ia).eq.2) xxx = -1.d0 269 : xperp1(1,ib)=xab*xxx+ xperp1(1,ib) 270 : xperp1(2,ib)=yab*xxx+ xperp1(2,ib) 271 : if(idim.eq.3) xperp1(3,ib)=zab*xxx+ xperp1(3,ib) 272 : xlo=sqrt (xperp1(1,ib)**2+xperp1(2,ib)**2) 273 : if(idim.eq.3) xlo=sqrt (xperp1(1,ib)**2+xperp1(2,ib)**2 274 : $ + xperp1(3,ib)**2) 275 : xperp1(1,ib)=xperp1(1,ib)/xlo 276 : xperp1(2,ib)=xperp1(2,ib)/xlo 277 : if(idim.eq.3)xperp1(3,ib)=xperp1(3,ib)/xlo 278 : endif 279 : iextr(ia)= iextr(ia)+1 280 : iextr(ib)= iextr(ib)+2 281 : 3658 continue 282 : c fin de boucle sur les elements 283 : if(ipt2.lisous(/1).ne.0) segdes meleme 284 : 3657 continue 285 : c fin de boucle sur les maillages elementaires 286 : 287 : segdes ipt2 288 : do 3659 i=1,nbpti 289 : if(iextr(i).eq.0)go to 3659 290 : if(abs(xperp1(1,i)).gt.2.e-2.OR.abs(xperp1(2,i)).gt.2.e-2) 291 : $ then 292 : xperp2(1,i)=-xperp1(2,i) 293 : xperp2(2,i)=xperp1(1,i) 294 : else 295 : xperp2(1,i)=-xperp1(3,i) 296 : xperp2(3,i)= xperp1(1,i) 297 : endif 298 : if(idim.eq.3) then 299 : xlo=sqrt (xperp2(1,i)**2+xperp2(2,i)**2 300 : $ + xperp2(3,i)**2) 301 : xperp2(1,i)=xperp2(1,i)/xlo 302 : xperp2(2,i)=xperp2(2,i)/xlo 303 : xperp2(3,i)=xperp2(3,i)/xlo 304 : xaa = xperp1(3,i)*xperp2(2,i)-xperp1(2,i)*xperp2(3,i) 305 : xbb = xperp1(1,i)*xperp2(3,i)-xperp1(3,i)*xperp2(1,i) 306 : xcc = xperp1(2,i)*xperp2(1,i)-xperp1(1,i)*xperp2(2,i) 307 : xperp1(1,i)=xaa 308 : xperp1(2,i)=xbb 309 : xperp1(3,i)=xcc 310 : * write(6,*) 'xaa xbb xcc',xaa,xbb,xcc 311 : else 312 : xperp1(1,i)=xperp2(1,i) 313 : xperp1(2,i)=xperp2(2,i) 314 : endif 315 : 3659 continue 316 : ENDIF 317 : *---- fin du traitement specifique accrochage glissant ---- 318 : C---------------------------------------------------------- 319 : 320 : IPT2 = IPP2 321 : SEGACT IPT2 322 : NODES = IPT2.NUM(/2) 323 : SEGDES IPT2 324 : C 325 : MELEME=MELEMI 326 : SEGACT MELEME 327 : IPP0 = MELEME 328 : NBSZ = MAX(1,LISOUS(/1)) 329 : IPP1 = MELEME 330 : SEGINI ICOSOM 331 : iaccro = 0 332 : 333 : C---- BOUCLE SUR LES SOUS ZONES ---- 334 : DO 22 IOB=1,MAX(1,LISOUS(/1)) 335 : IF (LISOUS(/1).NE.0) THEN 336 : IPP1=LISOUS(IOB) 337 : ENDIF 338 : C ICI LE ZONAGE 339 : C============================================ 339 : CALL ZONAG(IPP1,IPP2,ICC,IEXX,0) 340 : C============================================ 341 : C TABLE DES POINTEURS SUR LES TABLEAUX DE CORRESPONDANCE 342 : C WRITE(6,*) 'POINTEUR ZONE ',IOB,ICC 343 : IPTT(IOB) = ICC 344 : ICOUNT = ICC 345 : SEGACT ICOUNT 346 : iaccro = iaccro + MAXPZ(1) 347 : c do 779 i=1,MAXPZ(1) 348 : C write(6,*)' i ient ', i,IEINT(1,I),IEINT(2,I) 349 : c 779 continue 350 : SEGDES ICOUNT 351 : 22 CONTINUE 352 : C---- FIN DE LA BOUCLE SUR LES SOUS ZONES ---- 353 : 354 : C---- TABLEAU des POINT PAS VUS ---- 355 : NBOUB = NODES - IACCRO 356 : C write(6,*) 'nombre de points exterieurs ' ,nboub 357 : IPT2=IPP2 358 : SEGACT IPT2 359 : SEGINI IOUBLI 360 : SEGACT IEXCLU 361 : IUO = 0 362 : DO 21 I=1,NODES 363 : IP = IPT2.NUM(1,I) 364 : IF(IDEJVU(IP).EQ.0) THEN 365 : IUO = IUO + 1 366 : IPTOU(IUO,1)=IP 367 : ENDIF 368 : 21 CONTINUE 369 : 370 : C----------------------------------------------------------- 371 : C---- RATTRAPAGE EVENTUEL DE POINTS SUR LA LIMITE ---- 372 : IF(IUO.NE.0) THEN 373 : 374 : MELEME=IPP0 375 : SEGACT MELEME 376 : IPT1=MELEME 377 : krat = 0 378 : 379 : C---- boucle sur les sous zones ---- 380 : DO 1010 IOB=1,MAX(1,LISOUS(/1)) 381 : IF (LISOUS(/1).NE.0) THEN 382 : IPT1=LISOUS(IOB) 383 : segact ipt1 384 : ENDIF 385 : 386 : iele=ipt1.itypel 387 : kd = kdegre(iele) 388 : nso = nbsom(iele) 389 : iad=nspos(iele)-1 390 : segini inomil 391 : if(kd.eq.3) then 392 : C element quadratiques 393 : do 762 i=1,ipt1.num(/1) 394 : do 763 j=1,nso 395 : iso = ibsom(iad+j) 396 : if(i.eq.iso) goto 762 397 : 763 continue 398 : inomil(**)= i 399 : 762 continue 400 : elseif(kd.eq.2) then 401 : C elements lineaires 402 : do i=1,ipt1.num(/1) 403 : inomil(**)=i 404 : enddo 405 : endif 406 : NBNN = ipt1.num(/1) 407 : SEGINI MTRAV 408 : ICOUNT= IPTT(IOB) 409 : SEGACT ICOUNT*MOD 410 : 3333 FORMAT(I4,I6,4X,8I4) 411 : 3338 format(I8,3E12.5) 412 : C 413 : C...... boucle sur les elements de la sz ....... 414 : DO 1121 iem=1,ipt1.NUM(/2) 414 : call doxe(xcoor,idim,nbnn,ipt1.num,iem,xel) 415 : 416 : C.......boucle sur les points candidats au rattrapage ...... 417 : DO 1120 IPT = 1,IUO 418 : 419 : IP = IPTOU(IPT,1) 420 : if(ip.LT.0) GOTO 1120 421 : do jl = 1,ipt1.num(/1) 422 : if (ip.eq.ipt1.num(jl,iem)) goto 1120 423 : enddo 424 : C 425 : IREFP= (IP-1)*(IDIM+1)+1 426 : XPO(1)= XCOOR(IREFP) 427 : XPO(2)= XCOOR(IREFP+1) 428 : XPO(3)= 0.D0 429 : IF (IDIM.EQ.3) XPO(3)=XCOOR(IREFP+2) 430 : C coordonnes intrinseques et fonctions de forme : qsijs 431 : C write(6,5777 ) ((xel(ia,ib),ia=1,idim ),ib=1,nbnn) 432 : 5777 format(6e12.5) 432 : call qsijs(xel,iele,nbnn,idim,XPO,SHPP,qsi,iret) 433 : IF (IRET.NE.0) goto 1120 433 : C CALL ERREUR(797) 434 : C RETURN 435 : C ENDIF 436 : do i=1,inomil(/1) 437 : ilp= inomil(i) 438 : if( shpp(1,ilp).lt.ccrit) goto 1120 439 : enddo 440 : C ce point exterieur est rattrapable 441 : Krat = krat+1 442 : NNP = MAXPZ(1)+1 443 : MAXPZ(1)= NNP 444 : IEINT(1,NNP) = ip 445 : iptou(ipt,1) =-ip 446 : IEINT(2,NNP) = iem 447 : do i=1,idim 448 : QQQ(i,NNP) =qsi(i) 449 : enddo 450 : C write(6,2375)ip,iem,(XPO(i1),i1=1,idim),(qsi(i2),i2=1,idim) 451 : 452 : 1120 continue 453 : C.......fin de la boucle sur les points candidats au rattrapage ...... 454 : 455 : 2375 format(2i4,6e12.5) 456 : 1121 CONTINUE 457 : C...... fin de la boucle sur les elements enveloppe de la sz ....... 458 : SEGDES ICOUNT 459 : segsup inomil,MTRAV 460 : IF (LISOUS(/1).NE.0) SEGDES IPT1 461 : 462 : 1010 continue 463 : C---- fin de la boucle sur les zones ---- 464 : 465 : C comptage du nombre de points rattrapes 466 : k= 0 467 : do 777 i=1,IUO 468 : IF(IPTOU(I,1).GT.0) GOTO 777 469 : k= k+1 470 : 777 continue 471 : C 472 : C write(6,*) ' nombre de points rattrapes ' ,k 473 : C on donne le nombre de points non projetés 474 : C if (k.ne.iuo) then 475 : C INTERR(1)=iuo-k 475 : C CALL ERREUR(-322) 476 : C endif 477 : iaccro = iaccro + k 478 : 479 : ENDIF 480 : C---- FIN DE LA SEQUENCE DE RATTRAPAGE 481 : C----------------------------------------------------------- 482 : 483 : c MESSAGE : Nombre de points accrochés %i1 sur %i2 proposés 484 : INTERR(1)=iaccro 485 : INTERR(2)=NODES 485 : CALL ERREUR(-319) 486 : C 487 : C 488 : MELEME = IPP0 489 : SEGACT MELEME 490 : IPT1= MELEME 491 : NRIGEL = 0 492 : NBSZ = 0 493 : NBPTOT = 0 494 : if( igliss.eq.1.and .ifroca.eq.0) ninc=ninc-1 495 : 496 : C NOMBRE DE SOUS ZONES TOUCHÉES ET NOMBRE TOTAL DE NOEUDS CONCERNÉS 497 : DO 399 IOB=1,MAX(1,LISOUS(/1)) 498 : ICOUNT= IPTT(IOB) 499 : C write(6,*) ' segment active' ,IPTT(IOB) 500 : SEGACT ICOUNT 501 : IF(MAXPZ(1).EQ.0) GOTO 398 502 : IF(MAXPZ(1).EQ.0) GOTO 398 503 : NBSZ = NBSZ + 1 504 : NBPTOT = NBPTOT + MAXPZ(1) 505 : 398 SEGDES ICOUNT 506 : 399 CONTINUE 507 : 508 : 509 : C dans le cas ifroca ne 0 les multiplicateurs existent deja 510 : if(ifroca.eq.0) then 511 : NBPTS = NBPTI + NBPTOT*NINC 512 : SEGADJ MCOORD 513 : endif 514 : NRIGEL = NBSZ * NINC 515 : C INITIALISATION DE MRIGID 516 : NRIGE = 8 517 : SEGINI MRIGID 518 : IFORIG = IFOMOD 519 : MTYMAT ='RIGIDITE' 520 : ICZ = 0 521 : INRI=0 522 : C________________________________________________________________ 523 : C____ BOUCLE SUR LES SOUS ZONES __________________________ 524 : DO 400 IOB=1,MAX(1,LISOUS(/1)) 525 : 526 : ICOUNT = IPTT(IOB) 527 : SEGACT ICOUNT 528 : IF (MAXPZ(1).EQ.0) GOTO 395 529 : ICZ = ICZ + 1 530 : C WRITE(6,*) 'SOUS ZONE ',IOB ,'NOMBRE DE POINTS ',MAXPZ(1),ICOUNT 531 : IF (LISOUS(/1).NE.0) THEN 532 : IPT1=LISOUS(IOB) 533 : SEGACT IPT1 534 : ENDIF 535 : KEL = IPT1.ITYPEL 536 : NBNN = IPT1.NUM(/1) 537 : C 538 : SEGINI MTRAV 539 : C ce sera ici 540 : NBELEM = MAXPZ(1) 541 : NBMA1 = IPT1.NUM(/1) 542 : NBNN = NBMA1 + 2 543 : NBSOUS = 0 544 : NBREF = 0 545 : C creation des xmatri 546 : NELRIG = NBELEM 547 : NLIGRD = IPT1.NUM(/1)+1+1 548 : if(igliss.eq.1) nligrd = 1 + idim*(IPT1.NUM(/1)+1) 549 : NLIGRP = NLIGRD 550 : do iop=1,ninc 551 : SEGINI XMATRI 552 : imtt(iop)=Xmatri 553 : enddo 554 : C 555 : INRI =NINC *(ICZ-1)+1 556 : C write(6,*) ' ninc' , ninc 557 : c boucle sur les xmatri 558 : DO 201 I= 1,NINC 559 : SEGINI IPT3 560 : IPT3.ITYPEL= 22 561 : if( ifroca.NE.0.AND.i.eq.ninc) irigel(6,INRI+I-1)=2 562 : C write (6,*) ' ipt3 ',ipt3.num(/1),ipt3.num(/2),ipt1.num(/1), 563 : C > ipt1.num(/2) 564 : IMEL(I)= IPT3 565 : IRIGEL (1,INRI+I-1) = IPT3 566 : IRIGEL (4,INRI+I-1) = IMtt(i) 567 : COERIG(INRI+I-1)= 1.D0 568 : 201 CONTINUE 569 : NELRIG = 0 570 : 571 : C________________________________________________________________ 572 : 573 : C.... boucle sur les elements ....... 574 : DO 300 K = 1,NBELEM 575 : C NELRIG = NELRIG + 1 576 : IP = IEINT(1,K) 577 : IEL = IEINT(2,K) 578 : irefp=(ip-1)*(idim+1)+1 579 : xpo(1)=xcoor(irefp) 580 : xpo(2)=xcoor(irefp+1) 581 : xpo(3)= 0.D0 582 : if(idim.eq.3) xpo(3)=xcoor(irefp+2) 583 : qsi(1)=qqq(1,k) 584 : qsi(2)=qqq(2,k) 585 : qsi(3)=qqq(3,k) 586 : 586 : CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,IRET) 587 : 588 : c -boucle sur les xmatri- 589 : DO 59 IK =1,NINC 590 : 591 : ipt3=imel(ik) 592 : C write(6,*) ' IPT1.NUM(/1) ipt1;num(/2)ipt3.num(/1) ipt3.num(/2)' 593 : C write(6,*) IPT1.NUM(/1), ipt1.num(/2),ipt3.num(/1), ipt3.num(/2) 594 : if(ifroca.eq.0) then 595 : NBPTI = NBPTI + 1 596 : IDE =( NBPTI-1)*(IDIM+1 ) 597 : XCOOR(IDE+1) = XPO(1) 598 : XCOOR(IDE+2) = XPO(2) 599 : IF(IDIM.EQ.3) THEN 600 : XCOOR(IDE+3) = XPO(3) 601 : ENDIF 602 : else 603 : ibp=icpr(ip) 604 : nbpti= imul(ik,ibp) 605 : endif 606 : 607 : C REMPLISSAGE DU MELEME IPT3 608 : IPT3.NUM(1,K) = NBPTI 609 : IPT3.NUM(2,K) = IP 610 : DO 301 LI = 1,IPT1.NUM(/1) 611 : IPT3.NUM(2+LI,K) = IPT1.NUM(LI,IEL) 612 : 301 CONTINUE 613 : 614 : C REMPLISSAGE DU XMATRI 615 : xmatri=imtt(ik) 616 : CC(1) = 1.D0 617 : NHD = IPT1.NUM(/1)+1 618 : DO 305 KM=1,IPT1.NUM(/1) 619 : CC(1+KM) = -SHPP(1,KM) 620 : 305 CONTINUE 621 : 622 : c Cas glissant 623 : IF(igliss.eq.1) THEN 624 : xaa = xperp1(1,ip) 625 : xbb = xperp1(2,ip) 626 : if(idim.eq.3) xcc = xperp1(3,ip) 627 : if( ik.eq.2) then 628 : xaa = xperp2(1,ip) 629 : xbb=xperp2(2,ip) 630 : if(idim.eq.3) xcc= xperp2(3,ip) 631 : endif 632 : if(ifroca.ne.0.and.ik.eq.ninc) then 633 : C petit calcul du vercteur tangent 634 : if( idim.eq.2) then 635 : XAA = xperp1(2,ip) 636 : XBB =-xperp1(1,ip) 637 : else 638 : XAA = xperp1(2,ip)*xperp2(3,ip)-xperp1(3,ip)*xperp2(2,ip) 639 : XBB = xperp1(3,ip)*xperp2(1,ip)-xperp1(1,ip)*xperp2(3,ip) 640 : XCC = xperp1(1,ip)*xperp2(2,ip)-xperp1(2,ip)*xperp2(1,ip) 641 : endif 642 : endif 643 : DO 337 km=1,nhd 644 : ccC = -cc(KM) 645 : RE(1,1+(km-1)*IDIM +1,k)= -CCC*XAA 646 : RE(1,1+(km-1)*IDIM +2,k)= -CCC*XBB 647 : RE(1+(km-1)*IDIM +1,1,k)= -CCC*XAA 648 : RE(1+(km-1)*IDIM +2,1,k)= -CCC*XBB 649 : IF(IDIM.EQ.3) THEN 650 : RE(1,1+(km-1)*IDIM +3,k)= -CCC*XCC 651 : RE(1+(km-1)*IDIM +3,1,k)= -CCC*XCC 652 : ENDIF 653 : 337 CONTINUE 654 : 655 : c Cas non glissant 656 : ELSE 657 : DO 306 KM =1,NHD 658 : CCC = CC(KM) 659 : RE(1,KM+1,k)= CCC 660 : RE(1+KM,1,k)= CCC 661 : 306 CONTINUE 662 : ENDIF 663 : 664 : 59 continue 665 : c - fin de boucle sur les xmatri - 666 : 667 : 300 CONTINUE 668 : C.... fin de boucle sur les elements ....... 669 : SEGDES xMATRI 670 : do ik=1,ninc 671 : IPT3 = IMEL(IK) 672 : SEGDES IPT3 673 : enddo 674 : 675 : C identification des duals des ddl en question 676 : c NL = 20 677 : c bp 18.01.2011 : mieux vaut utiliser LNOMDD renseigné dans cchamp+bdata 678 : NL = LNOMDD 679 : DO 326 I = 1,MOTS(/2) 680 : DO 327 KM = 1,NL 681 : IF(NOMDD(KM).EQ.MOTS(I)) THEN 682 : ICOR(I) = KM 683 : GOTO 326 684 : ENDIF 685 : 327 CONTINUE 686 : 326 CONTINUE 687 : 688 : NLIGRD = IPT1.NUM(/1)+1+1 689 : if(igliss.eq.1) nligrd = 1 + idim*(IPT1.NUM(/1)+1) 690 : NLIGRP = NLIGRD 691 : c---- boucle sur les xmatri ---- 692 : DO 328 I=1,NINC 693 : 694 : NLIGRP = NLIGRD 695 : SEGINI DESCR 696 : LISINC(1)= 'LX' 697 : LISDUA(1)= 'FLX' 698 : NOELEP(1)= 1 699 : NOELED(1)= 1 700 : 701 : c - Cas non glissant - 702 : if(igliss.eq.0) then 703 : DO 307 KM1 = 1,IPT1.NUM(/1)+1 704 : LISINC(KM1+1)= MOTS(I) 705 : LISDUA(KM1+1)= NOMDU(ICOR(I)) 706 : NOELEP(KM1+1)= KM1+1 707 : NOELED(KM1+1)= KM1+1 708 : 307 CONTINUE 709 : 710 : c - Cas glissant - 711 : else 712 : nhd=IPT1.NUM(/1)+1 713 : do 347 km1=1, nhd 714 : lisinc(1+(km1-1)*idim +1)=MOTS(1) 715 : lisinc(1+(km1-1)*idim +2)=MOTS(2) 716 : LISDUA(1 +(km1-1)*idim +1)= NOMDU(ICOR(1)) 717 : LISDUA(1 +(km1-1)*idim +2)= NOMDU(ICOR(2)) 718 : NOELEP(1+(km1-1)*idim +1)= KM1+1 719 : NOELEP(1+(km1-1)*idim +2)= KM1+1 720 : NOELED(1+(km1-1)*idim +1)= KM1+1 721 : NOELED(1+(km1-1)*idim +2)= KM1+1 722 : if(idim.eq.3) then 723 : lisinc(1+(km1-1)*idim +3)=MOTS(3) 724 : LISDUA(1+(km1-1)*idim +3)= NOMDU(ICOR(3)) 725 : NOELEP(1+(km1-1)*idim +3)= KM1+1 726 : NOELED(1+(km1-1)*idim +3)= KM1+1 727 : endif 728 : 347 continue 729 : endif 730 : 731 : SEGDES DESCR 732 : IRIGEL(3,INRI+I-1) = DESCR 733 : IRIGEL(5,INRI+I-1) = NIFOUR 734 : C IRIGEL(6,INRI+I-1) = 0 735 : 328 CONTINUE 736 : c---- fin de boucle sur les xmatri ---- 737 : SEGSUP MTRAV 738 : 739 : IF (LISOUS(/1).NE.0) SEGDES IPT1 740 : 395 CONTINUE 741 : SEGDES ICOUNT 742 : C____ FIN DE BOUCLE SUR LES SOUS ZONES __________________________ 743 : 744 : 400 CONTINUE 745 : 746 : 747 : C**** MENAGE ET ECRITURE DE LA RIGIDTE RESULTATS ********************** 748 : SEGSUP ICOSOM 749 : SEGSUP IEXCLU 750 : 751 : SEGDES MRIGID 752 : SEGDES MELEME , IPT2 753 : SEGDES MLMOTS 753 : CALL ECROBJ('RIGIDITE',MRIGID) 754 : 755 : RETURN 756 : END 757 : 758 : 759 : 760 :
© Cast3M 2003 - All rights reserved.
Disclaimer