C SCND3 SOURCE PV090527 26/04/30 21:16:25 12529 C subroutine scnd3(mrigid,ri1,ri4) C---------------------------------------------------------------------- C Eliminer les inconnues dependantes des matrices elementaires C C Entrees : C --------- C MRIGID : pointeur MRIGID sur les rigidites elementaires a conserver. C ri1 : pointeur MRIGID sur les rigidites elementaires reecrites C sous forme de dependance : indique les noeuds a eliminer. C C Sortie : C --------- C ri4 : pointeur MRIGID sur les rigidites elementaires a conserver C ou ont ete eliminees les inconnues dependantes. C C Remarque : Le(s) ddl(s) dont dependantes les ddls a eliminer ne sont C pas eliminables dans la meme passe. C---------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C -INC PPARAM -INC CCOPTIO -INC CCHAMP -INC CCREEL -INC SMCOORD -INC SMRIGID -INC SMELEME C integer oooval C segment icpr(nbno) segment trav1 character*4 compd(nbincl),compdd(nbincl) endsegment segment trav2 integer irigd(nbnoc,nbincd),nelrigd(nbnoc,nbincd) endsegment segment trav3 integer cor1p(nligrp5),cor1d(nligrd5) endsegment pointeur trav3b.trav3 C segment lisde(3,nbd) segment ldesc(6,ndes) segment lrigi(nbrig) segment lelem(3,nelem) segment ltrav integer iporig(ncomp),ipnouv(ncomp),iligne(nlr,ncomp) integer ipinco(ncomp),ipnode(ncomp),innode(ncomp),nbrel(ncomp) real*8 xcoeff(nlr,ncomp) endsegment pointeur ltravp.ltrav,ltravd.ltrav segment linco(nnoeud) segment xrela(nligrm5) segment lincp(nligrp5) segment lincd(nligrd5) character*(LOCHPO) coelim,coelid logical bvu,bnew,bnsym,bok,ldiag integer nnhash(4) data nnhash/4*0/ C ================================================================= C Partitionner les IRIGEL selon LTRK et taille du XMATRI C ================================================================= LTRK = OOOVAL(1,4) IF (LTRK.EQ.0) LTRK = OOOVAL(1,1) LTRK = MAX(LTRK,2**24) C ================================================================= C Segments trav1 (noms de composantes) et icpr (numero de noeuds) C ================================================================= ldiag=.false. xzref=5.d-12 nde0=8 nbdeg=0 nbincd=0 nbincl=8 segini trav1 nbno=nbpts segini icpr nbnoc=0 nrigto=0 nligp1=0 C C Rigidite contenant les inconnues a eliminer segact ri1 nri1=ri1.irigel(/2) do 10 irig=1,nri1 descr=ri1.irigel(3,irig) segact descr C C Nombre max d'inco en relation avec le ddl elimine lpri=lisinc(/2) nligp1=max(nligp1,lpri) C C Composante a eliminer (lisdua(/2)=1) coelim=lisdua(1) do 50 inc=1,nbincd if (compd(inc).eq.coelim) goto 40 50 continue C C Ajout de la composante nbincd=nbincd+1 if(nbincd.gt.nbincl) then nbincl=nbincd+8 segadj trav1 endif compd(nbincd)=coelim C 40 continue C C Identifier les noeuds presents meleme=ri1.irigel(1,irig) segact meleme do j=1,num(/2) do k=1,num(/1) ip=num(k,j) if(icpr(ip).eq.0) then nbnoc=nbnoc+1 icpr(ip)=nbnoc endif enddo enddo C C Boucle sur les composantes dont peut dependre l'inc eliminee do 45 iligrd=1,lisinc(/2) do 55 inc=1,nbincd if (compd(inc).eq.lisinc(iligrd)) goto 45 55 continue C C Ajout de la composante nbincd=nbincd+1 inc=nbincd if(nbincd.gt.nbincl) then nbincl=nbincd+8 segadj trav1 endif compd(nbincd)=lisinc(iligrd) C 45 continue 10 continue * write (6,*) ' scnd3 - 1 ',(compd(i),i=1,nbincd) C C Rigidite contenant les matrices elementaires a conserver segact mrigid nbrig=mrigid.irigel(/2) nnoeud=0 nligrp5=0 nligrd5=0 do 110 irig=1,nbrig descr=irigel(3,irig) segact descr lpri=lisinc(/2) ldua=lisdua(/2) nligrp5=max(nligrp5,lpri) nligrd5=max(nligrd5,ldua) do 140 iprim=1,lpri do 150 ico=1,nbincd if (compd(ico).eq.lisinc(iprim)) goto 140 150 continue C C Ajout de la composante nbincd=nbincd+1 ico=nbincd if(nbincd.gt.nbincl) then nbincl=nbincd+8 segadj trav1 endif compd(nbincd)=lisinc(iprim) C 140 continue C C Identifier les noeuds presents meleme=irigel(1,irig) segact meleme nnoeud=max(nnoeud,num(/1)) do j=1,num(/2) do i=1,num(/1) ip=num(i,j) if(icpr(ip).eq.0) then nbnoc=nbnoc+1 icpr(ip)=nbnoc endif enddo enddo 110 continue nligrm5=max(nligrp5,nligrd5) segini,linco,xrela,lincp,lincd,trav3 * * La matrice ou l'on a elimine les inconnues peut etre plus grosse nligrp5=(nligrm5+nligp1)*2 nligrd5=(nligrm5+nligp1)*2 segini,trav3b * write (6,*) ' scnd3 - 2 ',(compd(i),i=1,nbincd) C C Noms des composantes duales : recuperer dans nomdu C eventuellement surcharger/completer par l'utilisateur (opti inco) do 300 i=1,nbincd do 310 j=1,lnomdd if (compd(i).eq.nomdd(j)) goto 320 310 continue C C write (6,*) 'Inco. non definie ',compd(i),' ',compdd(i) moterr(1:LOCHPO)=compd(i) call erreur(1150) return C 320 continue compdd(i)=nomdu(j) ** write (6,*) ' correspondance ', compd(i), compdd(i) 300 continue * write (6,*) ' scnd3 - 3 ',compd(/2),(compd(i),i=1,compd(/2)) * write (6,*) ' scnd3 - 3 ',compdd(/2),(compdd(i),i=1,compdd(/2)) C ================================================================= C Segment trav2 : identifier les numeros des noeuds a eliminer C ================================================================= segini trav2 do 200 irig=1,nri1 meleme=ri1.irigel(1,irig) descr =ri1.irigel(3,irig) do 250 ico=1,nbincd if (compd(ico).eq.lisdua(1)) goto 260 250 continue call erreur(5) 260 continue do 240 i=1,num(/2) ip=num(noeled(1),i) irigd(icpr(ip),ico)=irig nelrigd(icpr(ip),ico)=i 240 continue 200 continue C C Nombre max de relations dans lesquelles intervient un ddl donne nlr=0 do 201 irig=1,nri1 meleme=ri1.irigel(1,irig) descr =ri1.irigel(3,irig) lpri=noelep(/1) do 241 i=1,num(/2) do 242 j=1,lpri ip=num(noelep(j),i) do 252 ico=1,nbincd if (compd(ico).eq.lisinc(j)) goto 262 252 continue call erreur(5) 262 continue if (nelrigd(icpr(ip),ico).gt.0) goto 242 nelrigd(icpr(ip),ico)=nelrigd(icpr(ip),ico)-1 nlr=max(nlr,abs(nelrigd(icpr(ip),ico))) 242 continue 241 continue 201 continue C ================================================================= C Eliminer des matrices elementaires les inc qui doivent l'etre C ================================================================= ndes=nde0 nbd=nde0 segini,lrigi,ldesc,lisde krigel=0 kdescr=0 C ---------------------------- C Boucle sur les sous-matrices do 1000 irig=1,nbrig ipt4 =mrigid.irigel(1,irig) des4 =mrigid.irigel(3,irig) xmatr4=mrigid.irigel(4,irig) segact,ipt4,des4,xmatr4 C C La matrice est-elle symetrique bnsym=(xmatr4.symre.ne.0) C nelrig=xmatr4.re(/3) nligrp4=des4.lisinc(/2) nligrd4=des4.lisdua(/2) C C Correspondance lisinc, lisdua, numero ico dans trav1 do 410 iprim=1,nligrp4 cor1p(iprim)=0 do 420 ico=1,nbincd if (compd(ico).eq.des4.lisinc(iprim)) then cor1p(iprim)=ico goto 410 endif 420 continue ** write(6,*) 'cor1p pas de correspondance' ** call erreur(5) 410 continue C do 510 idua=1,nligrd4 cor1d(idua)=0 do 520 ico=1,nbincd if (compdd(ico).eq.des4.lisdua(idua)) then cor1d(idua)=ico goto 510 endif 520 continue ** call erreur(5) 510 continue C nnoeu=ipt4.num(/1) nelem=ipt4.num(/2) segini,lelem lrigi(irig)=lelem bvu=.false. C C =============================================================== C Boucle sur les elements C =============================================================== do 1040 iel=1,nelem C C Reinitialiser linco, lincp et lincd do inn=1,nnoeu linco(inn)=0 enddo do inn=1,nligrp4 lincp(inn)=0 enddo if (bnsym) then do inn=1,nligrd4 lincd(inn)=0 enddo endif ipass=1 C icpt=0 kelimp=0 ncor=0 do 1009 iligri=1,nligrp4 incp=cor1p(iligri) ip=ipt4.num(des4.noelep(iligri),iel) inn=des4.noelep(iligri) linco(inn)=linco(inn)+1 C C Ce ddl est-il a eliminer if (nelrigd(icpr(ip),incp).le.0) then icpt=icpt+1 lincp(iligri)=icpt goto 1009 endif C linco(inn)=linco(inn)-1 irig1=irigd(icpr(ip),incp) i1 =nelrigd(icpr(ip),incp) des1 =ri1.irigel(3,irig1) ncor=ncor+des1.lisinc(/2) kelimp=kelimp+1 1009 continue C C Matrice non symetrique : verifier si la duale de l'inconnue C eliminee n'est pas presente kelimd=kelimp if (bnsym) then icpt=0 kelimd=0 do 1010 iligrd=1,nligrd4 incd=cor1d(iligrd) ip=ipt4.num(des4.noeled(iligrd),iel) inn=des4.noeled(iligrd) linco(inn)=linco(inn)+1 C C L'inconnue primale associee est-il eliminee if (nelrigd(icpr(ip),incd).le.0) then icpt=icpt+1 lincd(iligrd)=icpt goto 1010 endif C linco(inn)=linco(inn)-1 irig1=irigd(icpr(ip),incd) i1 =nelrigd(icpr(ip),incd) des1 =ri1.irigel(3,irig1) ncor=ncor+des1.lisinc(/2) kelimd=kelimd+1 1010 continue endif C C ------------------------------------------------------------- C Pas de modification a apporter, passer a l'element suivant if ((kelimp+kelimd).eq.0) then C if (.not.bvu) then C nnhash(1)=nnhash(1)+1 call deshas(nbincd,trav3,nnoeu,des4,idhash) C C Tentative de reutilisation d'un descripteur mdes=0 do 2030 ides=kdescr,1,-1 idhtmp=lisde(2,ides) if (idhash.ne.idhtmp) then nnhash(2)=nnhash(2)+1 goto 2030 endif C des3=lisde(1,ides) if (nligrd4.ne.des3.noeled(/1)) goto 2029 if (nligrp4.ne.des3.noelep(/1)) goto 2029 do id1=1,nligrd4 if (des4.noeled(id1).ne.des3.noeled(id1)) goto 2029 if (des4.lisdua(id1).ne.des3.lisdua(id1)) goto 2029 enddo do id1=1,nligrp4 if (des4.noelep(id1).ne.des3.noelep(id1)) goto 2029 if (des4.lisinc(id1).ne.des3.lisinc(id1)) goto 2029 enddo C C Reutilisation d'un descripteur existant nnhash(4)=nnhash(4)+1 mdes=ides goto 2032 C 2029 continue nnhash(3)=nnhash(3)+1 2030 continue C C Nouveau descripteur kdescr=kdescr+1 mdes=kdescr if (lisde(/2).lt.kdescr) then nbd=kdescr+nde0 segadj,lisde endif lisde(1,mdes)=des4 lisde(2,mdes)=idhash lisde(3,mdes)=LTRK/(nligrd4*nligrp4) C 2032 continue bvu=.true. endif C C Rigidite avec memes descripteur et caracteristiques? iptd=lisde(1,mdes) do 2031 irie=krigel,1,-1 C if (ldesc(1,irie).ne.iptd) goto 2031 if (ldesc(3,irie).ne.mrigid.irigel(2,irig)) goto 2031 if (ldesc(4,irie).ne.mrigid.irigel(5,irig)) goto 2031 if (ldesc(5,irie).ne.mrigid.irigel(6,irig)) goto 2031 if (ldesc(6,irie).ne.mrigid.irigel(7,irig)) goto 2031 if (ldesc(2,irie).ge.lisde(3,mdes)) goto 2031 C C Nouvel element pour la rigidite mrig=irie goto 2034 C 2031 continue C C Nouvelle rigidite krigel=krigel+1 mrig=krigel if (ldesc(/2).lt.mrig) then ndes=mrig+nde0 segadj,ldesc endif ldesc(1,mrig)=lisde(1,mdes) ldesc(2,mrig)=0 ldesc(3,mrig)=mrigid.irigel(2,irig) ldesc(4,mrig)=mrigid.irigel(5,irig) ldesc(5,mrig)=mrigid.irigel(6,irig) ldesc(6,mrig)=mrigid.irigel(7,irig) C 2034 continue C lelem(1,iel)=ipt4 lelem(2,iel)=mrig lelem(3,iel)=-1*xmatr4 ldesc(2,mrig)=ldesc(2,mrig)+1 C goto 1040 endif C C ------------------------------------------------------------- C Modifications de la matrice elementraire C C Les noeuds sont-ils toujours presents nbn1=0 do ind=1,nnoeu if (linco(ind).gt.0) then nbn1=nbn1+1 linco(ind)=nbn1 endif enddo C ncomp=nligrp4+ncor segini,ltravp ltrav=ltravp nbcomp=nligrp4 ncompp=ncomp kelim=kelimp kno=0 C 1111 continue C if (ipass.eq.2) then ncomp=nligrd4+ncor segini,ltravd ltrav=ltravd nbcomp=nligrd4 ncompd=ncomp kelim=kelimd endif C do icc=1,nbcomp ltrav.iporig(icc)=icc if (ipass.eq.1) then ipn=des4.noelep(icc) incp=cor1p(icc) else ipn=des4.noeled(icc) incp=cor1d(icc) endif if (linco(ipn).ne.0) then ip=ipt4.num(ipn,iel) ltrav.innode(icc)=ip ltrav.ipnode(icc)=linco(ipn) ltrav.ipinco(icc)=incp C Ce ddl est-il a eliminer if (nelrigd(icpr(ip),incp).le.0) then if (ipass.eq.1) ltrav.ipnouv(icc)=lincp(icc) if (ipass.eq.2) ltrav.ipnouv(icc)=lincd(icc) endif endif enddo C kligrp4=0 kcomp=nbcomp do 1011 iligri=1,nbcomp C if (ipass.eq.1) then ipn=des4.noelep(iligri) incp=cor1p(iligri) else ipn=des4.noeled(iligri) incp=cor1d(iligri) endif ip=ipt4.num(ipn,iel) C C Ce ddl est-il a eliminer if (nelrigd(icpr(ip),incp).le.0) goto 1011 C C Le(s) ddl(s) lie(s) au ddl elimine exite(nt)-il(s) irig1=irigd(icpr(ip),incp) i1 =nelrigd(icpr(ip),incp) ipt1 =ri1.irigel(1,irig1) des1 =ri1.irigel(3,irig1) xmatr1=ri1.irigel(4,irig1) segact,xmatr1 C xvalm=0.d0 do icoel=1,des1.lisinc(/2) xvalm=max(xvalm,abs(xmatr1.re(1,icoel,i1))) enddo xvalm=xvalm*xzref C C ----------------------------------------------------------- do 1012 icoel=1,des1.lisinc(/2) C C Coefficient de la relation non significatif : iterer if (abs(xmatr1.re(1,icoel,i1)).le.xvalm) goto 1012 C ipl=ipt1.num(des1.noelep(icoel),i1) C do 1420 ico=1,nbincd if (compd(ico).eq.des1.lisinc(icoel)) then iprim=ico goto 1410 endif 1420 continue 1410 continue C C Noeuds deja presents dans l'element ino=0 ipo=0 do 1013 idq=1,ncomp if (ltrav.innode(idq).ne.ipl) goto 1013 ino=ltrav.ipnode(idq) if (ltrav.ipinco(idq).ne.iprim) goto 1013 ipo=idq goto 1014 1013 continue C C Il faut ajouter un noeud au segment descripteur if (ino.eq.0) then kno=kno+1 ino=nbn1+kno endif C C Il faut completer le segment descripteur kligrp4=kligrp4+1 C 1014 continue C if (ipo.eq.0) then C Conserver les informations du nouveau ddl kcomp=kcomp+1 ltrav.innode(kcomp)=ipl ltrav.ipnode(kcomp)=ino ltrav.ipinco(kcomp)=iprim ltrav.ipnouv(kcomp)=nbcomp-kelim+kligrp4 ipo=kcomp endif C ltrav.nbrel(ipo)=ltrav.nbrel(ipo)+1 ltrav.xcoeff(ltrav.nbrel(ipo),ipo)=xmatr1.re(1,icoel,i1) ltrav.iligne(ltrav.nbrel(ipo),ipo)=ltrav.iporig(iligri) C 1012 continue C ----------------------------------------------------------- C 1011 continue C C Cas d'une matrice vide if (ipass.eq.1) then nligrp=nligrp4-kelimp+kligrp4 if (nligrp.eq.0) goto 1040 C C Matrice non-symetrique : refaire le travail sur les duales if (bnsym) then ipass=2 goto 1111 endif C C Matrice symetrique : infos duals identiques ltravd=ltravp ncompd=ncompp nligrd=nligrp else nligrd=nligrd4-kelimd+kligrp4 if (nligrd.eq.0) goto 1040 endif C C Matrice de relation vide car tout a ete elimine ityp4=28 if ((ipt4.itypel.eq.49).or.(ipt4.itypel.eq.22)) then iid=1 if (ipt4.itypel.eq.49) iid=2 if (nligrp.le.iid) then nbdeg=nbdeg+1 goto 1040 endif ityp4=ipt4.itypel endif C C Creation du nouveau maillage, descripteur, matrice C nbnn=nbn1+kno nbelem=1 nbref=0 nbsous=0 segini,ipt2 ipt2.itypel=ityp4 C segini,des2 C nelrig=1 rigrel=0 segini,xmatr2 xmatr2.symre=xmatr4.symre C C Renseigner le maillage et le descripteur do 1021 ic=1,ncompp inew=ltravp.ipnouv(ic) if (inew.ne.0) then des2.noelep(inew)=ltravp.ipnode(ic) des2.lisinc(inew)=compd(ltravp.ipinco(ic)) trav3b.cor1p(inew)=ltravp.ipinco(ic) ipt2.num(ltravp.ipnode(ic),1)=ltravp.innode(ic) endif 1021 continue do 1022 ic=1,ncompd inew=ltravd.ipnouv(ic) if (inew.ne.0) then des2.noeled(inew)=ltravd.ipnode(ic) des2.lisdua(inew)=compdd(ltravd.ipinco(ic)) trav3b.cor1d(inew)=ltravd.ipinco(ic) ipt2.num(ltravd.ipnode(ic),1)=ltravd.innode(ic) endif 1022 continue C C Renseigner la matrice do 1024 ic=1,ncompp iori=ltravp.iporig(ic) inew=ltravp.ipnouv(ic) if (inew.ne.0) then C C Reinitialiser xrela if (iori.ne.0) then do il=1,nligrd4 xrela(il)=xmatr4.re(il,iori,iel) enddo else do il=1,nligrd4 xrela(il)=0.D0 enddo endif C do inr=1,ltravp.nbrel(ic) xtmp=ltravp.xcoeff(inr,ic) ilig=ltravp.iligne(inr,ic) do il=1,nligrd4 xval=xmatr4.re(il,ilig,iel)*xtmp xrela(il)=xrela(il)+xval enddo enddo C do ie=1,ncompd knew=ltravd.ipnouv(ie) if (knew.ne.0) then xval=0.D0 if (ie.le.nligrd4) xval=xrela(ie) do inr=1,ltravd.nbrel(ie) xtmp=ltravd.xcoeff(inr,ie) ilig=ltravd.iligne(inr,ie) xval=xval+xtmp*xrela(ilig) enddo xmatr2.re(knew,inew,1)=xmatr2.re(knew,inew,1)+xval endif enddo C endif 1024 continue C C Suppression des termes negligeables de la relation if (ipt2.itypel.eq.49.or.ipt2.itypel.eq.22) then C renorm=0.d0 do ir=1,nligrp4 do jr=1,nligrd4 re1=xmatr4.re(ir,jr,1) renorm=renorm+abs(re1) enddo enddo renorm=max(xpetit,renorm*xzref) C iid=2 if (ipt2.itypel.eq.49) iid=3 C idec=0 do ii=iid,nligrd if (abs(xmatr2.re(ii,1,1)).lt.renorm.and. > abs(xmatr2.re(1,ii,1)).lt.renorm) then idec=idec+1 elseif(idec.ne.0) then xmatr2.re(ii-idec,1,1)=xmatr2.re(ii,1,1) if(iid.eq.3) xmatr2.re(ii-idec,2,1)=xmatr2.re(ii,2,1) xmatr2.re(1,ii-idec,1)=xmatr2.re(1,ii,1) if(iid.eq.3) xmatr2.re(2,ii-idec,1)=xmatr2.re(2,ii,1) des2.lisinc(ii-idec)=des2.lisinc(ii) des2.lisdua(ii-idec)=des2.lisdua(ii) des2.noelep(ii-idec)=des2.noelep(ii) des2.noeled(ii-idec)=des2.noeled(ii) trav3b.cor1p(ii-idec)=trav3b.cor1p(ii) trav3b.cor1d(ii-idec)=trav3b.cor1d(ii) endif enddo C nligrd=nligrd-idec if(nligrd.lt.iid) then ** write (6,*) 'scnd3 matrice supprimee LX',nligrd,itypel,meleme nbdeg=nbdeg+1 segsup,ipt2,des2,xmatr2 goto 1040 endif C C Ajuster les segments et retirer les noeuds non references if (idec.ne.0) then C nligrp=nligrp-idec nelrig=1 rigrel=0 segadj,des2,xmatr2 C do ii=1,nligrp inp=des2.noelep(ii) ipt2.num(inp,1)=-1*abs(ipt2.num(inp,1)) enddo do ii=1,nligrd ind=des2.noeled(ii) ipt2.num(ind,1)=-1*abs(ipt2.num(ind,1)) enddo C iden=0 do ii=1,nbnn if (ipt2.num(ii,1).gt.0) then iden=iden+1 else if (iden.ne.0) then do no=1,nligrp if (des2.noelep(no).eq.ii) then des2.noelep(no)=des2.noelep(no)-iden endif enddo do no=1,nligrd if (des2.noeled(no).eq.ii) then des2.noeled(no)=des2.noeled(no)-iden endif enddo endif ipt2.num(ii-iden,1)=abs(ipt2.num(ii,1)) endif enddo nbnn=nbnn-iden nbelem=1 nbsous=0 nbref=0 segadj,ipt2 endif C endif C C Matrice symetrique : forcer la symetrie de xmatri if(xmatr2.symre.eq.0) then ninc=32 do ib=1,nligrd,ninc do jr=1,nligrd do ir=ib,min(jr-1,ib+ninc-1) re1=xmatr2.re(ir,jr,1) re2=xmatr2.re(jr,ir,1) rem=(re1+re2)/2.d0 xmatr2.re(ir,jr,1)=rem xmatr2.re(jr,ir,1)=rem enddo enddo enddo endif C C Supprimer les segments de travail segsup,ltravp if (bnsym) segsup,ltravd C C Tentative de reutilisation d'un descripteur nnhash(1)=nnhash(1)+1 call deshas(nbincd,trav3b,nbnn,des2,idhash) jdes=0 do 1030 ides=kdescr,1,-1 idhtmp=lisde(2,ides) if (idhash.ne.idhtmp) then nnhash(2)=nnhash(2)+1 goto 1030 endif C des3=lisde(1,ides) if (nligrd.ne.des3.noeled(/1)) goto 1029 if (nligrp.ne.des3.noelep(/1)) goto 1029 do id1=1,nligrd if (des2.noeled(id1).ne.des3.noeled(id1)) goto 1029 if (des2.lisdua(id1).ne.des3.lisdua(id1)) goto 1029 enddo do id1=1,nligrp if (des2.noelep(id1).ne.des3.noelep(id1)) goto 1029 if (des2.lisinc(id1).ne.des3.lisinc(id1)) goto 1029 enddo C C Descripteur existant nnhash(4)=nnhash(4)+1 jdes=ides segsup,des2 goto 1032 C 1029 continue nnhash(3)=nnhash(3)+1 1030 continue C C Nouveau descripteur nnhash(3)=nnhash(3)+1 kdescr=kdescr+1 jdes=kdescr if (lisde(/2).lt.kdescr) then nbd=kdescr+nde0 segadj,lisde endif lisde(1,jdes)=des2 lisde(2,jdes)=idhash lisde(3,jdes)=LTRK/(nligrd*nligrp) goto 1033 C 1032 continue C C Rigidite avec memes descripteur et caracteristiques? iptd=lisde(1,jdes) do 1031 irie=krigel,1,-1 C if (ldesc(1,irie).ne.iptd) goto 1031 if (ldesc(3,irie).ne.mrigid.irigel(2,irig)) goto 1031 if (ldesc(4,irie).ne.mrigid.irigel(5,irig)) goto 1031 if (ldesc(5,irie).ne.mrigid.irigel(6,irig)) goto 1031 if (ldesc(6,irie).ne.mrigid.irigel(7,irig)) goto 1031 if (ldesc(2,irie).ge.lisde(3,jdes)) goto 1031 C jrig=irie goto 1034 C 1031 continue C 1033 continue C C Nouvelle rigidite krigel=krigel+1 jrig=krigel if (ldesc(/2).lt.jrig) then ndes=jrig+nde0 segadj,ldesc endif ldesc(1,jrig)=lisde(1,jdes) ldesc(2,jrig)=0 ldesc(3,jrig)=mrigid.irigel(2,irig) ldesc(4,jrig)=mrigid.irigel(5,irig) ldesc(5,jrig)=mrigid.irigel(6,irig) ldesc(6,jrig)=mrigid.irigel(7,irig) C 1034 continue C lelem(1,iel)=ipt2 lelem(2,iel)=jrig lelem(3,iel)=xmatr2 ldesc(2,jrig)=ldesc(2,jrig)+1 C 1040 continue C =============================================================== C segdes xmatr4 1000 continue C Fin de boucle sur les sous-matrices C ----------------------------------- segsup,linco,xrela,lincp,lincd,trav3,trav3b if (ldiag) then write(ioimp,*)'Nb calculs hash =',nnhash(1) write(ioimp,*)'Nb comparaisons raccourcies hash=',nnhash(2) write(ioimp,*)'Nb comparaisons raccourcies nohash=',nnhash(3) write(ioimp,*)'Nb comparaisons completes=',nnhash(4) write(ioimp,*)'Nb sous rigidite initiale=',mrigid.irigel(/2) write(ioimp,*)'Nb sous rigidite finale =',krigel write(ioimp,*)'Nb de descripteurs =',kdescr endif C ================================================================= C Construire la rigidite ou les inconnues ont ete eliminees C ================================================================= nrigel=krigel segini,ri4 ri4.mtymat='TEMPORAI' ri4.iforig=mrigid.iforig do ibr=1,nbrig lelem=lrigi(ibr) ipt4 =mrigid.irigel(1,ibr) C C Renseigner les matrices elementaires do 1042 iel=1,ipt4.num(/2) ipt6=lelem(1,iel) if (ipt6.eq.0) goto 1042 kdes=lelem(2,iel) xmatr6=lelem(3,iel) C itmp=1 bnew=.true. if (xmatr6.lt.0) then bnew=.false. xmatr6=abs(xmatr6) itmp=iel endif segact xmatr6 C des3=ldesc(1,kdes) nelt=ldesc(2,kdes) ipt5=ri4.irigel(1,kdes) C C Creation de la sous-matrice if (ipt5.eq.0) then ldesc(2,kdes)=1 C nligrp=des3.noelep(/1) nligrd=des3.noeled(/1) nelrig=nelt C bok=(nelrig.eq.xmatr6.re(/3).and.bnew) if (bok) then xmatr5=xmatr6 ipt5=ipt6 goto 1041 endif C ** write(6,*) 'segini xmatr5',nligrp,nligrd,nelrig rigrel=0 segini,xmatr5 xmatr5.symre=xmatr6.symre C nbnn=ipt6.num(/1) nbelem=nelt nbref=0 nbsous=0 segini,ipt5 if (bnew) then ipt5.itypel=28 if ((ipt6.itypel.eq.49).or.(ipt6.itypel.eq.22)) then ipt5.itypel=ipt6.itypel endif else ipt5.itypel=ipt6.itypel endif C 1041 continue C ri4.coerig(kdes)=mrigid.coerig(ibr) ri4.irigel(1,kdes)=ipt5 ri4.irigel(2,kdes)=mrigid.irigel(2,ibr) ri4.irigel(3,kdes)=des3 ri4.irigel(4,kdes)=xmatr5 ri4.irigel(5,kdes)=mrigid.irigel(5,ibr) ri4.irigel(6,kdes)=mrigid.irigel(6,ibr) ri4.irigel(7,kdes)=mrigid.irigel(7,ibr) ri4.irigel(8,kdes)=mrigid.irigel(8,ibr) C if (bok) goto 1042 endif C ielmt=ldesc(2,kdes) do inn=1,ipt6.num(/1) ipt5.num(inn,ielmt)=ipt6.num(inn,itmp) enddo C xcoer=mrigid.coerig(ibr)/ri4.coerig(kdes) xmatr5=ri4.irigel(4,kdes) segact xmatr5*mod do ilc=1,xmatr5.re(/2) do ili=1,xmatr5.re(/1) xmatr5.re(ili,ilc,ielmt)=xcoer*xmatr6.re(ili,ilc,itmp) enddo enddo segdes xmatr5,xmatr6 C ldesc(2,kdes)=ldesc(2,kdes)+1 if (bnew) segsup,ipt6,xmatr6 1042 continue segsup,lelem enddo segsup,icpr,trav1,trav2,lrigi,ldesc C CCC WRITE(*,*) 'MRIGID' CCC call prrigi(MRIGID,0) CCC WRITE(*,*) 'RI1' CCC call prrigi(ri1,0) CCC WRITE(*,*) 'RI4' CCC call prrigi(ri4,0) C ================================================================= C Desactivation generale C ================================================================= segact mrigid do irig=1,nbrig xmatri=irigel(4,irig) meleme=irigel(1,irig) descr=irigel(3,irig) segdes xmatri,descr enddo segact ri4 do irig=1,ri4.irigel(/2) xmatri=ri4.irigel(4,irig) meleme=ri4.irigel(1,irig) descr=ri4.irigel(3,irig) segdes xmatri,descr enddo segact ri1 do irig=1,nri1 xmatri=ri1.irigel(4,irig) meleme=ri1.irigel(1,irig) descr=ri1.irigel(3,irig) segdes xmatri,descr enddo if(iimpi.ne.0.and.nbdeg.ne.0) > write (6,*) ' nombre de relations degenerees'// > ' apres elimination ',nbdeg * write(6,*) 'nbincl nligrp5 nnoeud ndes nlr', * > nbincl,nligrp5,nnoeud,ndes,nlr end