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    MB234859  25/09/30    21:15:14     12374          
  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.           segini,xmatr2
  658.           xmatr2.symre=xmatr4.symre
  659. C
  660. C         Renseigner le maillage et le descripteur
  661.           do 1021 ic=1,ncompp
  662.             inew=ltravp.ipnouv(ic)
  663.             if (inew.ne.0) then
  664.               des2.noelep(inew)=ltravp.ipnode(ic)
  665.               des2.lisinc(inew)=compd(ltravp.ipinco(ic))
  666.               trav3b.cor1p(inew)=ltravp.ipinco(ic)
  667.               ipt2.num(ltravp.ipnode(ic),1)=ltravp.innode(ic)
  668.             endif
  669.  1021     continue
  670.           do 1022 ic=1,ncompd
  671.             inew=ltravd.ipnouv(ic)
  672.             if (inew.ne.0) then
  673.               des2.noeled(inew)=ltravd.ipnode(ic)
  674.               des2.lisdua(inew)=compdd(ltravd.ipinco(ic))
  675.               trav3b.cor1d(inew)=ltravd.ipinco(ic)
  676.               ipt2.num(ltravd.ipnode(ic),1)=ltravd.innode(ic)
  677.             endif
  678.  1022     continue
  679. C
  680. C         Renseigner la matrice
  681.           do 1024 ic=1,ncompp
  682.             iori=ltravp.iporig(ic)
  683.             inew=ltravp.ipnouv(ic)
  684.             if (inew.ne.0) then
  685. C
  686. C             Reinitialiser xrela
  687.               if (iori.ne.0) then
  688.                 do il=1,nligrd4
  689.                   xrela(il)=xmatr4.re(il,iori,iel)
  690.                 enddo
  691.               else
  692.                 do il=1,nligrd4
  693.                   xrela(il)=0.D0
  694.                 enddo
  695.               endif
  696. C
  697.               do inr=1,ltravp.nbrel(ic)
  698.                 xtmp=ltravp.xcoeff(inr,ic)
  699.                 ilig=ltravp.iligne(inr,ic)
  700.                 do il=1,nligrd4
  701.                   xval=xmatr4.re(il,ilig,iel)*xtmp
  702.                   xrela(il)=xrela(il)+xval
  703.                 enddo
  704.               enddo
  705. C
  706.               do ie=1,ncompd
  707.                 knew=ltravd.ipnouv(ie)
  708.                 if (knew.ne.0) then
  709.                   xval=0.D0
  710.                   if (ie.le.nligrd4) xval=xrela(ie)
  711.                   do inr=1,ltravd.nbrel(ie)
  712.                     xtmp=ltravd.xcoeff(inr,ie)
  713.                     ilig=ltravd.iligne(inr,ie)
  714.                     xval=xval+xtmp*xrela(ilig)
  715.                   enddo
  716.                   xmatr2.re(knew,inew,1)=xmatr2.re(knew,inew,1)+xval
  717.                 endif
  718.               enddo
  719. C
  720.             endif
  721.  1024     continue
  722. C
  723. C         Suppression des termes negligeables de la relation
  724.           if (ipt2.itypel.eq.49.or.ipt2.itypel.eq.22) then
  725. C
  726.             renorm=0.d0
  727.             do ir=1,nligrp4
  728.               do jr=1,nligrd4
  729.                 re1=xmatr4.re(ir,jr,1)
  730.                 renorm=renorm+abs(re1)
  731.               enddo
  732.             enddo
  733.             renorm=max(xpetit,renorm*xzref)
  734. C
  735.             iid=2
  736.             if (ipt2.itypel.eq.49) iid=3
  737. C
  738.             idec=0
  739.             do ii=iid,nligrd
  740.               if (abs(xmatr2.re(ii,1,1)).lt.renorm.and.
  741.      >            abs(xmatr2.re(1,ii,1)).lt.renorm) then
  742.                 idec=idec+1
  743.               elseif(idec.ne.0) then
  744.                 xmatr2.re(ii-idec,1,1)=xmatr2.re(ii,1,1)
  745.                 if(iid.eq.3) xmatr2.re(ii-idec,2,1)=xmatr2.re(ii,2,1)
  746.                 xmatr2.re(1,ii-idec,1)=xmatr2.re(1,ii,1)
  747.                 if(iid.eq.3) xmatr2.re(2,ii-idec,1)=xmatr2.re(2,ii,1)
  748.                 des2.lisinc(ii-idec)=des2.lisinc(ii)
  749.                 des2.lisdua(ii-idec)=des2.lisdua(ii)
  750.                 des2.noelep(ii-idec)=des2.noelep(ii)
  751.                 des2.noeled(ii-idec)=des2.noeled(ii)
  752.                 trav3b.cor1p(ii-idec)=trav3b.cor1p(ii)
  753.                 trav3b.cor1d(ii-idec)=trav3b.cor1d(ii)
  754.               endif
  755.             enddo
  756. C
  757.             nligrd=nligrd-idec
  758.             if(nligrd.lt.iid) then
  759. **        write (6,*) 'scnd3 matrice supprimee LX',nligrd,itypel,meleme
  760.               nbdeg=nbdeg+1
  761.               segsup,ipt2,des2,xmatr2
  762.               goto 1040
  763.             endif
  764. C
  765. C           Ajuster les segments et retirer les noeuds non references
  766.             if (idec.ne.0) then
  767. C
  768.               nligrp=nligrp-idec
  769.               nelrig=1
  770.               segadj,des2,xmatr2
  771. C
  772.               do ii=1,nligrp
  773.                 inp=des2.noelep(ii)
  774.                 ipt2.num(inp,1)=-1*abs(ipt2.num(inp,1))
  775.               enddo
  776.               do ii=1,nligrd
  777.                 ind=des2.noeled(ii)
  778.                 ipt2.num(ind,1)=-1*abs(ipt2.num(ind,1))
  779.               enddo
  780. C
  781.               iden=0
  782.               do ii=1,nbnn
  783.                 if (ipt2.num(ii,1).gt.0) then
  784.                   iden=iden+1
  785.                 else
  786.                   if (iden.ne.0) then
  787.                     do no=1,nligrp
  788.                       if (des2.noelep(no).eq.ii) then
  789.                         des2.noelep(no)=des2.noelep(no)-iden
  790.                       endif
  791.                     enddo
  792.                     do no=1,nligrd
  793.                       if (des2.noeled(no).eq.ii) then
  794.                         des2.noeled(no)=des2.noeled(no)-iden
  795.                       endif
  796.                     enddo
  797.                   endif
  798.                   ipt2.num(ii-iden,1)=abs(ipt2.num(ii,1))
  799.                 endif
  800.               enddo
  801.               nbnn=nbnn-iden
  802.               nbelem=1
  803.               nbsous=0
  804.               nbref=0
  805.               segadj,ipt2
  806.             endif
  807. C
  808.           endif
  809. C
  810. C         Matrice symetrique : forcer la symetrie de xmatri
  811.           if(xmatr2.symre.eq.0) then
  812.             ninc=32
  813.             do ib=1,nligrd,ninc
  814.               do jr=1,nligrd
  815.                 do ir=ib,min(jr-1,ib+ninc-1)
  816.                   re1=xmatr2.re(ir,jr,1)
  817.                   re2=xmatr2.re(jr,ir,1)
  818.                   rem=(re1+re2)/2.d0
  819.                   xmatr2.re(ir,jr,1)=rem
  820.                   xmatr2.re(jr,ir,1)=rem
  821.                 enddo
  822.               enddo
  823.             enddo
  824.           endif
  825. C
  826. C         Supprimer les segments de travail
  827.           segsup,ltravp
  828.           if (bnsym) segsup,ltravd
  829. C
  830. C         Tentative de reutilisation d'un descripteur
  831.           nnhash(1)=nnhash(1)+1
  832.           call deshas(nbincd,trav3b,nbnn,des2,idhash)
  833.           jdes=0
  834.           do 1030 ides=kdescr,1,-1
  835.             idhtmp=lisde(2,ides)
  836.             if (idhash.ne.idhtmp) then
  837.               nnhash(2)=nnhash(2)+1
  838.               goto 1030
  839.             endif
  840. C
  841.             des3=lisde(1,ides)
  842.             if (nligrd.ne.des3.noeled(/1)) goto 1029
  843.             if (nligrp.ne.des3.noelep(/1)) goto 1029
  844.             do id1=1,nligrd
  845.               if (des2.noeled(id1).ne.des3.noeled(id1)) goto 1029
  846.               if (des2.lisdua(id1).ne.des3.lisdua(id1)) goto 1029
  847.             enddo
  848.             do id1=1,nligrp
  849.               if (des2.noelep(id1).ne.des3.noelep(id1)) goto 1029
  850.               if (des2.lisinc(id1).ne.des3.lisinc(id1)) goto 1029
  851.             enddo
  852. C
  853. C           Descripteur existant
  854.             nnhash(4)=nnhash(4)+1
  855.             jdes=ides
  856.             segsup,des2
  857.             goto 1032
  858. C
  859.  1029       continue
  860.             nnhash(3)=nnhash(3)+1
  861.  1030     continue
  862. C
  863. C         Nouveau descripteur
  864.           nnhash(3)=nnhash(3)+1
  865.           kdescr=kdescr+1
  866.           jdes=kdescr
  867.           if (lisde(/2).lt.kdescr) then
  868.             nbd=kdescr+nde0
  869.             segadj,lisde
  870.           endif
  871.           lisde(1,jdes)=des2
  872.           lisde(2,jdes)=idhash
  873.           lisde(3,jdes)=LTRK/(nligrd*nligrp)
  874.           goto 1033
  875. C
  876.  1032     continue
  877. C
  878. C         Rigidite avec memes descripteur et caracteristiques?
  879.           iptd=lisde(1,jdes)
  880.           do 1031 irie=krigel,1,-1
  881. C
  882.             if (ldesc(1,irie).ne.iptd)                  goto 1031
  883.             if (ldesc(3,irie).ne.mrigid.irigel(2,irig)) goto 1031
  884.             if (ldesc(4,irie).ne.mrigid.irigel(5,irig)) goto 1031
  885.             if (ldesc(5,irie).ne.mrigid.irigel(6,irig)) goto 1031
  886.             if (ldesc(6,irie).ne.mrigid.irigel(7,irig)) goto 1031
  887.             if (ldesc(2,irie).ge.lisde(3,jdes))         goto 1031
  888. C
  889.             jrig=irie
  890.             goto 1034
  891. C
  892.  1031     continue
  893. C
  894.  1033     continue
  895. C
  896. C         Nouvelle rigidite
  897.           krigel=krigel+1
  898.           jrig=krigel
  899.           if (ldesc(/2).lt.jrig) then
  900.             ndes=jrig+nde0
  901.             segadj,ldesc
  902.           endif
  903.           ldesc(1,jrig)=lisde(1,jdes)
  904.           ldesc(2,jrig)=0
  905.           ldesc(3,jrig)=mrigid.irigel(2,irig)
  906.           ldesc(4,jrig)=mrigid.irigel(5,irig)
  907.           ldesc(5,jrig)=mrigid.irigel(6,irig)
  908.           ldesc(6,jrig)=mrigid.irigel(7,irig)
  909. C
  910.  1034     continue
  911. C
  912.           lelem(1,iel)=ipt2
  913.           lelem(2,iel)=jrig
  914.           lelem(3,iel)=xmatr2
  915.           ldesc(2,jrig)=ldesc(2,jrig)+1
  916. C
  917.  1040   continue
  918. C       ===============================================================
  919. C
  920.         segdes xmatr4
  921.  1000 continue
  922. C     Fin de boucle sur les sous-matrices
  923. C     -----------------------------------
  924.       segsup,linco,xrela,lincp,lincd,trav3,trav3b
  925.       if (ldiag) then
  926.         write(ioimp,*)'Nb calculs hash          =',nnhash(1)
  927.         write(ioimp,*)'Nb comparaisons raccourcies hash=',nnhash(2)
  928.         write(ioimp,*)'Nb comparaisons raccourcies nohash=',nnhash(3)
  929.         write(ioimp,*)'Nb comparaisons completes=',nnhash(4)
  930.         write(ioimp,*)'Nb sous rigidite initiale=',mrigid.irigel(/2)
  931.         write(ioimp,*)'Nb sous rigidite finale  =',krigel
  932.         write(ioimp,*)'Nb de descripteurs       =',kdescr
  933.       endif
  934. C     =================================================================
  935. C     Construire la rigidite ou les inconnues ont ete eliminees
  936. C     =================================================================
  937.       nrigel=krigel
  938.       segini,ri4
  939.       ri4.mtymat='TEMPORAI'
  940.       ri4.iforig=mrigid.iforig
  941.       do ibr=1,nbrig
  942.         lelem=lrigi(ibr)
  943.         ipt4 =mrigid.irigel(1,ibr)
  944. C
  945. C       Renseigner les matrices elementaires
  946.         do 1042 iel=1,ipt4.num(/2)
  947.           ipt6=lelem(1,iel)
  948.           if (ipt6.eq.0) goto 1042
  949.           kdes=lelem(2,iel)
  950.           xmatr6=lelem(3,iel)
  951. C
  952.           itmp=1
  953.           bnew=.true.
  954.           if (xmatr6.lt.0) then
  955.             bnew=.false.
  956.             xmatr6=abs(xmatr6)
  957.             itmp=iel
  958.           endif
  959.           segact xmatr6
  960. C
  961.           des3=ldesc(1,kdes)
  962.           nelt=ldesc(2,kdes)
  963.           ipt5=ri4.irigel(1,kdes)
  964. C
  965. C         Creation de la sous-matrice
  966.           if (ipt5.eq.0) then
  967.             ldesc(2,kdes)=1
  968. C
  969.             nligrp=des3.noelep(/1)
  970.             nligrd=des3.noeled(/1)
  971.             nelrig=nelt
  972. C
  973.             bok=(nelrig.eq.xmatr6.re(/3).and.bnew)
  974.             if (bok) then
  975.               xmatr5=xmatr6
  976.               ipt5=ipt6
  977.               goto 1041
  978.             endif
  979. C
  980. **          write(6,*) 'segini xmatr5',nligrp,nligrd,nelrig
  981.             segini,xmatr5
  982.             xmatr5.symre=xmatr6.symre
  983. C
  984.             nbnn=ipt6.num(/1)
  985.             nbelem=nelt
  986.             nbref=0
  987.             nbsous=0
  988.             segini,ipt5
  989.             if (bnew) then
  990.               ipt5.itypel=28
  991.               if ((ipt6.itypel.eq.49).or.(ipt6.itypel.eq.22)) then
  992.                 ipt5.itypel=ipt6.itypel
  993.               endif
  994.             else
  995.               ipt5.itypel=ipt6.itypel
  996.             endif
  997. C
  998.  1041       continue
  999. C
  1000.             ri4.coerig(kdes)=mrigid.coerig(ibr)
  1001.             ri4.irigel(1,kdes)=ipt5
  1002.             ri4.irigel(2,kdes)=mrigid.irigel(2,ibr)
  1003.             ri4.irigel(3,kdes)=des3
  1004.             ri4.irigel(4,kdes)=xmatr5
  1005.             ri4.irigel(5,kdes)=mrigid.irigel(5,ibr)
  1006.             ri4.irigel(6,kdes)=mrigid.irigel(6,ibr)
  1007.             ri4.irigel(7,kdes)=mrigid.irigel(7,ibr)
  1008.             ri4.irigel(8,kdes)=mrigid.irigel(8,ibr)
  1009. C
  1010.             if (bok) goto 1042
  1011.           endif
  1012. C
  1013.           ielmt=ldesc(2,kdes)
  1014.           do inn=1,ipt6.num(/1)
  1015.             ipt5.num(inn,ielmt)=ipt6.num(inn,itmp)
  1016.           enddo
  1017. C
  1018.           xcoer=mrigid.coerig(ibr)/ri4.coerig(kdes)
  1019.           xmatr5=ri4.irigel(4,kdes)
  1020.           segact xmatr5*mod
  1021.           do ilc=1,xmatr5.re(/2)
  1022.             do ili=1,xmatr5.re(/1)
  1023.               xmatr5.re(ili,ilc,ielmt)=xcoer*xmatr6.re(ili,ilc,itmp)
  1024.             enddo
  1025.           enddo
  1026.           segdes xmatr5,xmatr6
  1027. C
  1028.           ldesc(2,kdes)=ldesc(2,kdes)+1
  1029.           if (bnew) segsup,ipt6,xmatr6
  1030.  1042   continue
  1031.         segsup,lelem
  1032.       enddo
  1033.       segsup,icpr,trav1,trav2,lrigi,ldesc
  1034. C
  1035. CCC   WRITE(*,*) 'MRIGID'
  1036. CCC   call prrigi(MRIGID,0)
  1037. CCC   WRITE(*,*) 'RI1'
  1038. CCC   call prrigi(ri1,0)
  1039. CCC   WRITE(*,*) 'RI4'
  1040. CCC   call prrigi(ri4,0)
  1041. C     =================================================================
  1042. C     Desactivation generale
  1043. C     =================================================================
  1044.       segact mrigid
  1045.       do irig=1,nbrig
  1046.         xmatri=irigel(4,irig)
  1047.         meleme=irigel(1,irig)
  1048.         descr=irigel(3,irig)
  1049.         segdes xmatri,descr
  1050.       enddo
  1051.       segact ri4
  1052.       do irig=1,ri4.irigel(/2)
  1053.         xmatri=ri4.irigel(4,irig)
  1054.         meleme=ri4.irigel(1,irig)
  1055.         descr=ri4.irigel(3,irig)
  1056.         segdes xmatri,descr
  1057.       enddo
  1058.       segact ri1
  1059.       do irig=1,nri1
  1060.         xmatri=ri1.irigel(4,irig)
  1061.         meleme=ri1.irigel(1,irig)
  1062.         descr=ri1.irigel(3,irig)
  1063.         segdes xmatri,descr
  1064.       enddo
  1065.       if(iimpi.ne.0.and.nbdeg.ne.0)
  1066.      >  write (6,*) ' nombre de relations degenerees'//
  1067.      >  ' apres elimination ',nbdeg
  1068. *     write(6,*) 'nbincl nligrp5 nnoeud ndes nlr',
  1069. *    >            nbincl,nligrp5,nnoeud,ndes,nlr
  1070.       end
  1071.  
  1072.  
  1073.  
  1074.  

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