Télécharger la plaquette Cast3M


fr en

Télécharger scnd3.eso

Retour à la liste

Numérotation des lignes : 

scnd3

  1. C SCND3     SOURCE    PV090527  26/04/30    21:16:25     12529          
  2. C
  3.       subroutine scnd3(mrigid,ri1,ri4)
  4. C----------------------------------------------------------------------
  5. C Eliminer les inconnues dependantes des matrices elementaires
  6. C
  7. C Entrees :
  8. C ---------
  9. C   MRIGID : pointeur MRIGID sur les rigidites elementaires a conserver.
  10. C   ri1    : pointeur MRIGID sur les rigidites elementaires reecrites
  11. C            sous forme de dependance : indique les noeuds a eliminer.
  12. C
  13. C Sortie  :
  14. C ---------
  15. C   ri4    : pointeur MRIGID sur les rigidites elementaires a conserver
  16. C            ou ont ete eliminees les inconnues dependantes.
  17. C
  18. C Remarque : Le(s) ddl(s) dont dependantes les ddls a eliminer ne sont
  19. C            pas eliminables dans la meme passe.
  20. C----------------------------------------------------------------------
  21.       IMPLICIT INTEGER(I-N)
  22.       IMPLICIT REAL*8 (A-H,O-Z)
  23. C
  24. -INC PPARAM
  25. -INC CCOPTIO
  26. -INC CCHAMP
  27. -INC CCREEL
  28. -INC SMCOORD
  29. -INC SMRIGID
  30. -INC SMELEME
  31. C
  32.       integer oooval
  33. C
  34.       segment icpr(nbno)
  35.       segment trav1
  36.         character*4 compd(nbincl),compdd(nbincl)
  37.       endsegment
  38.       segment trav2
  39.         integer irigd(nbnoc,nbincd),nelrigd(nbnoc,nbincd)
  40.       endsegment
  41.       segment trav3
  42.         integer cor1p(nligrp5),cor1d(nligrd5)
  43.       endsegment
  44.       pointeur trav3b.trav3
  45. C
  46.       segment lisde(3,nbd)
  47.       segment ldesc(6,ndes)
  48.       segment lrigi(nbrig)
  49.       segment lelem(3,nelem)
  50.       segment ltrav
  51.         integer iporig(ncomp),ipnouv(ncomp),iligne(nlr,ncomp)
  52.         integer ipinco(ncomp),ipnode(ncomp),innode(ncomp),nbrel(ncomp)
  53.         real*8  xcoeff(nlr,ncomp)
  54.       endsegment
  55.       pointeur ltravp.ltrav,ltravd.ltrav
  56.       segment linco(nnoeud)
  57.       segment xrela(nligrm5)
  58.       segment lincp(nligrp5)
  59.       segment lincd(nligrd5)
  60.       character*(LOCHPO) coelim,coelid
  61.       logical bvu,bnew,bnsym,bok,ldiag
  62.       integer nnhash(4)
  63.       data nnhash/4*0/
  64. C     =================================================================
  65. C     Partitionner les IRIGEL selon LTRK et taille du XMATRI
  66. C     =================================================================
  67.       LTRK = OOOVAL(1,4)
  68.       IF (LTRK.EQ.0) LTRK = OOOVAL(1,1)
  69.       LTRK = MAX(LTRK,2**24)
  70. C     =================================================================
  71. C     Segments trav1 (noms de composantes) et icpr (numero de noeuds)
  72. C     =================================================================
  73.       ldiag=.false.
  74.       xzref=5.d-12
  75.       nde0=8
  76.       nbdeg=0
  77.       nbincd=0
  78.       nbincl=8
  79.       segini trav1
  80.       nbno=nbpts
  81.       segini icpr
  82.       nbnoc=0
  83.       nrigto=0
  84.       nligp1=0
  85. C
  86. C     Rigidite contenant les inconnues a eliminer
  87.       segact ri1
  88.       nri1=ri1.irigel(/2)
  89.       do 10 irig=1,nri1
  90.         descr=ri1.irigel(3,irig)
  91.         segact descr
  92. C
  93. C       Nombre max d'inco en relation avec le ddl elimine
  94.         lpri=lisinc(/2)
  95.         nligp1=max(nligp1,lpri)
  96. C
  97. C       Composante a eliminer (lisdua(/2)=1)
  98.         coelim=lisdua(1)
  99.         do 50 inc=1,nbincd
  100.           if (compd(inc).eq.coelim) goto 40
  101.  50     continue
  102. C
  103. C       Ajout de la composante
  104.         nbincd=nbincd+1
  105.         if(nbincd.gt.nbincl) then
  106.           nbincl=nbincd+8
  107.           segadj trav1
  108.         endif
  109.         compd(nbincd)=coelim
  110. C
  111.  40     continue
  112. C
  113. C       Identifier les noeuds presents
  114.         meleme=ri1.irigel(1,irig)
  115.         segact meleme
  116.         do j=1,num(/2)
  117.           do k=1,num(/1)
  118.             ip=num(k,j)
  119.             if(icpr(ip).eq.0) then
  120.               nbnoc=nbnoc+1
  121.               icpr(ip)=nbnoc
  122.             endif
  123.           enddo
  124.         enddo
  125. C
  126. C       Boucle sur les composantes dont peut dependre l'inc eliminee
  127.         do 45 iligrd=1,lisinc(/2)
  128.           do 55 inc=1,nbincd
  129.             if (compd(inc).eq.lisinc(iligrd)) goto 45
  130.  55       continue
  131. C
  132. C         Ajout de la composante
  133.           nbincd=nbincd+1
  134.           inc=nbincd
  135.           if(nbincd.gt.nbincl) then
  136.             nbincl=nbincd+8
  137.             segadj trav1
  138.           endif
  139.           compd(nbincd)=lisinc(iligrd)
  140. C
  141.  45     continue
  142.  10   continue
  143. *     write (6,*) ' scnd3 - 1 ',(compd(i),i=1,nbincd)
  144. C
  145. C     Rigidite contenant les matrices elementaires a conserver
  146.       segact mrigid
  147.       nbrig=mrigid.irigel(/2)
  148.       nnoeud=0
  149.       nligrp5=0
  150.       nligrd5=0
  151.       do 110 irig=1,nbrig
  152.         descr=irigel(3,irig)
  153.         segact descr
  154.         lpri=lisinc(/2)
  155.         ldua=lisdua(/2)
  156.         nligrp5=max(nligrp5,lpri)
  157.         nligrd5=max(nligrd5,ldua)
  158.         do 140 iprim=1,lpri
  159.           do 150 ico=1,nbincd
  160.             if (compd(ico).eq.lisinc(iprim)) goto 140
  161.  150      continue
  162. C
  163. C         Ajout de la composante
  164.           nbincd=nbincd+1
  165.           ico=nbincd
  166.           if(nbincd.gt.nbincl) then
  167.             nbincl=nbincd+8
  168.             segadj trav1
  169.           endif
  170.           compd(nbincd)=lisinc(iprim)
  171. C
  172.  140    continue
  173. C
  174. C       Identifier les noeuds presents
  175.         meleme=irigel(1,irig)
  176.         segact meleme
  177.         nnoeud=max(nnoeud,num(/1))
  178.         do j=1,num(/2)
  179.           do i=1,num(/1)
  180.             ip=num(i,j)
  181.             if(icpr(ip).eq.0) then
  182.               nbnoc=nbnoc+1
  183.               icpr(ip)=nbnoc
  184.             endif
  185.           enddo
  186.         enddo
  187.  110  continue
  188.       nligrm5=max(nligrp5,nligrd5)
  189.       segini,linco,xrela,lincp,lincd,trav3
  190. *
  191. *     La matrice ou l'on a elimine les inconnues peut etre plus grosse
  192.       nligrp5=(nligrm5+nligp1)*2
  193.       nligrd5=(nligrm5+nligp1)*2
  194.       segini,trav3b
  195. *     write (6,*) ' scnd3 - 2 ',(compd(i),i=1,nbincd)
  196. C
  197. C     Noms des composantes duales : recuperer dans nomdu
  198. C     eventuellement surcharger/completer par l'utilisateur (opti inco)
  199.       do 300 i=1,nbincd
  200.         do 310 j=1,lnomdd
  201.           if (compd(i).eq.nomdd(j)) goto 320
  202.  310    continue
  203. C
  204. C       write (6,*) 'Inco. non definie ',compd(i),' ',compdd(i)
  205.         moterr(1:LOCHPO)=compd(i)
  206.         call erreur(1150)
  207.         return
  208. C
  209.  320    continue
  210.         compdd(i)=nomdu(j)
  211. **      write (6,*) ' correspondance ', compd(i), compdd(i)
  212.  300  continue
  213. *     write (6,*) ' scnd3 - 3 ',compd(/2),(compd(i),i=1,compd(/2))
  214. *     write (6,*) ' scnd3 - 3 ',compdd(/2),(compdd(i),i=1,compdd(/2))
  215. C     =================================================================
  216. C     Segment trav2 : identifier les numeros des noeuds a eliminer
  217. C     =================================================================
  218.       segini trav2
  219.       do 200 irig=1,nri1
  220.         meleme=ri1.irigel(1,irig)
  221.         descr =ri1.irigel(3,irig)
  222.         do 250 ico=1,nbincd
  223.           if (compd(ico).eq.lisdua(1)) goto 260
  224.  250    continue
  225.         call erreur(5)
  226.  260    continue
  227.         do 240 i=1,num(/2)
  228.           ip=num(noeled(1),i)
  229.           irigd(icpr(ip),ico)=irig
  230.           nelrigd(icpr(ip),ico)=i
  231.  240    continue
  232.  200  continue
  233. C
  234. C     Nombre max de relations dans lesquelles intervient un ddl donne
  235.       nlr=0
  236.       do 201 irig=1,nri1
  237.         meleme=ri1.irigel(1,irig)
  238.         descr =ri1.irigel(3,irig)
  239.         lpri=noelep(/1)
  240.         do 241 i=1,num(/2)
  241.           do 242 j=1,lpri
  242.             ip=num(noelep(j),i)
  243.             do 252 ico=1,nbincd
  244.               if (compd(ico).eq.lisinc(j)) goto 262
  245.  252        continue
  246.             call erreur(5)
  247.  262        continue
  248.             if (nelrigd(icpr(ip),ico).gt.0) goto 242
  249.             nelrigd(icpr(ip),ico)=nelrigd(icpr(ip),ico)-1
  250.             nlr=max(nlr,abs(nelrigd(icpr(ip),ico)))
  251.  242      continue
  252.  241    continue
  253.  201  continue
  254. C     =================================================================
  255. C     Eliminer des matrices elementaires les inc qui doivent l'etre
  256. C     =================================================================
  257.       ndes=nde0
  258.       nbd=nde0
  259.       segini,lrigi,ldesc,lisde
  260.       krigel=0
  261.       kdescr=0
  262. C     ----------------------------
  263. C     Boucle sur les sous-matrices
  264.       do 1000 irig=1,nbrig
  265.         ipt4  =mrigid.irigel(1,irig)
  266.         des4  =mrigid.irigel(3,irig)
  267.         xmatr4=mrigid.irigel(4,irig)
  268.         segact,ipt4,des4,xmatr4
  269. C
  270. C       La matrice est-elle symetrique
  271.         bnsym=(xmatr4.symre.ne.0)
  272. C
  273.         nelrig=xmatr4.re(/3)
  274.         nligrp4=des4.lisinc(/2)
  275.         nligrd4=des4.lisdua(/2)
  276. C
  277. C       Correspondance lisinc, lisdua, numero ico dans trav1
  278.         do 410 iprim=1,nligrp4
  279.           cor1p(iprim)=0
  280.           do 420 ico=1,nbincd
  281.             if (compd(ico).eq.des4.lisinc(iprim)) then
  282.               cor1p(iprim)=ico
  283.               goto 410
  284.             endif
  285.  420      continue
  286. **        write(6,*) 'cor1p pas de correspondance'
  287. **        call erreur(5)
  288.  410    continue
  289. C
  290.         do 510 idua=1,nligrd4
  291.           cor1d(idua)=0
  292.           do 520 ico=1,nbincd
  293.             if (compdd(ico).eq.des4.lisdua(idua)) then
  294.               cor1d(idua)=ico
  295.               goto 510
  296.             endif
  297.  520      continue
  298. **        call erreur(5)
  299.  510    continue
  300. C
  301.         nnoeu=ipt4.num(/1)
  302.         nelem=ipt4.num(/2)
  303.         segini,lelem
  304.         lrigi(irig)=lelem
  305.         bvu=.false.
  306. C
  307. C       ===============================================================
  308. C       Boucle sur les elements
  309. C       ===============================================================
  310.         do 1040 iel=1,nelem
  311. C
  312. C         Reinitialiser linco, lincp et lincd
  313.           do inn=1,nnoeu
  314.             linco(inn)=0
  315.           enddo
  316.           do inn=1,nligrp4
  317.             lincp(inn)=0
  318.           enddo
  319.           if (bnsym) then
  320.             do inn=1,nligrd4
  321.               lincd(inn)=0
  322.             enddo
  323.           endif
  324.           ipass=1
  325. C
  326.           icpt=0
  327.           kelimp=0
  328.           ncor=0
  329.           do 1009 iligri=1,nligrp4
  330.             incp=cor1p(iligri)
  331.             ip=ipt4.num(des4.noelep(iligri),iel)
  332.             inn=des4.noelep(iligri)
  333.             linco(inn)=linco(inn)+1
  334. C
  335. C           Ce ddl est-il a eliminer
  336.             if (nelrigd(icpr(ip),incp).le.0) then
  337.               icpt=icpt+1
  338.               lincp(iligri)=icpt
  339.               goto 1009
  340.             endif
  341. C
  342.             linco(inn)=linco(inn)-1
  343.             irig1=irigd(icpr(ip),incp)
  344.             i1   =nelrigd(icpr(ip),incp)
  345.             des1 =ri1.irigel(3,irig1)
  346.             ncor=ncor+des1.lisinc(/2)
  347.             kelimp=kelimp+1
  348.  1009     continue
  349. C
  350. C         Matrice non symetrique : verifier si la duale de l'inconnue
  351. C         eliminee n'est pas presente
  352.           kelimd=kelimp
  353.           if (bnsym) then
  354.             icpt=0
  355.             kelimd=0
  356.             do 1010 iligrd=1,nligrd4
  357.               incd=cor1d(iligrd)
  358.               ip=ipt4.num(des4.noeled(iligrd),iel)
  359.               inn=des4.noeled(iligrd)
  360.               linco(inn)=linco(inn)+1
  361. C
  362. C             L'inconnue primale associee est-il eliminee
  363.               if (nelrigd(icpr(ip),incd).le.0) then
  364.                 icpt=icpt+1
  365.                 lincd(iligrd)=icpt
  366.                 goto 1010
  367.               endif
  368. C
  369.               linco(inn)=linco(inn)-1
  370.               irig1=irigd(icpr(ip),incd)
  371.               i1   =nelrigd(icpr(ip),incd)
  372.               des1 =ri1.irigel(3,irig1)
  373.               ncor=ncor+des1.lisinc(/2)
  374.               kelimd=kelimd+1
  375.  1010       continue
  376.           endif
  377. C
  378. C         -------------------------------------------------------------
  379. C         Pas de modification a apporter, passer a l'element suivant
  380.           if ((kelimp+kelimd).eq.0) then
  381. C
  382.             if (.not.bvu) then
  383. C
  384.               nnhash(1)=nnhash(1)+1
  385.               call deshas(nbincd,trav3,nnoeu,des4,idhash)
  386. C
  387. C             Tentative de reutilisation d'un descripteur
  388.               mdes=0
  389.               do 2030 ides=kdescr,1,-1
  390.                 idhtmp=lisde(2,ides)
  391.                 if (idhash.ne.idhtmp) then
  392.                   nnhash(2)=nnhash(2)+1
  393.                   goto 2030
  394.                 endif
  395. C
  396.                 des3=lisde(1,ides)
  397.                 if (nligrd4.ne.des3.noeled(/1)) goto 2029
  398.                 if (nligrp4.ne.des3.noelep(/1)) goto 2029
  399.                 do id1=1,nligrd4
  400.                   if (des4.noeled(id1).ne.des3.noeled(id1)) goto 2029
  401.                   if (des4.lisdua(id1).ne.des3.lisdua(id1)) goto 2029
  402.                 enddo
  403.                 do id1=1,nligrp4
  404.                   if (des4.noelep(id1).ne.des3.noelep(id1)) goto 2029
  405.                   if (des4.lisinc(id1).ne.des3.lisinc(id1)) goto 2029
  406.                 enddo
  407. C
  408. C               Reutilisation d'un descripteur existant
  409.                 nnhash(4)=nnhash(4)+1
  410.                 mdes=ides
  411.                 goto 2032
  412. C
  413.  2029           continue
  414.                 nnhash(3)=nnhash(3)+1
  415.  2030         continue
  416. C
  417. C             Nouveau descripteur
  418.               kdescr=kdescr+1
  419.               mdes=kdescr
  420.               if (lisde(/2).lt.kdescr) then
  421.                 nbd=kdescr+nde0
  422.                 segadj,lisde
  423.               endif
  424.               lisde(1,mdes)=des4
  425.               lisde(2,mdes)=idhash
  426.               lisde(3,mdes)=LTRK/(nligrd4*nligrp4)
  427. C
  428.  2032         continue
  429.               bvu=.true.
  430.             endif
  431. C
  432. C           Rigidite avec memes descripteur et caracteristiques?
  433.             iptd=lisde(1,mdes)
  434.             do 2031 irie=krigel,1,-1
  435. C
  436.               if (ldesc(1,irie).ne.iptd)                  goto 2031
  437.               if (ldesc(3,irie).ne.mrigid.irigel(2,irig)) goto 2031
  438.               if (ldesc(4,irie).ne.mrigid.irigel(5,irig)) goto 2031
  439.               if (ldesc(5,irie).ne.mrigid.irigel(6,irig)) goto 2031
  440.               if (ldesc(6,irie).ne.mrigid.irigel(7,irig)) goto 2031
  441.               if (ldesc(2,irie).ge.lisde(3,mdes))         goto 2031
  442. C
  443. C             Nouvel element pour la rigidite
  444.               mrig=irie
  445.               goto 2034
  446. C
  447.  2031       continue
  448. C
  449. C           Nouvelle rigidite
  450.             krigel=krigel+1
  451.             mrig=krigel
  452.             if (ldesc(/2).lt.mrig) then
  453.               ndes=mrig+nde0
  454.               segadj,ldesc
  455.             endif
  456.             ldesc(1,mrig)=lisde(1,mdes)
  457.             ldesc(2,mrig)=0
  458.             ldesc(3,mrig)=mrigid.irigel(2,irig)
  459.             ldesc(4,mrig)=mrigid.irigel(5,irig)
  460.             ldesc(5,mrig)=mrigid.irigel(6,irig)
  461.             ldesc(6,mrig)=mrigid.irigel(7,irig)
  462. C
  463.  2034       continue
  464. C
  465.             lelem(1,iel)=ipt4
  466.             lelem(2,iel)=mrig
  467.             lelem(3,iel)=-1*xmatr4
  468.             ldesc(2,mrig)=ldesc(2,mrig)+1
  469. C
  470.             goto 1040
  471.           endif
  472. C
  473. C         -------------------------------------------------------------
  474. C         Modifications de la matrice elementraire
  475. C
  476. C         Les noeuds sont-ils toujours presents
  477.           nbn1=0
  478.           do ind=1,nnoeu
  479.             if (linco(ind).gt.0) then
  480.               nbn1=nbn1+1
  481.               linco(ind)=nbn1
  482.             endif
  483.           enddo
  484. C
  485.           ncomp=nligrp4+ncor
  486.           segini,ltravp
  487.           ltrav=ltravp
  488.           nbcomp=nligrp4
  489.           ncompp=ncomp
  490.           kelim=kelimp
  491.           kno=0
  492. C
  493.  1111     continue
  494. C
  495.           if (ipass.eq.2) then
  496.             ncomp=nligrd4+ncor
  497.             segini,ltravd
  498.             ltrav=ltravd
  499.             nbcomp=nligrd4
  500.             ncompd=ncomp
  501.             kelim=kelimd
  502.           endif
  503. C
  504.           do icc=1,nbcomp
  505.             ltrav.iporig(icc)=icc
  506.             if (ipass.eq.1) then
  507.               ipn=des4.noelep(icc)
  508.               incp=cor1p(icc)
  509.             else
  510.               ipn=des4.noeled(icc)
  511.               incp=cor1d(icc)
  512.             endif
  513.             if (linco(ipn).ne.0) then
  514.               ip=ipt4.num(ipn,iel)
  515.               ltrav.innode(icc)=ip
  516.               ltrav.ipnode(icc)=linco(ipn)
  517.               ltrav.ipinco(icc)=incp
  518. C             Ce ddl est-il a eliminer
  519.               if (nelrigd(icpr(ip),incp).le.0) then
  520.                 if (ipass.eq.1) ltrav.ipnouv(icc)=lincp(icc)
  521.                 if (ipass.eq.2) ltrav.ipnouv(icc)=lincd(icc)
  522.               endif
  523.             endif
  524.           enddo
  525. C
  526.           kligrp4=0
  527.           kcomp=nbcomp
  528.           do 1011 iligri=1,nbcomp
  529. C
  530.             if (ipass.eq.1) then
  531.               ipn=des4.noelep(iligri)
  532.               incp=cor1p(iligri)
  533.             else
  534.               ipn=des4.noeled(iligri)
  535.               incp=cor1d(iligri)
  536.             endif
  537.             ip=ipt4.num(ipn,iel)
  538. C
  539. C           Ce ddl est-il a eliminer
  540.             if (nelrigd(icpr(ip),incp).le.0) goto 1011
  541. C
  542. C           Le(s) ddl(s) lie(s) au ddl elimine exite(nt)-il(s)
  543.             irig1=irigd(icpr(ip),incp)
  544.             i1    =nelrigd(icpr(ip),incp)
  545.             ipt1  =ri1.irigel(1,irig1)
  546.             des1  =ri1.irigel(3,irig1)
  547.             xmatr1=ri1.irigel(4,irig1)
  548.             segact,xmatr1
  549. C
  550.             xvalm=0.d0
  551.             do icoel=1,des1.lisinc(/2)
  552.               xvalm=max(xvalm,abs(xmatr1.re(1,icoel,i1)))
  553.             enddo
  554.             xvalm=xvalm*xzref
  555. C
  556. C           -----------------------------------------------------------
  557.             do 1012 icoel=1,des1.lisinc(/2)
  558. C
  559. C             Coefficient de la relation non significatif : iterer
  560.               if (abs(xmatr1.re(1,icoel,i1)).le.xvalm) goto 1012
  561. C
  562.               ipl=ipt1.num(des1.noelep(icoel),i1)
  563. C
  564.               do 1420 ico=1,nbincd
  565.                 if (compd(ico).eq.des1.lisinc(icoel)) then
  566.                   iprim=ico
  567.                   goto 1410
  568.                 endif
  569.  1420         continue
  570.  1410         continue
  571. C
  572. C             Noeuds deja presents dans l'element
  573.               ino=0
  574.               ipo=0
  575.               do 1013 idq=1,ncomp
  576.                 if (ltrav.innode(idq).ne.ipl) goto 1013
  577.                 ino=ltrav.ipnode(idq)
  578.                 if (ltrav.ipinco(idq).ne.iprim) goto 1013
  579.                 ipo=idq
  580.                 goto 1014
  581.  1013         continue
  582. C
  583. C             Il faut ajouter un noeud au segment descripteur
  584.               if (ino.eq.0) then
  585.                 kno=kno+1
  586.                 ino=nbn1+kno
  587.               endif
  588. C
  589. C             Il faut completer le segment descripteur
  590.               kligrp4=kligrp4+1
  591. C
  592.  1014         continue
  593. C
  594.               if (ipo.eq.0) then
  595. C               Conserver les informations du nouveau ddl
  596.                 kcomp=kcomp+1
  597.                 ltrav.innode(kcomp)=ipl
  598.                 ltrav.ipnode(kcomp)=ino
  599.                 ltrav.ipinco(kcomp)=iprim
  600.                 ltrav.ipnouv(kcomp)=nbcomp-kelim+kligrp4
  601.                 ipo=kcomp
  602.               endif
  603. C
  604.               ltrav.nbrel(ipo)=ltrav.nbrel(ipo)+1
  605.               ltrav.xcoeff(ltrav.nbrel(ipo),ipo)=xmatr1.re(1,icoel,i1)
  606.               ltrav.iligne(ltrav.nbrel(ipo),ipo)=ltrav.iporig(iligri)
  607. C
  608.  1012       continue
  609. C           -----------------------------------------------------------
  610. C
  611.  1011     continue
  612. C
  613. C         Cas d'une matrice vide
  614.           if (ipass.eq.1) then
  615.             nligrp=nligrp4-kelimp+kligrp4
  616.             if (nligrp.eq.0) goto 1040
  617. C
  618. C           Matrice non-symetrique : refaire le travail sur les duales
  619.             if (bnsym) then
  620.               ipass=2
  621.               goto 1111
  622.             endif
  623. C
  624. C           Matrice symetrique : infos duals identiques
  625.             ltravd=ltravp
  626.             ncompd=ncompp
  627.             nligrd=nligrp
  628.           else
  629.             nligrd=nligrd4-kelimd+kligrp4
  630.             if (nligrd.eq.0) goto 1040
  631.           endif
  632. C
  633. C         Matrice de relation vide car tout a ete elimine
  634.           ityp4=28
  635.           if ((ipt4.itypel.eq.49).or.(ipt4.itypel.eq.22)) then
  636.             iid=1
  637.             if (ipt4.itypel.eq.49) iid=2
  638.             if (nligrp.le.iid) then
  639.               nbdeg=nbdeg+1
  640.               goto 1040
  641.             endif
  642.             ityp4=ipt4.itypel
  643.           endif
  644. C
  645. C         Creation du nouveau maillage, descripteur, matrice
  646. C
  647.           nbnn=nbn1+kno
  648.           nbelem=1
  649.           nbref=0
  650.           nbsous=0
  651.           segini,ipt2
  652.           ipt2.itypel=ityp4
  653. C
  654.           segini,des2
  655. C
  656.           nelrig=1
  657.           rigrel=0
  658.           segini,xmatr2
  659.           xmatr2.symre=xmatr4.symre
  660. C
  661. C         Renseigner le maillage et le descripteur
  662.           do 1021 ic=1,ncompp
  663.             inew=ltravp.ipnouv(ic)
  664.             if (inew.ne.0) then
  665.               des2.noelep(inew)=ltravp.ipnode(ic)
  666.               des2.lisinc(inew)=compd(ltravp.ipinco(ic))
  667.               trav3b.cor1p(inew)=ltravp.ipinco(ic)
  668.               ipt2.num(ltravp.ipnode(ic),1)=ltravp.innode(ic)
  669.             endif
  670.  1021     continue
  671.           do 1022 ic=1,ncompd
  672.             inew=ltravd.ipnouv(ic)
  673.             if (inew.ne.0) then
  674.               des2.noeled(inew)=ltravd.ipnode(ic)
  675.               des2.lisdua(inew)=compdd(ltravd.ipinco(ic))
  676.               trav3b.cor1d(inew)=ltravd.ipinco(ic)
  677.               ipt2.num(ltravd.ipnode(ic),1)=ltravd.innode(ic)
  678.             endif
  679.  1022     continue
  680. C
  681. C         Renseigner la matrice
  682.           do 1024 ic=1,ncompp
  683.             iori=ltravp.iporig(ic)
  684.             inew=ltravp.ipnouv(ic)
  685.             if (inew.ne.0) then
  686. C
  687. C             Reinitialiser xrela
  688.               if (iori.ne.0) then
  689.                 do il=1,nligrd4
  690.                   xrela(il)=xmatr4.re(il,iori,iel)
  691.                 enddo
  692.               else
  693.                 do il=1,nligrd4
  694.                   xrela(il)=0.D0
  695.                 enddo
  696.               endif
  697. C
  698.               do inr=1,ltravp.nbrel(ic)
  699.                 xtmp=ltravp.xcoeff(inr,ic)
  700.                 ilig=ltravp.iligne(inr,ic)
  701.                 do il=1,nligrd4
  702.                   xval=xmatr4.re(il,ilig,iel)*xtmp
  703.                   xrela(il)=xrela(il)+xval
  704.                 enddo
  705.               enddo
  706. C
  707.               do ie=1,ncompd
  708.                 knew=ltravd.ipnouv(ie)
  709.                 if (knew.ne.0) then
  710.                   xval=0.D0
  711.                   if (ie.le.nligrd4) xval=xrela(ie)
  712.                   do inr=1,ltravd.nbrel(ie)
  713.                     xtmp=ltravd.xcoeff(inr,ie)
  714.                     ilig=ltravd.iligne(inr,ie)
  715.                     xval=xval+xtmp*xrela(ilig)
  716.                   enddo
  717.                   xmatr2.re(knew,inew,1)=xmatr2.re(knew,inew,1)+xval
  718.                 endif
  719.               enddo
  720. C
  721.             endif
  722.  1024     continue
  723. C
  724. C         Suppression des termes negligeables de la relation
  725.           if (ipt2.itypel.eq.49.or.ipt2.itypel.eq.22) then
  726. C
  727.             renorm=0.d0
  728.             do ir=1,nligrp4
  729.               do jr=1,nligrd4
  730.                 re1=xmatr4.re(ir,jr,1)
  731.                 renorm=renorm+abs(re1)
  732.               enddo
  733.             enddo
  734.             renorm=max(xpetit,renorm*xzref)
  735. C
  736.             iid=2
  737.             if (ipt2.itypel.eq.49) iid=3
  738. C
  739.             idec=0
  740.             do ii=iid,nligrd
  741.               if (abs(xmatr2.re(ii,1,1)).lt.renorm.and.
  742.      >            abs(xmatr2.re(1,ii,1)).lt.renorm) then
  743.                 idec=idec+1
  744.               elseif(idec.ne.0) then
  745.                 xmatr2.re(ii-idec,1,1)=xmatr2.re(ii,1,1)
  746.                 if(iid.eq.3) xmatr2.re(ii-idec,2,1)=xmatr2.re(ii,2,1)
  747.                 xmatr2.re(1,ii-idec,1)=xmatr2.re(1,ii,1)
  748.                 if(iid.eq.3) xmatr2.re(2,ii-idec,1)=xmatr2.re(2,ii,1)
  749.                 des2.lisinc(ii-idec)=des2.lisinc(ii)
  750.                 des2.lisdua(ii-idec)=des2.lisdua(ii)
  751.                 des2.noelep(ii-idec)=des2.noelep(ii)
  752.                 des2.noeled(ii-idec)=des2.noeled(ii)
  753.                 trav3b.cor1p(ii-idec)=trav3b.cor1p(ii)
  754.                 trav3b.cor1d(ii-idec)=trav3b.cor1d(ii)
  755.               endif
  756.             enddo
  757. C
  758.             nligrd=nligrd-idec
  759.             if(nligrd.lt.iid) then
  760. **        write (6,*) 'scnd3 matrice supprimee LX',nligrd,itypel,meleme
  761.               nbdeg=nbdeg+1
  762.               segsup,ipt2,des2,xmatr2
  763.               goto 1040
  764.             endif
  765. C
  766. C           Ajuster les segments et retirer les noeuds non references
  767.             if (idec.ne.0) then
  768. C
  769.               nligrp=nligrp-idec
  770.               nelrig=1
  771.               rigrel=0
  772.               segadj,des2,xmatr2
  773. C
  774.               do ii=1,nligrp
  775.                 inp=des2.noelep(ii)
  776.                 ipt2.num(inp,1)=-1*abs(ipt2.num(inp,1))
  777.               enddo
  778.               do ii=1,nligrd
  779.                 ind=des2.noeled(ii)
  780.                 ipt2.num(ind,1)=-1*abs(ipt2.num(ind,1))
  781.               enddo
  782. C
  783.               iden=0
  784.               do ii=1,nbnn
  785.                 if (ipt2.num(ii,1).gt.0) then
  786.                   iden=iden+1
  787.                 else
  788.                   if (iden.ne.0) then
  789.                     do no=1,nligrp
  790.                       if (des2.noelep(no).eq.ii) then
  791.                         des2.noelep(no)=des2.noelep(no)-iden
  792.                       endif
  793.                     enddo
  794.                     do no=1,nligrd
  795.                       if (des2.noeled(no).eq.ii) then
  796.                         des2.noeled(no)=des2.noeled(no)-iden
  797.                       endif
  798.                     enddo
  799.                   endif
  800.                   ipt2.num(ii-iden,1)=abs(ipt2.num(ii,1))
  801.                 endif
  802.               enddo
  803.               nbnn=nbnn-iden
  804.               nbelem=1
  805.               nbsous=0
  806.               nbref=0
  807.               segadj,ipt2
  808.             endif
  809. C
  810.           endif
  811. C
  812. C         Matrice symetrique : forcer la symetrie de xmatri
  813.           if(xmatr2.symre.eq.0) then
  814.             ninc=32
  815.             do ib=1,nligrd,ninc
  816.               do jr=1,nligrd
  817.                 do ir=ib,min(jr-1,ib+ninc-1)
  818.                   re1=xmatr2.re(ir,jr,1)
  819.                   re2=xmatr2.re(jr,ir,1)
  820.                   rem=(re1+re2)/2.d0
  821.                   xmatr2.re(ir,jr,1)=rem
  822.                   xmatr2.re(jr,ir,1)=rem
  823.                 enddo
  824.               enddo
  825.             enddo
  826.           endif
  827. C
  828. C         Supprimer les segments de travail
  829.           segsup,ltravp
  830.           if (bnsym) segsup,ltravd
  831. C
  832. C         Tentative de reutilisation d'un descripteur
  833.           nnhash(1)=nnhash(1)+1
  834.           call deshas(nbincd,trav3b,nbnn,des2,idhash)
  835.           jdes=0
  836.           do 1030 ides=kdescr,1,-1
  837.             idhtmp=lisde(2,ides)
  838.             if (idhash.ne.idhtmp) then
  839.               nnhash(2)=nnhash(2)+1
  840.               goto 1030
  841.             endif
  842. C
  843.             des3=lisde(1,ides)
  844.             if (nligrd.ne.des3.noeled(/1)) goto 1029
  845.             if (nligrp.ne.des3.noelep(/1)) goto 1029
  846.             do id1=1,nligrd
  847.               if (des2.noeled(id1).ne.des3.noeled(id1)) goto 1029
  848.               if (des2.lisdua(id1).ne.des3.lisdua(id1)) goto 1029
  849.             enddo
  850.             do id1=1,nligrp
  851.               if (des2.noelep(id1).ne.des3.noelep(id1)) goto 1029
  852.               if (des2.lisinc(id1).ne.des3.lisinc(id1)) goto 1029
  853.             enddo
  854. C
  855. C           Descripteur existant
  856.             nnhash(4)=nnhash(4)+1
  857.             jdes=ides
  858.             segsup,des2
  859.             goto 1032
  860. C
  861.  1029       continue
  862.             nnhash(3)=nnhash(3)+1
  863.  1030     continue
  864. C
  865. C         Nouveau descripteur
  866.           nnhash(3)=nnhash(3)+1
  867.           kdescr=kdescr+1
  868.           jdes=kdescr
  869.           if (lisde(/2).lt.kdescr) then
  870.             nbd=kdescr+nde0
  871.             segadj,lisde
  872.           endif
  873.           lisde(1,jdes)=des2
  874.           lisde(2,jdes)=idhash
  875.           lisde(3,jdes)=LTRK/(nligrd*nligrp)
  876.           goto 1033
  877. C
  878.  1032     continue
  879. C
  880. C         Rigidite avec memes descripteur et caracteristiques?
  881.           iptd=lisde(1,jdes)
  882.           do 1031 irie=krigel,1,-1
  883. C
  884.             if (ldesc(1,irie).ne.iptd)                  goto 1031
  885.             if (ldesc(3,irie).ne.mrigid.irigel(2,irig)) goto 1031
  886.             if (ldesc(4,irie).ne.mrigid.irigel(5,irig)) goto 1031
  887.             if (ldesc(5,irie).ne.mrigid.irigel(6,irig)) goto 1031
  888.             if (ldesc(6,irie).ne.mrigid.irigel(7,irig)) goto 1031
  889.             if (ldesc(2,irie).ge.lisde(3,jdes))         goto 1031
  890. C
  891.             jrig=irie
  892.             goto 1034
  893. C
  894.  1031     continue
  895. C
  896.  1033     continue
  897. C
  898. C         Nouvelle rigidite
  899.           krigel=krigel+1
  900.           jrig=krigel
  901.           if (ldesc(/2).lt.jrig) then
  902.             ndes=jrig+nde0
  903.             segadj,ldesc
  904.           endif
  905.           ldesc(1,jrig)=lisde(1,jdes)
  906.           ldesc(2,jrig)=0
  907.           ldesc(3,jrig)=mrigid.irigel(2,irig)
  908.           ldesc(4,jrig)=mrigid.irigel(5,irig)
  909.           ldesc(5,jrig)=mrigid.irigel(6,irig)
  910.           ldesc(6,jrig)=mrigid.irigel(7,irig)
  911. C
  912.  1034     continue
  913. C
  914.           lelem(1,iel)=ipt2
  915.           lelem(2,iel)=jrig
  916.           lelem(3,iel)=xmatr2
  917.           ldesc(2,jrig)=ldesc(2,jrig)+1
  918. C
  919.  1040   continue
  920. C       ===============================================================
  921. C
  922.         segdes xmatr4
  923.  1000 continue
  924. C     Fin de boucle sur les sous-matrices
  925. C     -----------------------------------
  926.       segsup,linco,xrela,lincp,lincd,trav3,trav3b
  927.       if (ldiag) then
  928.         write(ioimp,*)'Nb calculs hash          =',nnhash(1)
  929.         write(ioimp,*)'Nb comparaisons raccourcies hash=',nnhash(2)
  930.         write(ioimp,*)'Nb comparaisons raccourcies nohash=',nnhash(3)
  931.         write(ioimp,*)'Nb comparaisons completes=',nnhash(4)
  932.         write(ioimp,*)'Nb sous rigidite initiale=',mrigid.irigel(/2)
  933.         write(ioimp,*)'Nb sous rigidite finale  =',krigel
  934.         write(ioimp,*)'Nb de descripteurs       =',kdescr
  935.       endif
  936. C     =================================================================
  937. C     Construire la rigidite ou les inconnues ont ete eliminees
  938. C     =================================================================
  939.       nrigel=krigel
  940.       segini,ri4
  941.       ri4.mtymat='TEMPORAI'
  942.       ri4.iforig=mrigid.iforig
  943.       do ibr=1,nbrig
  944.         lelem=lrigi(ibr)
  945.         ipt4 =mrigid.irigel(1,ibr)
  946. C
  947. C       Renseigner les matrices elementaires
  948.         do 1042 iel=1,ipt4.num(/2)
  949.           ipt6=lelem(1,iel)
  950.           if (ipt6.eq.0) goto 1042
  951.           kdes=lelem(2,iel)
  952.           xmatr6=lelem(3,iel)
  953. C
  954.           itmp=1
  955.           bnew=.true.
  956.           if (xmatr6.lt.0) then
  957.             bnew=.false.
  958.             xmatr6=abs(xmatr6)
  959.             itmp=iel
  960.           endif
  961.           segact xmatr6
  962. C
  963.           des3=ldesc(1,kdes)
  964.           nelt=ldesc(2,kdes)
  965.           ipt5=ri4.irigel(1,kdes)
  966. C
  967. C         Creation de la sous-matrice
  968.           if (ipt5.eq.0) then
  969.             ldesc(2,kdes)=1
  970. C
  971.             nligrp=des3.noelep(/1)
  972.             nligrd=des3.noeled(/1)
  973.             nelrig=nelt
  974. C
  975.             bok=(nelrig.eq.xmatr6.re(/3).and.bnew)
  976.             if (bok) then
  977.               xmatr5=xmatr6
  978.               ipt5=ipt6
  979.               goto 1041
  980.             endif
  981. C
  982. **          write(6,*) 'segini xmatr5',nligrp,nligrd,nelrig
  983.             rigrel=0
  984.             segini,xmatr5
  985.             xmatr5.symre=xmatr6.symre
  986. C
  987.             nbnn=ipt6.num(/1)
  988.             nbelem=nelt
  989.             nbref=0
  990.             nbsous=0
  991.             segini,ipt5
  992.             if (bnew) then
  993.               ipt5.itypel=28
  994.               if ((ipt6.itypel.eq.49).or.(ipt6.itypel.eq.22)) then
  995.                 ipt5.itypel=ipt6.itypel
  996.               endif
  997.             else
  998.               ipt5.itypel=ipt6.itypel
  999.             endif
  1000. C
  1001.  1041       continue
  1002. C
  1003.             ri4.coerig(kdes)=mrigid.coerig(ibr)
  1004.             ri4.irigel(1,kdes)=ipt5
  1005.             ri4.irigel(2,kdes)=mrigid.irigel(2,ibr)
  1006.             ri4.irigel(3,kdes)=des3
  1007.             ri4.irigel(4,kdes)=xmatr5
  1008.             ri4.irigel(5,kdes)=mrigid.irigel(5,ibr)
  1009.             ri4.irigel(6,kdes)=mrigid.irigel(6,ibr)
  1010.             ri4.irigel(7,kdes)=mrigid.irigel(7,ibr)
  1011.             ri4.irigel(8,kdes)=mrigid.irigel(8,ibr)
  1012. C
  1013.             if (bok) goto 1042
  1014.           endif
  1015. C
  1016.           ielmt=ldesc(2,kdes)
  1017.           do inn=1,ipt6.num(/1)
  1018.             ipt5.num(inn,ielmt)=ipt6.num(inn,itmp)
  1019.           enddo
  1020. C
  1021.           xcoer=mrigid.coerig(ibr)/ri4.coerig(kdes)
  1022.           xmatr5=ri4.irigel(4,kdes)
  1023.           segact xmatr5*mod
  1024.           do ilc=1,xmatr5.re(/2)
  1025.             do ili=1,xmatr5.re(/1)
  1026.               xmatr5.re(ili,ilc,ielmt)=xcoer*xmatr6.re(ili,ilc,itmp)
  1027.             enddo
  1028.           enddo
  1029.           segdes xmatr5,xmatr6
  1030. C
  1031.           ldesc(2,kdes)=ldesc(2,kdes)+1
  1032.           if (bnew) segsup,ipt6,xmatr6
  1033.  1042   continue
  1034.         segsup,lelem
  1035.       enddo
  1036.       segsup,icpr,trav1,trav2,lrigi,ldesc
  1037. C
  1038. CCC   WRITE(*,*) 'MRIGID'
  1039. CCC   call prrigi(MRIGID,0)
  1040. CCC   WRITE(*,*) 'RI1'
  1041. CCC   call prrigi(ri1,0)
  1042. CCC   WRITE(*,*) 'RI4'
  1043. CCC   call prrigi(ri4,0)
  1044. C     =================================================================
  1045. C     Desactivation generale
  1046. C     =================================================================
  1047.       segact mrigid
  1048.       do irig=1,nbrig
  1049.         xmatri=irigel(4,irig)
  1050.         meleme=irigel(1,irig)
  1051.         descr=irigel(3,irig)
  1052.         segdes xmatri,descr
  1053.       enddo
  1054.       segact ri4
  1055.       do irig=1,ri4.irigel(/2)
  1056.         xmatri=ri4.irigel(4,irig)
  1057.         meleme=ri4.irigel(1,irig)
  1058.         descr=ri4.irigel(3,irig)
  1059.         segdes xmatri,descr
  1060.       enddo
  1061.       segact ri1
  1062.       do irig=1,nri1
  1063.         xmatri=ri1.irigel(4,irig)
  1064.         meleme=ri1.irigel(1,irig)
  1065.         descr=ri1.irigel(3,irig)
  1066.         segdes xmatri,descr
  1067.       enddo
  1068.       if(iimpi.ne.0.and.nbdeg.ne.0)
  1069.      >  write (6,*) ' nombre de relations degenerees'//
  1070.      >  ' apres elimination ',nbdeg
  1071. *     write(6,*) 'nbincl nligrp5 nnoeud ndes nlr',
  1072. *    >            nbincl,nligrp5,nnoeud,ndes,nlr
  1073.       end
  1074.  
  1075.  
  1076.  
  1077.  
  1078.  
  1079.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales