C SUPRI SOURCE PV090527 26/04/30 21:16:33 12529 SUBROUTINE SUPRI c==================================================================== c sous routine utilisee par l'operateur super option 'rigidite' c c on lit une rigiditee des noeuds maitres(geo ou rigi ou chpoi) c il sort un objet de type superelement c c appelee par super c==================================================================== IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC SMSUPER -INC SMCHPOI -INC SMRIGID -INC SMELEME -INC SMCOORD -INC PPARAM -INC CCOPTIO -INC CCGEOME SEGMENT ICPR(nbpts) SEGMENT JCPR(nbpts) SEGMENT NOINC(NNIN,ITA) SEGMENT NNOEU(NNIN) SEGMENT SNOMIN CHARACTER*(LOCOMP) NOMIN(0) ENDSEGMENT SEGMENT SNOMDU CHARACTER*(LOCOMP) NOMDU(0) ENDSEGMENT SEGMENT ITRANS(LISINC(/2)) SEGMENT INOE(ITA) c segment/xvl/(xva(nligra)*d) c segment ipass(nligra) SEGMENT ISIM(ISIMU) character*4 mcle(1) data mcle/'NOMU'/ logical symetr c * option NOMUltiplicateur nomu=0 call lirmot(mcle,1,nomu,0) CALL LIROBJ ('RIGIDITE',MRIGID,1,IRETOU) segact mrigid IF(IERR.NE.0) RETURN NEWKEQ=1 c c *** recuperation des noeuds maitres c c_______________________________cas du chpo___________________________________ c c CALL LIROBJ ('CHPOINT ',MCHPOI,0,IRETOU) IF(IRETOU.EQ.0) GO TO 1000 c c on vient de lirobj un chpoint on travaille a partir de lui c c creation des inconnues dans nomin c creation de icpr c creation de noinc(i,j)=1 si inconnues i existe pour le jeme noeud c SEGACT MCHPOI SEGINI SNOMIN,ICPR ITA=0 NOMIN(**)='LX ' NNIN=2 DO 1001 I=1,IPCHP(/1) MSOUPO=IPCHP(I) SEGACT MSOUPO MELEME=IGEOC SEGACT MELEME IF(I.EQ.1) THEN NOMIN(**)=NOCOMP(1) ENDIF DO 1002 J=1,NOCOMP(/2) DO 1003 K=1,NOMIN(/2) IF(NOMIN(K).EQ.NOCOMP(J)) GO TO 1002 1003 CONTINUE NNIN=NNIN+1 NOMIN(**)=NOCOMP(J) 1002 CONTINUE DO 1004 J=1,NUM(/2) ICPR(NUM(1,J))=ITA+J 1004 CONTINUE ITA =ITA + NUM(/2) 1001 CONTINUE SEGINI NOINC NTPMAI=ITA ITA=0 DO 1006 I=1,IPCHP(/1) MSOUPO=IPCHP(I) MELEME=IGEOC DO 1007 J=1,NOCOMP(/2) DO 1008 K=1,NNIN IF(NOMIN(K).EQ.NOCOMP(J)) GO TO 1009 1008 CONTINUE 1009 CONTINUE KK=K DO 1010 K=1,NUM(/2) NOINC(KK,K+ITA)=1 1010 CONTINUE 1007 CONTINUE ITA =ITA+NUM(/2) SEGDES MELEME,MSOUPO 1006 CONTINUE SEGDES MCHPOI SEGACT MRIGID c c ** initialisation de nomdu et verifi que chaque noeud et chaque inco c existe DO 1013 I=1,ICPR(/1) ICPR(I)=-ICPR(I) 1013 CONTINUE SEGINI SNOMDU nomdu(**)='FLX ' DO 1014 I=2,NOMIN(/2) NOMDU(**)=' ' 1014 CONTINUE DO 1015 I=1,IRIGEL(/2) MELEME=IRIGEL(1,I) DESCR=IRIGEL(3,I) SEGACT MELEME DO J=1,NUM(/2) DO K=1,NUM(/1) IP=NUM(K,J) ICPR(IP)=ABS(ICPR(IP)) ENDDO ENDDO SEGDES MELEME SEGACT DESCR DO 1017 K=1,NOMIN(/2) DO 1018 J=1,LISINC(/2) IF(NOMIN(K).NE.LISINC(J)) GO TO 1018 NOMDU(K)=LISDUA(J) GO TO 1017 1018 CONTINUE 1017 CONTINUE SEGDES DESCR 1015 CONTINUE DO 1019 I=1,ICPR(/1) IF(ICPR(I).GE.0) GO TO 1019 CALL ERREUR (293) RETURN 1019 CONTINUE DO 1020 I=1,NOMDU(/2) IF(NOMDU(I).NE.' ') GO TO 1020 CALL ERREUR (293) RETURN 1020 CONTINUE NBNN=ITA NBSOUS=0 NBREF=0 NBELEM=1 SEGINI MELEME ITYPEL=1 IMELE=MELEME DO 1021 I=1,ICPR(/1) IF(ICPR(I).EQ.0) GO TO 1021 NUM(ICPR(I),1)=I 1021 CONTINUE SEGDES MELEME GO TO 1011 c c_____________________________cas de rigi___________________________________ c 1000 CONTINUE CALL LIROBJ ('RIGIDITE',MRIG,0,IRETOU) IF(IRETOU.EQ.0) GO TO 1500 ITA=0 NNIN=0 c RI1=MRIG SEGACT RI1 do ir=1,ri1.irigel(/2) enddo c SEGINI SNOMIN,SNOMDU,ICPR nomin(**)='LX ' nomdu(**)='FLX ' c DO 1501 I=1,RI1.IRIGEL(/2) MELEME=RI1.IRIGEL(1,I) SEGACT MELEME DESCR=RI1.IRIGEL(3,I) SEGACT DESCR DO 1502 J=1,LISINC(/2) IF(LISINC(J).EQ.'LX '.AND.J.LE.1) GO TO 1502 DO 1503 K=1,NUM(/2) IP=NUM(NOELEP(J),K) IF(ICPR(IP).NE.0) GO TO 1503 ITA=ITA+1 ICPR(IP)=ITA 1503 CONTINUE IF(NOMIN(/2).EQ.0) THEN NOMIN(**)=LISINC(J) NOMDU(**)=LISDUA(J) ELSE DO 1504 K=1,NOMIN(/2) IF(LISINC(J).EQ.NOMIN(K)) GO TO 1505 1504 CONTINUE NOMIN(**)=LISINC(J) NOMDU(**)=LISDUA(J) 1505 CONTINUE ENDIF 1502 CONTINUE SEGDES MELEME,DESCR 1501 CONTINUE c NNIN=NOMIN(/2) NTPMAI=ITA c SEGINI NOINC c DO 1506 I=1,RI1.IRIGEL(/2) MELEME=RI1.IRIGEL(1,I) DESCR=RI1.IRIGEL(3,I) c SEGACT MELEME,DESCR DO 1507 J=1,LISINC(/2) IF(LISINC(J).EQ.'LX '.AND.J.LE.1) GO TO 1507 IP=NOELEP(J) DO 1509 KK=1,NOMIN(/2) IF(NOMIN(KK).EQ.LISINC(J)) GO TO 1510 1509 CONTINUE 1510 CONTINUE DO 1508 K=1,NUM(/2) IPP=ICPR(NUM(IP,K)) NOINC(KK,IPP)=1 1508 CONTINUE 1507 CONTINUE SEGDES MELEME,DESCR 1506 CONTINUE SEGDES RI1 c NBNN=ITA NBSOUS=0 NBREF=0 NBELEM=1 SEGINI MELEME c ITYPEL=1 IMELE=MELEME c DO 1511 I=1,ICPR(/1) IF(ICPR(I).EQ.0) GO TO 1511 NUM(ICPR(I),1)=I 1511 CONTINUE 1512 CONTINUE c SEGDES MELEME GO TO 1011 c c____________________________cas de geo_______________________________ c 1500 CONTINUE c CALL LIROBJ ('POINT ',MELEME,0,IRETOU) IF(IRETOU.NE.0) CALL CRELEM(MELEME) IF(IRETOU.EQ.0) CALL LIROBJ ('MAILLAGE',MELEME,1,IRETOU) IF(IERR.NE.0) RETURN CALL CHANGE(MELEME,1) SEGINI,IPT1=MELEME SEGDES,IPT1 c c ** on fabrique une numerotation interne uniquement les noeuds c ** maitres. c SEGINI ICPR ITE=0 DO 1 I = 1,NUM(/2) IP= NUM(1,I) IF(ICPR(IP).NE.0) GO TO 1 ITE=ITE+1 ICPR(IP)=ITE 1 CONTINUE NTPMAI=ITE ITA=ITE c c___________on cherche la liste des inconnues pour chaque noeuds___________ c SEGDES MELEME SEGACT MRIGID DESCR=IRIGEL(3,1) SEGACT DESCR SEGINI SNOMIN SEGINI SNOMDU nomin(**)='LX ' nomdu(**)='FLX ' DO 2 I=1,IRIGEL(/2) DESCR=IRIGEL(3,I) SEGACT DESCR DO 4 J=1,LISINC(/2) NO=NOMIN(/2) DO 5 K=1,NO IF(NOMIN(K).EQ.LISINC(J)) GO TO 4 5 CONTINUE NOMIN(**)=LISINC(J) NOMDU(**)=LISDUA(J) 4 CONTINUE SEGDES DESCR 2 CONTINUE NNIN=NOMIN(/2) c c ** on cree le tableau noinc(i,j)=1 si la ieme inconnue existe pour c ** le jeme noeud c ** itrans donne pour chaque descr la correspondance entre lisinc et c *** la liste des inconnues c SEGINI NOINC KN=NOMIN(/2) DO 6 I=1,IRIGEL(/2) DESCR=IRIGEL(3,I) SEGACT DESCR SEGINI ITRANS DO 10 L=1,LISINC(/2) DO 11 M=1,NNIN IF( LISINC(L).NE.NOMIN(M)) GO TO 11 * itrans relie lisinc � nomin ITRANS(L)=M GO TO 10 11 CONTINUE 10 CONTINUE MELEME=IRIGEL(1,I) SEGACT MELEME DO J=1,NUM(/2) DO K = 1,NUM(/1) IP=NUM(K,J) IF(ICPR(IP).NE.0) THEN IT=ICPR(IP) DO 8 L=1,LISINC(/2) IF(NOELEP(L).EQ.K) NOINC(ITRANS(L),IT)=1 8 CONTINUE ENDIF ENDDO ENDDO SEGDES MELEME SEGDES DESCR SEGSUP ITRANS 6 CONTINUE * SEGACT IPT1 NBELEM=1 NBNN=IPT1.NUM(/2) NBREF=0 NBSOUS=0 SEGINI MELEME c do 51 i=1,nbnn c num(i,1)=ipt1.num(1,i) c 51 continue ICOLOR(1)=IDCOUL ITYPEL=28 SEGDES IPT1 DO 52 I=1,ICPR(/1) IF(ICPR(I).EQ.0) GO TO 52 NUM(ICPR(I),1)=I 52 CONTINUE SEGDES MELEME c c ** verification tous les points maitres c DO 16 I=1,ITA DO 17 J=1,NNIN IF(NOINC(J,I).NE.0) GO TO 16 17 CONTINUE WRITE(*,*) 'CAS DE GEO' CALL ERREUR(293) RETURN 16 CONTINUE c 1011 CONTINUE c c ______________________partie commune___________________________ c c * en l'absence du mot cle NOMU, on va rajouter dans les noeuds maitres * les multiplicateurs de lagrange des relations portant sur ceux-ci * ca permet ulterieurement d'utiliser un chargement sur ces * multiplicateurs * * En presence du mot cle NOMU, on ne modifie pas la liste des noeuds * maitres. * c c ** on trie la rigidite initiale pour mettre de cote les bloquages c ** et relations qui concernent uniquement les noeuds maitres. c ** ri4 contiendra la rigidite sans ces bloquages, ri5 contiendra c ** uniquement ces bloquages. on fait deux passages pour pouvoir c ** dimensionner irigel c ** on effectue un traitement particulier pour les ddls de lagrange c ** maitres c segact mrigid segini,ri4=mrigid segini,ri5=mrigid nrig=irigel(/2) symetr=.true. do 800 ir=1,nrig if (irigel(7,ir).ne.0) symetr=.false. ipt3=irigel(1,ir) segact ipt3 if (ipt3.itypel.ne.22) then ri5.irigel(1,ir)=0 goto 800 endif segini,ipt4=ipt3 ri4.irigel(1,ir)=ipt4 segini,ipt5=ipt3 ri5.irigel(1,ir)=ipt5 descr=irigel(3,ir) segact descr xmatri=irigel(4,ir) segact xmatri segini,xmatr4=xmatri ri4.irigel(4,ir)=xmatr4 segini,xmatr5=xmatri ri5.irigel(4,ir)=xmatr5 iel4=0 iel5=0 do 810 iel=1,ipt3.num(/2) iaf=0 ir5=1 do 820 ipt=2,noelep(/1) if (icpr(ipt3.num(noelep(ipt),iel)).ne.0) then do k=1,nomin(/2) if (nomin(k).eq.lisinc(ipt)) goto 821 enddo ir5=0 goto 820 821 continue if (noinc(k,icpr(ipt3.num(noelep(ipt),iel))).eq.1) iaf=1 if (noinc(k,icpr(ipt3.num(noelep(ipt),iel))).eq.0) ir5=0 else ir5=0 endif 820 continue if (iaf.eq.1.and.nomu.eq.0.and.ir5.eq.0) then * il faut rajouter le mult de lagrange dans les noeuds maitres if (icpr(ipt3.num(noelep(1),iel)).eq.0) then ita=ita+1 icpr(ipt3.num(noelep(1),iel))=ita segadj noinc endif * write (6,*) ' mult transforme maitre ', * > icpr(ipt3.num(noelep(1),iel)) noinc(1,icpr(ipt3.num(noelep(1),iel)))=1 endif *** ir5=0 if (ir5.eq.1) then iel5=iel5+1 do ip=1,ipt3.num(/1) ipt5.num(ip,iel5)=ipt3.num(ip,iel) enddo do io=1,re(/2) do iu=1,re(/1) xmatr5.re(iu,io,iel5)=re(iu,io,iel) enddo enddo * imatr5.imattt(iel5)=imattt(iel) else iel4=iel4+1 do ip=1,ipt3.num(/1) ipt4.num(ip,iel4)=ipt3.num(ip,iel) enddo do io=1,re(/2) do iu=1,re(/1) xmatr4.re(iu,io,iel4)=re(iu,io,iel) enddo enddo * imatr4.imattt(iel4)=imattt(iel) endif 810 continue nbnn=ipt3.num(/1) nbsous=0 nbref=0 nbelem=iel5 segadj ipt5 nbelem=iel4 segadj ipt4 nelrig=iel5 segact xmatr5*mod nligrp=xmatr5.re(/2) nligrd=xmatr5.re(/1) rigrel=0 segadj xmatr5 segact xmatr4*mod nligrp=xmatr4.re(/2) nligrd=xmatr4.re(/1) nelrig=iel4 rigrel=0 segadj xmatr4 segdes xmatri,xmatr4,xmatr5 800 continue c c calcul de la rigidite equivalente c * write (6,*) ' ri4 avant dbblx ' * call prrigi(ri4,0) call dbblx(ri4) lagdua=ri4.imlag * write (6,*) ' ri4 avant calkeq ' * call prrigi(ri4,0) * write (6,*) ' ri5 avant calkeq ' * call prrigi(ri5,1) * segdes xmatri * decompte du nombre total de noeuds nbelem=0 segini jcpr segact mrigid do iel=1,irigel(/2) ipt8=irigel(1,iel) segact ipt8 do j=1,ipt8.num(/2) do i=1,ipt8.num(/1) ipt=ipt8.num(i,j) if (jcpr(ipt).eq.0) then nbelem=nbelem+1 jcpr(ipt)=nbelem endif enddo enddo enddo itb=nbelem * ita nombre de noeuds maitre itb nombre de noeuds total * write(6,*) 'ita itb avant calkeq',ita,itb iok=1 if (ita**2. .gt.(50.*itb)**1.3333333333) iok=-1 if(.not.symetr) iok=-2 if(iok.eq.1) > CALL CALKEQ(ri4,NOINC,SNOMIN,ICPR,XMATR1,DES1,ICROUT,iok) c if (ierr.ne.0) return if (iok.le.0) then * pas interessant de calculer le super element * write(6,*) ' pas de super element iok:',iok SEGINI MSUPER MSURAI= MRIGID MRIGTO= MRIGID MCROUT=0 nbnn=1 nbsous=0 nbref=0 segini ipt7 do i=1,nbpts if (jcpr(i).ne.0) ipt7.num(1,jcpr(i))=i enddo segsup jcpr MSUPEL=IPT7 CALL ECROBJ ('SUPERELE',MSUPER) RETURN ENDIF segsup jcpr SEGACT,NOINC,icpr NLIGRA=0 DO I=1,ITA DO J=1,NNIN NLIGRA=NLIGRA+NOINC(J,I) ENDDO ENDDO c c creation de la raideur c c c creation des blocages a ajouter a mrigto c nligrp=1 nligrd=1 nelrig=1 rigrel=0 segini xmatr2 segdes xmatr2 * segini imatr2 * imatr2.imattt(1)=xmatr2 * segdes imatr2 SEGACT RI4 SEGACT DES1 IP=RI4.IRIGEL(/2) NRIGEL=IP+NLIGRA NRIGE=MAX(8,IRIGEL(/1)) SEGINI RI2 ri2.MTYMAT='RIGIDITE' ri2.IFORIG=IFOUR * write (6,*) ' ita ',ita NBNN=ITA NBSOUS=0 NBREF=0 NBELEM=1 SEGINI MELEME c ITYPEL=28 c DO 1611 I=1,ICPR(/1) IF(ICPR(I).EQ.0) GO TO 1611 NUM(ICPR(I),1)=I 1611 CONTINUE * call ecmail(meleme,0) segact meleme imele=meleme NBREF=0 NBSOUS=0 NBNN=1 NBELEM=1 DO 15 I=1,NLIGRA SEGINI IPT1 IPT1.ITYPEL=1 * write (6,*) ' des1.noelep ',des1.noelep(i) IPT1.NUM(1,1)=NUM(DES1.NOELEP(I),1) ri2.irigel(1,i)=ipt1 segdes ipt1 segini des2 des2.lisinc(1)=DES1.LISINC(I) des2.lisdua(1)=DES1.LISdua(I) des2.noelep(1)=1 des2.noeled(1)=1 segdes des2 ri2.irigel(3,i)=des2 ri2.irigel(4,i)=xmatr2 ri2.irigel(5,i)=nifour 15 CONTINUE iplus=0 DO 25 I=1,IP ipt4=ri4.irigel(1,i) segact ipt4 if (ipt4.num(/2).ne.0) then iplus=iplus+1 DO 26 J=1,ri4.irigel(/1) RI2.IRIGEL(J,Iplus+NLIGRA)=ri4.IRIGEL(J,I) 26 CONTINUE RI2.COERIG(Iplus+NLIGRA)=ri4.COERIG(I) endif 25 CONTINUE nrigel=iplus+nligra segadj ri2 c * write (6,*) ' ri4 modifie ' * call prrigi(ri4,1) SEGINI MSUPER MBLOQU=NLIGRA MRIGTO=ri2 SEGDES ri2 NRIGE=8 NELRIG=1 NRIGEL=1 SEGINI MRIGID COERIG(1)=1.D0 * SEGINI IMATRI * IMATTT(1)=XMATR1 MTYMAT='RIGIDITE' IFORIG=IFOUR IRIGEL(1,1)=MELEME IRIGEL(2,1)=0 IRIGEL(3,1)=DES1 IRIGEL(4,1)=xMATR1 IRIGEL(5,1)=NIFOUR IRIGEL(6,1)=0 segact ri5 nrigel5=ri5.irigel(/2) nrigel=1+nrigel5 segadj mrigid iplus=0 do ir=1,ri5.irigel(/2) meleme=ri5.irigel(1,ir) if (meleme.ne.0) then segact meleme if (num(/2).ne.0) then iplus=iplus+1 do id=1,ri5.irigel(/1) irigel(id,1+iplus)=ri5.irigel(id,ir) enddo coerig(1+iplus)=ri5.coerig(ir) endif segdes meleme endif enddo nrigel=1+iplus segadj mrigid c si des inconnues maitres ont ete normalis�e il faut modifier la c la matrice condens�e if (ierr.ne.0) return * segdes xmatri CALL SUPNRM(ICROUT,MRIGID,norr) mdnorr=norr * segdes xmatri segact mrigid*mod c MCROUT=ICROUT msuper.islag=lagdua mrigid.imlag=lagdua * write (6,*) ' msurai lagdua dans supri ',mrigid,lagdua * ici on rajoute meleme=imele SEGDES MELEME * SEGDES xMATRI SEGDES MRIGID SEGSUP ICPR,SNOMIN,SNOMDU,NOINC SEGDES DES1 MSURAI=MRIGID * write (6,*) ' mrigto *************************' * call prrigi(mrigto,0) * write (6,*) ' msurai *************************' * call prrigi(msurai,0) MSUPEL=imele SEGDES MSUPER * CALL ECROBJ ('SUPERELE',MSUPER) * RETURN END