C SCND2 SOURCE PV 11/11/25 21:15:00 7191 * subroutine scnd2(mrigid,ri1,ri4) * * entree * mrigid rigidité sans sous matrices de dependances * ri1 rigidite de dependance * * sortie * ri4 rigidite condensee les noeuds dependants sont * elimines * * on suppose que la premiere inconnue de la raideur de dependance * (celle qu'on elimine) n'apparait que dans cette dependance * Cela va permettre de constituer les tableau : * noeu-inconnue (primale ou duale) ==> numéro raideur de dépendance * au moment de la création des matrices, on va s'efforcer de factoriser * les descripteurs * IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * -INC SMRIGID -INC SMELEME -INC CCOPTIO -INC SMCOORD -INC CCHAMP segment trav1 character*4 compd(nbincd),compdd(nbincd) endsegment segment trav2 integer irigd(nbno,nbincd),nelrigd(nbno,nbincd) endsegment segment trav3 integer cor1p(nligrp),cor1d(nligrd) endsegment segment imatno(lmatno) segment idesno(ldesno) segment igeomno(lgeomno) segment coerno(lcoerno) segment irg7no(lrg7no) segment irg6no(lrg6no) segment irg5no(lrg5no) segment irg2no(lrg2no) segment imatt(nelrig) logical dupli,descar character*4 co1,co1d segment svu logical vu(lvu) endsegment * * write (6,*) ' raideur initiale ' * call prrigi(mrigid,0) * write (6,*) ' raideur dependance ' * call prrigi(ri1,1) lvu=0 lmatno=1000 ldesno=1000 lgeomno=1000 lcoerno=1000 lrg7no=1000 lrg6no=1000 lrg5no=1000 lrg2no=1000 kmatno=0 kdesno=0 kgeomno=0 kcoerno=0 krg7no=0 krg6no=0 krg5no=0 krg2no=0 segini imatno,idesno,igeomno,coerno,irg7no,irg5no,svu segini irg2no,irg6no * remplissage du segment trav1 (liste des noms de composantes) nbincd=0 segini trav1 nbno=xcoor(/1)/(idim+1) segact ri1 do 10 irig=1,ri1.irigel(/2) descr=ri1.irigel(3,irig) segact descr do 40 nligrd=1,lisdua(/2) do 50 inc=1,nbincd if (compd(inc).eq.lisdua(nligrd)) goto 40 50 continue nbincd=nbincd+1 segadj trav1 compd(nbincd)=lisdua(nligrd) 40 continue do 45 nligrd=1,lisinc(/2) do 55 inc=1,nbincd if (compd(inc).eq.lisinc(nligrd)) goto 45 55 continue nbincd=nbincd+1 segadj trav1 compd(nbincd)=lisinc(nligrd) 45 continue 10 continue * write (6,*) ' scnd2 - 1 ',(compd(i),i=1,nbincd) segact mrigid do 110 irig=1,irigel(/2) descr=irigel(3,irig) segact descr do 140 nligrd=1,lisinc(/2) do 150 inc=1,nbincd if (compd(inc).eq.lisinc(nligrd)) goto 140 150 continue nbincd=nbincd+1 segadj trav1 compd(nbincd)=lisinc(nligrd) 140 continue 110 continue * write (6,*) ' scnd2 - 2 ',(compd(i),i=1,nbincd) * au tour de trav2 (inconnue => raideur de dépendance) segini trav2 do 200 irig=1,ri1.irigel(/2) meleme=ri1.irigel(1,irig) descr=ri1.irigel(3,irig) xmatri=ri1.irigel(4,irig) segact meleme,descr,xmatri do 250 incd=1,nbincd if (compd(incd).eq.lisdua(1)) goto 260 250 continue call erreur(5) 260 continue do 240 i=1,num(/2) ip=num(noeled(1),i) irigd(ip,incd)=irig nelrigd(ip,incd)=i 240 continue 200 continue * rechercher les noms d'inconnues duales do 300 i=1,compd(/2) do 310 j=1,lnomdd * write (6,*) ' i j compd nomdd ',i,j,compd(i),nomdd(j) if (compd(i).eq.nomdd(j)) goto 320 310 continue * write (6,*) ' pas de nom dual trouvé ' compdd(i)=compd(i) goto 300 320 continue compdd(i)=nomdu(j) * write (6,*) ' correspondance ', compd(i), compdd(i) 300 continue * write (6,*) ' scnd2 - 3 ',(compdd(i),i=1,compdd(/2)) * * maintenant, on balaye les raideurs. Pour chacune d'entre elle, on * balaye les inconnues pour trouver une dépendance. Si oui on construit * le résultat de l'opération segini,ri4=mrigid ri4.ichole=0 ri4.imgeo1=0 ri4.imgeo2=0 ri4.isupeq=0 * boucle raideur en cours de condensation do 1000 irig=1,ri4.irigel(/2) * write(ioimp,*) 'Balayage raideur ',irig meleme=ri4.irigel(1,irig) segini,ipt4=meleme segdes meleme ri4.irigel(1,irig)=ipt4 des4=ri4.irigel(3,irig) Xmatri=ri4.irigel(4,irig) segini,Xmatr4=Xmatri nelrig=xmatr4.re(/3) segini imatt segdes Xmatri ri4.irigel(4,irig)=Xmatr4 * On regarde le descripteur. S'il est symétrique, on pourra * faire le travail en une passe sinon deux. * Correspondance lisinc, lisdua, numéro incd dans trav1 nligrp=des4.lisinc(/2) nligrd=des4.lisdua(/2) segini trav3 do 410 iligrp=1,nligrp do 420 incd=1,nbincd if (compd(incd).eq.des4.lisinc(iligrp)) then cor1p(iligrp)=incd goto 410 endif 420 continue 410 continue do 510 iligrd=1,nligrd do 520 incd=1,nbincd if (compdd(incd).eq.des4.lisdua(iligrd)) then cor1d(iligrd)=incd goto 510 endif 520 continue 510 continue * write(ioimp,*) 'nligrp=',nligrp * write(ioimp,*) 'nligrd=',nligrd descar=.false. if (nligrp.ne.nligrd) goto 600 do iligr=1,nligrp if (cor1p(iligr).ne.cor1d(iligr)) goto 600 if (des4.noelep(iligr).ne.des4.noeled(iligr)) goto 600 enddo descar=.true. 600 continue * WRITE(IOIMP,194)(des4.NOELEP(J),des4.LISINC(J),cor1p(j),J=1 * $ ,nligrp) * WRITE(IOIMP,194)(des4.NOELED(J),des4.LISDUA(J),cor1d(j),J=1 * $ ,nligrd) * 194 FORMAT( I6,9X,A4,9X,I6) * if (descar) then * write(ioimp,*) 'le descripteur est carre' * else * write(ioimp,*) 'le descripteur n''est pas carre' * endif npass=2 if (descar) npass=1 do 1040 i=1,ipt4.num(/2) dupli=.false. * On fait deux passes : une pour la primale et une pour la duale * sauf si le descripteur est carré do 1045 ipass=1,npass if (ipass.eq.1) then nligr4=des4.lisinc(/2) else nligr4=des4.lisdua(/2) endif do 1010 nligri=1,nligr4 if (ipass.eq.1) then incd=cor1p(nligri) else incd=cor1d(nligri) endif if (incd.eq.0) goto 1010 if (ipass.eq.1) then ip=ipt4.num(des4.noelep(nligri),i) else ip=ipt4.num(des4.noeled(nligri),i) endif * y a t'il une dependance concernant ce noeud if (nelrigd(ip,incd).eq.0) goto 1010 * ici il faut eliminer l'inconnue * on commence par dupliquer la la matrice si ce n'est pas deja fait if (.not.dupli) then ird=irigel(/1) krg7no=krg7no+1 if (krg7no.gt.irg7no(/1)) then lrg7no=irg7no(/1)+1000 segadj irg7no endif if (ird.ge.7) then irg7no(krg7no)=ri4.irigel(7,irig) else irg7no(krg7no)=0 endif krg6no=krg6no+1 if (krg6no.gt.irg6no(/1)) then lrg6no=irg6no(/1)+1000 segadj irg6no endif if (ird.ge.6) then irg6no(krg6no)=ri4.irigel(6,irig) else irg6no(krg6no)=0 endif krg5no=krg5no+1 if (krg5no.gt.irg5no(/1)) then lrg5no=irg5no(/1)+1000 segadj irg5no endif if (ird.ge.5) then irg5no(krg5no)=ri4.irigel(5,irig) else irg5no(krg5no)=0 endif krg2no=krg2no+1 if (krg2no.gt.irg2no(/1)) then lrg2no=irg2no(/1)+1000 segadj irg2no endif irg2no(krg2no)=ri4.irigel(2,irig) segini,descr=des4 kdesno=kdesno+1 if (kdesno.gt.idesno(/1)) then ldesno=idesno(/1)+1000 segadj idesno endif idesno(kdesno)=descr nbnn=ipt4.num(/1) nbelem=1 nbref=0 nbsous=0 segini meleme itypel=28 * si c'est une matrice de relation, il faut garder l'attribut if (ipt4.itypel.eq.22) itypel=22 do i1=1,nbnn num(i1,1)=ipt4.num(i1,i) enddo kgeomno=kgeomno+1 if (kgeomno.gt.igeomno(/1)) then lgeomno=igeomno(/1)+1000 segadj igeomno endif igeomno(kgeomno)=meleme * xmatr4=imatr4.imattt(i) * segini,xmatri=xmatr4 nligrp=xmatr4.re(/2) nligrd=xmatr4.re(/1) nelrig=1 segini xmatri do io=1,nligrp do iu=1,nligrd re(iu,io,1)=xmatr4.re(iu,io,i) enddo enddo kmatno=kmatno+1 if (kmatno.gt.imatno(/1))then lmatno=imatno(/1)+1000 segadj imatno endif imatno(kmatno)=xmatri kcoerno=kcoerno+1 if (kcoerno.gt.coerno(/1))then lcoerno=coerno(/1)+1000 segadj coerno endif coerno(kcoerno)=ri4.coerig(irig) dupli=.true. * imatr4.imattt(i)=0 imatt(i)=1 endif * on elimine dans la matrice dupliquee irig1=irigd(ip,incd) i1=nelrigd(ip,incd) ipt1=ri1.irigel(1,irig1) des1=ri1.irigel(3,irig1) xmatr1=ri1.irigel(4,irig1) segact ipt1,des1,xmatr1 * nligrd1 =des1.lisdua(/2) * nligrd1=1 est l'inconnue à éliminer (de la matrice de dependance) nligrd1=1 nligrp1=des1.lisinc(/2) * write (6,*) ' scnd2 nligrp1 nligrd1 i1 ',nligrp1,nligrd1,i1, * > ip,incp * nlig1 : inconnue primale à remettre dans la matrice (comme primale) do 1050 nlig1=1,nligrp1 ip1=ipt1.num(des1.noelep(nlig1),i1) co1=des1.lisinc(nlig1) if (ipass.eq.2.or.descar) then do incd=1,nbincd if (compd(incd).eq.co1) then co1d=compdd(incd) goto 1055 endif enddo write (6,*) ' scnd2 2 pas de correspondance ' $ ,co1 call erreur(5) return 1055 continue endif * write (6,*) ' ip1 co1 nligrpi ',ip1,co1,nligrpi * le ddl est il deja dans la raideur nbno=0 if (ipass.eq.1) then do 1060 nligr=1,lisinc(/2) if (noelep(nligr).eq.0) goto 1060 if (num(noelep(nligr),1).ne.ip1) goto 1060 * on note le numéro local du noeud pour le reutiliser nbno=noelep(nligr) if (lisinc(nligr).ne.co1) goto 1060 goto 1070 1060 continue else do 2060 nligr=1,lisdua(/2) if (noeled(nligr).eq.0) goto 2060 if (num(noeled(nligr),1).ne.ip1) goto 2060 * on note le numéro local du noeud pour le reutiliser nbno=noeled(nligr) if (lisdua(nligr).ne.co1d) goto 2060 goto 1070 2060 continue endif * on n'a pas trouve le ddl, on le rajoute * write (6,*) ' ajout de ddl' if (descar) then nligr =lisinc(/2)+1 nligrp=lisinc(/2)+1 nligrd=lisdua(/2)+1 else if (ipass.eq.1) then nligr =lisinc(/2)+1 nligrp=lisinc(/2)+1 nligrd=lisdua(/2) else nligr =lisdua(/2)+1 nligrp=lisinc(/2) nligrd=lisdua(/2)+1 endif endif segadj descr segadj xmatri if (nbno.eq.0) then nbsous=0 nbelem=1 nbnn=num(/1)+1 nbref=0 segadj meleme nbno=nbnn endif num(nbno,1)=ip1 if (ipass.eq.1) then noelep(nligrp)=nbno lisinc(nligrp)=co1 endif if (ipass.eq.2.or.descar) then lisdua(nligrd)=co1d noeled(nligrd)=nbno endif 1070 continue * on met a jour la raideur nligrp=lisinc(/2) nligrd=lisdua(/2) * xmatr1=imatr1.imattt(i1) * segact xmatr1 if (ipass.eq.2.or.descar) then do j=1,nligrp re(nligr,j,1)=re(nligr,j,1)+re(nligri,j,1) $ *xmatr1.re(1,nlig1,i1) enddo endif if (ipass.eq.1) then do j=1,nligrd re(j,nligr,1)=re(j,nligr,1)+re(j,nligri,1) $ *xmatr1.re(1,nlig1,i1) enddo if (ipt4.itypel.eq.22) then do ii=3,nligrp if (abs(re(1,ii,1)).gt.1d10) goto 1051 enddo * write (6,*) ' scnd2 relation trop petite ' * write (6,*) 'scnd2 ',(re(1,ii),ii=1,nligrp) * write (6,*) 'scnd2 ',(lisinc(ii),ii=1,nligrp) 1051 continue endif endif * segdes xmatr1 1050 continue * on note les inconnues à supprimer (primale et duale) * write (6,*) ' il faut eliminer ',nligrpi,lisinc(nligrpi) * write (6,*) ' lisinc contient ',(lisinc(iou),iou=1,lisinc(/2)) if (ipass.eq.1) noelep(nligri)=0 if (ipass.eq.2.or.descar) noeled(nligri)=0 1010 continue * * write(ioimp,*) 'Avant suppression des inconnues flaguees' * WRITE(IOIMP,194)(NOELEP(J),LISINC(J),J=1,LISINC(/2)) * WRITE(IOIMP,194)(NOELED(J),LISDUA(J),J=1,LISDUA(/2)) * 194 FORMAT( I6,9X,A4) * WRITE(IOIMP,195) ((RE(L,K,1),K=1,RE(/2)),L=1,RE(/1)) * 195 FORMAT(1X,6E12.5) * on supprime les inconnues flaguées 1045 continue if (dupli) then idecp=0 do 1080 nligr=1,lisinc(/2) if (noelep(nligr).eq.0) then idecp=idecp+1 else lisinc(nligr-idecp)=lisinc(nligr) noelep(nligr-idecp)=noelep(nligr) do nligd=1,lisdua(/2) * re(nligr-idecp,nligd,1)=re(nligr,nligd,1) re(nligd,nligr-idecp,1)=re(nligd,nligr,1) enddo endif 1080 continue idecd=0 do 1090 nligd=1,lisdua(/2) if (noeled(nligd).eq.0) then idecd=idecd+1 else lisdua(nligd-idecd)=lisdua(nligd) noeled(nligd-idecd)=noeled(nligd) do nligr=1,lisinc(/2) * re(nligr,nligd-idecd,1)=re(nligr,nligd,1) re(nligd-idecd,nligr,1)=re(nligd,nligr,1) enddo endif 1090 continue nligrp=lisinc(/2)-idecp nligrd=lisdua(/2)-idecd * on supprime les matrices ne contenant que les multiplicateurs de Lagrange * ou les autres termes non significatifs if (itypel.eq.22) then *** write (6,*) (lisinc(iou),iou=1,lisinc(/2)) if (nligrd.eq.1) then * write (6,*) ' scnd2 matrice supprimee 1 LX' nligrd=0 nligrp=0 endif if (nligrd.eq.2.and.lisinc(2).eq.'LX') then * write (6,*) ' scnd2 matrice supprimee 2 LX' nligrd=0 nligrp=0 endif iid=2 if (nligrd.ge.2.and.lisinc(2).eq.'LX') iid=3 do ii=iid,nligrd if (abs(re(1,ii)).gt.1D-6) goto 1091 enddo * write (6,*) ' scnd2 matrice supprimee LX ',iid-1 nligrd=0 nligrp=0 1091 continue endif segadj xmatri segdes xmatri segadj descr * la matrice est fabriquée, on elimine maintenant les noeuds * non references lvu=num(/1) segadj svu do in=1,vu(/1) vu(in)=.false. enddo do in=1,noelep(/1) vu(noelep(in))=.true. enddo if (.not.descar) then do in=1,noeled(/1) vu(noeled(in))=.true. enddo endif idec=0 do 1095 in=1,num(/1) if (.not.vu(in)) then * noeud a eliminer idec=idec+1 else do no=1,noelep(/1) if (noelep(no).eq.in) then noelep(no)=noelep(no)-idec if (descar) noeled(no)=noeled(no)-idec endif enddo if (.not.descar) then do no=1,noeled(/1) if (noeled(no).eq.in) then noeled(no)=noeled(no)-idec endif enddo endif num(in-idec,1)=num(in,1) endif 1095 continue nbnn=num(/1)-idec nbelem=1 nbsous=0 nbref=0 segadj meleme if (nbnn.eq.0) imatno(kmatno)=0 endif 1040 continue * compression de ipt4 et imatr4 idec=0 * do 2000 i=1,imatr4.imattt(/1) do 2000 i=1,xmatr4.re(/3) * if (imatr4.imattt(i).eq.0) then if(imatt(i).ne.0) then idec=idec+1 else do io=1,xmatr4.re(/2) do iu=1,xmatr4.re(/1) xmatr4.re(iu,io,i-idec)=xmatr4.re(iu,io,i) enddo enddo * imatr4.imattt(i-idec)=imatr4.imattt(i) do in=1,ipt4.num(/1) ipt4.num(in,i-idec)=ipt4.num(in,i) enddo endif 2000 continue * nelrig=imatr4.imattt(/1)-idec nelrig=xmatr4.re(/3)-idec nligrp=xmatr4.re(/2) nligrd=xmatr4.re(/1) segadj Xmatr4 segdes xmatr4 * write (6,*) ' nelrig dans scnd2 ',nelrig if (nelrig.eq.0) ri4.irigel(4,irig)=0 nbnn=ipt4.num(/1) nbelem=ipt4.num(/2)-idec nbref=0 nbsous=0 segadj ipt4 segsup imatt segsup trav3 1000 continue * il ne reste plus qu'a ranger dans ri4 les nouvelles raideurs engendrées nrigeo=ri4.irigel(/2) nrigel=nrigeo+kmatno * nrige=max(ri4.irigel(/1),7) segadj ri4 do 3000 irig=nrigeo+1,nrigel ri4.coerig(irig)=coerno(irig-nrigeo) ri4.irigel(1,irig)=igeomno(irig-nrigeo) ri4.irigel(3,irig)=idesno(irig-nrigeo) ri4.irigel(7,irig)=irg7no(irig-nrigeo) ri4.irigel(5,irig)=irg5no(irig-nrigeo) ri4.irigel(2,irig)=irg2no(irig-nrigeo) ri4.irigel(6,irig)=irg6no(irig-nrigeo) * nelrig=1 * segini imatri * imattt(1)=imatno(irig-nrigeo) * ri4.irigel(4,irig)=imatri ri4.irigel(4,irig)=imatno(irig-nrigeo) * if (imattt(1).eq.0) ri4.irigel(4,irig)=0 3000 continue * compression de irigel idec=0 do 2100 irig=1,ri4.irigel(/2) if (ri4.irigel(4,irig).eq.0) then idec=idec+1 else do in=1,ri4.irigel(/1) ri4.irigel(in,irig-idec)=ri4.irigel(in,irig) enddo ri4.coerig(irig-idec)=ri4.coerig(irig) endif 2100 continue * nrige=ri4.irigel(/1) nrigel=ri4.irigel(/2)-idec * write (6,*) ' nrigel dans scnd2 ',nrigel segadj ri4 * write (6,*) ' raideur resultante ' * call prrigi(ri4) segsup trav1,trav2,imatno,idesno,igeomno,coerno,irg7no $ ,irg5no,irg2no,irg6no,svu * desactivation generale segact mrigid do irig=1,irigel(/2) xmatri=irigel(4,irig) meleme=irigel(1,irig) descr=irigel(3,irig) segdes xmatri,meleme,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,meleme,descr enddo segact ri1 do irig=1,ri1.irigel(/2) xmatri=ri1.irigel(4,irig) meleme=ri1.irigel(1,irig) descr=ri1.irigel(3,irig) segdes xmatri,meleme,descr enddo end