Télécharger la plaquette Cast3M


fr en

Télécharger scnd2.eso

Retour à la liste

Numérotation des lignes : 

scnd2

  1. C SCND2     SOURCE    OF166741  24/05/02    21:15:04     11927          
  2. *
  3.       subroutine scnd2(mrigid,ri1,ri4)
  4. *
  5. *    entree
  6. *         mrigid rigidité  sans sous matrices de dependances
  7. *         ri1 rigidite  de dependance
  8. *
  9. *    sortie
  10. *         ri4     rigidite condensee les noeuds dependants sont
  11. *                 elimines
  12. *
  13. *    on suppose que la premiere inconnue de la raideur de dependance
  14. *    (celle qu'on elimine) n'apparait que dans cette dependance
  15. *    Cela va permettre de constituer les tableau :
  16. *    noeu-inconnue (primale ou duale) ==> numéro raideur de dépendance
  17. *    au moment de la création des matrices, on va s'efforcer de factoriser
  18. *    les descripteurs
  19. *
  20.       IMPLICIT INTEGER(I-N)
  21.       IMPLICIT REAL*8 (A-H,O-Z)
  22. *
  23. -INC SMRIGID
  24. -INC SMELEME
  25.  
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC SMCOORD
  29. -INC CCHAMP
  30. -INC CCREEL
  31.       segment icpr(nbno)
  32.       segment trav1
  33.       character*4 compd(nbincl),compdd(nbincl)
  34.       endsegment
  35.       segment trav2
  36.       integer irigd(nbnoc,nbincd),nelrigd(nbnoc,nbincd)
  37.       endsegment
  38.       segment trav3
  39.       integer cor1p(nligrp),cor1d(nligrd)
  40.       integer imatt(nelrig)
  41.       endsegment
  42.       segment trav4
  43.       integer lisdes(klisd)
  44.       endsegment
  45.       segment trav5
  46.       real*8 xmatmpv(nligpv)
  47.       integer imatpv(nligpv)
  48.       endsegment
  49.       segment imatno(lmatno)
  50.       segment idesno(ldesno)
  51.       segment igeomno(lgeomno)
  52.       segment coerno(lcoerno)
  53.       segment irg7no(lrg7no)
  54.       segment irg6no(lrg6no)
  55.       segment irg5no(lrg5no)
  56.       segment irg2no(lrg2no)
  57.       logical dupli,descar,dupli4
  58.       character*4 co1,co1d
  59.       segment svu
  60.       logical vu(lrvu)
  61.       endsegment
  62.       integer lvu,symr
  63. *
  64. *      write (6,*) ' raideur initiale '
  65. *      call prrigi(mrigid,0)
  66. *      write (6,*) ' raideur dependance '
  67. *      call prrigi(ri1,1)
  68.  
  69. *  meme valeur de test que dans chole
  70.       xzref=5.d-12        
  71.       lrvu=100
  72.       lmatno=1000
  73.       ldesno=1000
  74.       lgeomno=1000
  75.       lcoerno=1000
  76.       lrg7no=1000
  77.       lrg6no=1000
  78.       lrg5no=1000
  79.       lrg2no=1000
  80.       kmatno=0
  81.       kdesno=0
  82.       kgeomno=0
  83.       kcoerno=0
  84.       krg7no=0
  85.       krg6no=0
  86.       krg5no=0
  87.       krg2no=0
  88.       segini imatno,idesno,igeomno,coerno,irg7no,irg5no,svu
  89.       lvu=0
  90.       segini irg2no,irg6no
  91.       nbdeg=0
  92. *  remplissage du segment trav1 (liste des noms de composantes)
  93.       nbincd=0
  94.       nbincl=16
  95.       segini trav1
  96.       nbno=nbpts
  97.       segact ri1
  98.       do 10 irig=1,ri1.irigel(/2)
  99.          descr=ri1.irigel(3,irig)
  100.          segact descr
  101.          do 40 nligrd=1,lisdua(/2)
  102.             do 50 inc=1,nbincd
  103.                if (compd(inc).eq.lisdua(nligrd)) goto 40
  104.  50         continue
  105.             nbincd=nbincd+1
  106.             if(nbincd.gt.nbincl) then
  107.              nbincl=2*nbincd
  108.              segadj trav1
  109.             endif
  110.             compd(nbincd)=lisdua(nligrd)
  111.             compdd(nbincd)=' '
  112. **          write(6,*) 'inc-1 compdd ',inc,compd(inc),compdd(inc)
  113.  40      continue
  114.          ldua=lisdua(/2)
  115.          do 45 nligrd=1,lisinc(/2)
  116.             do 55 inc=1,nbincd
  117.                if (compd(inc).eq.lisinc(nligrd)) goto 46
  118.  55         continue
  119.             nbincd=nbincd+1
  120.             inc=nbincd
  121.             if(nbincd.gt.nbincl) then
  122.              nbincl=2*nbincd
  123.              segadj trav1
  124.             endif
  125.             compd(nbincd)=lisinc(nligrd)
  126.  46      continue
  127.  
  128. ** si la matrice n'est pas carre, avec quoi on rempli?
  129.             compdd(inc)=lisdua(min(ldua,nligrd))
  130. **          write(6,*) 'inc-2 compdd ',inc,compd(inc),compdd(inc)
  131.  45      continue
  132.  10   continue
  133. *      write (6,*) ' scnd2 - 1 ',(compd(i),i=1,nbincd)
  134.       segact mrigid
  135.       do 110 irig=1,irigel(/2)
  136.          descr=irigel(3,irig)
  137.          segact descr
  138. **       write(6,*) 'nligrp,nligrd ',lisinc(/2),lisdua(/2)
  139.          ldua=lisdua(/2)
  140.          do 140 nligrd=1,lisinc(/2)
  141.             do 150 inc=1,nbincd
  142.                if (compd(inc).eq.lisinc(nligrd)) goto 141
  143.  150        continue
  144.             nbincd=nbincd+1
  145.             inc=nbincd
  146.             if(nbincd.gt.nbincl) then
  147.              nbincl=2*nbincd
  148.              segadj trav1
  149.             endif
  150.             compd(nbincd)=lisinc(nligrd)
  151.  141     continue
  152. *  au  cas ou la primale n'est pas dans nomdd, on sauve la duale
  153. **            write(6,*) 'inc lisdua descr',inc,lisdua(nligrd),descr
  154. ** si la matrice n'est pas carre, avec quoi on rempli la duale?
  155.             compdd(inc)=lisdua(min(ldua,nligrd))
  156. **          write(6,*) 'inc-3 compdd ',inc,compd(inc),compdd(inc)
  157.  140     continue
  158.  110  continue
  159. *      write (6,*) ' scnd2 - 2 ',(compd(i),i=1,nbincd)
  160. *  au tour de trav2 (inconnue => raideur de dépendance)
  161. *  mais d'abord icpr
  162.       segini icpr
  163.       nbnoc=0
  164.       do irig=1,ri1.irigel(/2)
  165.          meleme=ri1.irigel(1,irig)
  166.          descr=ri1.irigel(3,irig)
  167.          segact meleme,descr
  168.          do  j=1,num(/2)
  169.             ip=num(noeled(1),j)
  170.             if(icpr(ip).eq.0) then
  171.               nbnoc=nbnoc+1
  172.               icpr(ip)=nbnoc
  173.             endif
  174.          enddo
  175.       enddo
  176.       do irig=1,irigel(/2)
  177.          meleme=irigel(1,irig)
  178.          segact meleme
  179.          do  j=1,num(/2)
  180.          do  i=1,num(/1)
  181.             ip=num(i,j) 
  182.             if(icpr(ip).eq.0) then
  183.               nbnoc=nbnoc+1
  184.               icpr(ip)=nbnoc
  185.             endif
  186.          enddo
  187.          enddo
  188.       enddo
  189.       segini trav2
  190.       do 200 irig=1,ri1.irigel(/2)
  191.          meleme=ri1.irigel(1,irig)
  192.          descr=ri1.irigel(3,irig)
  193.          xmatri=ri1.irigel(4,irig)
  194.          segact xmatri
  195.          do 250 incd=1,nbincd
  196.             if (compd(incd).eq.lisdua(1)) goto 260
  197.  250     continue
  198.          call erreur(5)
  199.  260     continue
  200.          do 240 i=1,num(/2)
  201.             ip=num(noeled(1),i)
  202.             irigd(icpr(ip),incd)=irig
  203.             nelrigd(icpr(ip),incd)=i
  204.  240     continue
  205.  200  continue
  206. *  rechercher les noms d'inconnues duales
  207.       do 300 i=1,compd(/2)
  208.          do 310 j=1,lnomdd
  209.             if (compd(i).eq.nomdd(j)) goto 320
  210.  310     continue
  211. **          write (6,*) ' pas de nom dual trouvé ',compd(i),compdd(i)
  212.          goto 300
  213.  320     continue
  214.             compdd(i)=nomdu(j)
  215. **        write (6,*) ' correspondance ', compd(i), compdd(i)
  216.  300  continue
  217. *        write (6,*) ' scnd2 - 3 ',compdd(/2),(compdd(i),i=1,compdd(/2))
  218. *
  219. *  maintenant, on balaye les raideurs. Pour chacune d'entre elle, on
  220. *  balaye les inconnues pour trouver une dépendance. Si oui on construit
  221. *  le résultat de l'opération
  222.       segini,ri4=mrigid
  223. *     write(6,*) 'ri4 mrigid ',ri4,mrigid
  224.       ri4.ichole=0
  225.       ri4.imgeo1=0
  226.       ri4.imgeo2=0
  227.       ri4.isupeq=0
  228. *  boucle raideur en cours de condensation
  229. *  on range dans des2 le descripteur precedent pour le reutiliser si possible
  230.       des2=0
  231.       isup=0
  232.       klisd=100
  233.       llisd=0
  234.       segini trav4
  235. *     call prrigi(ri4,1)
  236.       do 1000 irig=1,ri4.irigel(/2)
  237. *         write(ioimp,*) 'Balayage raideur ',irig
  238. ** on ne dupliquera meleme et xmatri que si necessaire
  239.          meleme=ri4.irigel(1,irig)
  240.          ipt4=meleme
  241.          dupli4=.false.
  242.          ri4.irigel(1,irig)=ipt4
  243.          des4=ri4.irigel(3,irig)
  244.          segact des4
  245.          Xmatri=ri4.irigel(4,irig)
  246.          segact xmatri
  247. *  au cas ou symre ait ete oublie
  248.          if (symre.ne.ri4.irigel(7,irig)) then
  249.            write(6,*) 'incoherence symre irigel(7) pour ',ri4,xmatri
  250.            segact xmatri*mod
  251.            symre = max(symre,ri4.irigel(7,irig))
  252.            ri4.irigel(7,irig)=symre
  253.          endif
  254.          Xmatr4=Xmatri
  255. *        write(6,*) 'scnd2 xmatri xmatr4 dims ',xmatri,xmatr4,
  256. *    >     xmatr4.re(/1),xmatr4.re(/2),xmatr4.re(/3)
  257.          nelrig=xmatr4.re(/3)
  258.          ri4.irigel(4,irig)=Xmatr4
  259. * On regarde le descripteur. S'il est symétrique, on pourra
  260. * faire le travail en une passe sinon deux.
  261. *   Correspondance lisinc, lisdua, numéro incd dans trav1
  262.          nligrp=des4.lisinc(/2)
  263.          nligrd=des4.lisdua(/2)
  264.          segini trav3
  265.          nligpv=nligrp*2+256
  266.          segini trav5
  267.          do 410 iligrp=1,nligrp
  268.             do 420 incd=1,nbincd
  269.                if (compd(incd).eq.des4.lisinc(iligrp)) then
  270.                   cor1p(iligrp)=incd
  271.                   goto 410
  272.                endif
  273.  420        continue
  274. **          write(6,*) 'cor1p pas de correspondance'
  275. **          call erreur(5)
  276.  410     continue
  277.          do 510 iligrd=1,nligrd
  278.             do 520 incd=1,nbincd
  279. **             write(6,*) 'compdd ',compdd(incd),des4.lisdua(iligrd)
  280.                if (compdd(incd).eq.des4.lisdua(iligrd)) then
  281.                   cor1d(iligrd)=incd
  282.                   goto 510
  283.                endif
  284.  520        continue
  285. **          write(6,*) 'cor1d pas de correspondance descar',descar
  286. **          call erreur(5)
  287.  510     continue
  288. *         write(ioimp,*) 'nligrp=',nligrp
  289. *         write(ioimp,*) 'nligrd=',nligrd
  290.          descar=.false.
  291.          if (nligrp.ne.nligrd) goto 600
  292.          do iligr=1,nligrp
  293.             if (cor1p(iligr).ne.cor1d(iligr)) goto 600
  294.             if (des4.noelep(iligr).ne.des4.noeled(iligr)) goto 600
  295.          enddo
  296.          if (xmatr4.symre.eq.0) descar=.true.
  297.  600     continue
  298. *         WRITE(IOIMP,194)(des4.NOELEP(J),des4.LISINC(J),cor1p(j),J=1
  299. *     $        ,nligrp)
  300. *         WRITE(IOIMP,194)(des4.NOELED(J),des4.LISDUA(J),cor1d(j),J=1
  301. *     $        ,nligrd)
  302. *  194     FORMAT( I6,9X,A4,9X,I6)
  303. *         if (descar) then
  304. *            write(ioimp,*) 'le descripteur est carre'
  305. *         else
  306. *            write(ioimp,*) 'le descripteur n''est pas carre'
  307. *         endif
  308.          npass=2
  309.          idsup=0
  310.          do 1040 i=1,ipt4.num(/2)
  311.             dupli=.false.
  312.             renorm=0.d0
  313. *    On fait deux passes : une pour la primale et une pour la duale
  314.             do 1045 ipass=1,npass
  315.                if (ipass.eq.1) then
  316.                   nligr4=des4.lisinc(/2)
  317.                else
  318.                   nligr4=des4.lisdua(/2)
  319.                endif
  320.                do 1010 nligri=1,nligr4
  321.                   if (ipass.eq.1) then
  322.                      incd=cor1p(nligri)
  323.                   else
  324.                      incd=cor1d(nligri)
  325.                   endif
  326. ***               write(6,*) 'ipass incd ',ipass,incd
  327.                   if (incd.eq.0) goto 1010
  328.                   if (ipass.eq.1) then
  329.                      ip=ipt4.num(des4.noelep(nligri),i)
  330.                   else
  331.                      ip=ipt4.num(des4.noeled(nligri),i)
  332.                   endif
  333. *  y a t'il une dependance concernant ce noeud
  334. *        write(6,*) ' descar',descar
  335. ***                 write(6,*) 'ipass ip icpr incd nelrig',ipass,ip,
  336. **   >              icpr(ip),incd,nelrigd(icpr(ip),incd)
  337.                   if (nelrigd(icpr(ip),incd).eq.0) then
  338.                     goto 1010
  339.                   endif
  340. *  ici il faut eliminer l'inconnue
  341. *  on commence par dupliquer la la matrice si ce n'est pas deja fait
  342.                   if (.not.dupli) then
  343.                      ird=irigel(/1)
  344.                      krg7no=krg7no+1
  345.                      if (krg7no.gt.irg7no(/1)) then
  346.                        lrg7no=irg7no(/1)+1000
  347.                        segadj irg7no
  348.                      endif
  349.                      if (ird.ge.7) then
  350.                         irg7no(krg7no)=ri4.irigel(7,irig)
  351.                      else
  352.                         irg7no(krg7no)=0
  353.                      endif
  354.                      krg6no=krg6no+1
  355.                      if (krg6no.gt.irg6no(/1)) then
  356.                        lrg6no=irg6no(/1)+1000
  357.                        segadj irg6no
  358.                      endif
  359.                      if (ird.ge.6) then
  360.                         irg6no(krg6no)=ri4.irigel(6,irig)
  361.                      else
  362.                         irg6no(krg6no)=0
  363.                      endif
  364.                      krg5no=krg5no+1
  365.                      if (krg5no.gt.irg5no(/1)) then
  366.                        lrg5no=irg5no(/1)+1000
  367.                        segadj irg5no
  368.                      endif
  369.                      if (ird.ge.5) then
  370.                         irg5no(krg5no)=ri4.irigel(5,irig)
  371.                      else
  372.                         irg5no(krg5no)=0
  373.                      endif
  374.                      krg2no=krg2no+1
  375.                      if (krg2no.gt.irg2no(/1)) then
  376.                        lrg2no=irg2no(/1)+1000
  377.                        segadj irg2no
  378.                      endif
  379.                      irg2no(krg2no)=ri4.irigel(2,irig)
  380.                      kdesno=kdesno+1
  381.                      if (kdesno.gt.idesno(/1)) then
  382.                        ldesno=idesno(/1)+1000
  383.                        segadj idesno
  384.                      endif
  385.                      nbnn=ipt4.num(/1)
  386.                      nbelem=1
  387.                      nbref=0
  388.                      nbsous=0
  389.                      segini meleme
  390.                      itypel=28
  391. *  si c'est une matrice de relation, il faut garder l'attribut
  392. **                   if(ipt4.itypel.eq.49) write(6,*)
  393. **   >                  'ipt4 dans scnd2 ',ipt4
  394.                      if (ipt4.itypel.eq.49) itypel=49
  395.                      if (ipt4.itypel.eq.22) itypel=22
  396.                      do i1=1,nbnn
  397.                         num(i1,1)=ipt4.num(i1,i)
  398.                      enddo
  399.                      kgeomno=kgeomno+1
  400.                      if (kgeomno.gt.igeomno(/1)) then
  401.                        lgeomno=igeomno(/1)+1000
  402.                        segadj igeomno
  403.                      endif
  404.                      igeomno(kgeomno)=meleme
  405.                      kligrp=xmatr4.re(/2)
  406.                      kligrd=xmatr4.re(/1)
  407.                   irig1=irigd(icpr(ip),incd)
  408.                   des1=ri1.irigel(3,irig1)
  409.                   segact des1
  410. *  on agrandit xmatri du nb d'inc. primales presentent dans la relation
  411.                   nligrp1=max(des1.lisinc(/2),0)
  412.                      nligrp=kligrp+nligrp1
  413.                      nligrd=kligrd+nligrp1
  414.                      nelrig=1
  415.                      segini xmatri
  416.                      xmatri.symre=xmatr4.symre
  417. *  re est plus grand que xmatr4.re donc pas d'appel a scndc possible
  418.                      renorm=0.d0
  419.                      ii=1
  420.                      if(itypel.eq.22) ii=2
  421.                      if(itypel.eq.49) ii=3
  422.                      do io=1,kligrp
  423.                         do iu=1,kligrd
  424.                            retmp=xmatr4.re(iu,io,i)
  425.                            re(iu,io,1)=retmp
  426.                            if (io.ge.ii.or.iu.ge.ii) renorm=renorm+
  427.      >                                                abs(retmp)
  428.                         enddo
  429.                      enddo
  430.                      renorm=max(xpetit,renorm*xzref)
  431.                      segini,descr
  432.                      idesno(kdesno)=descr
  433.                      do io=1,des4.noelep(/1)
  434.                         noelep(io)=des4.noelep(io)
  435.                         lisinc(io)=des4.lisinc(io)
  436.                      enddo
  437.                      do io=1,des4.noeled(/1)
  438.                         noeled(io)=des4.noeled(io)
  439.                         lisdua(io)=des4.lisdua(io)
  440.                      enddo
  441.  
  442.                      kmatno=kmatno+1
  443.                      if (kmatno.gt.imatno(/1))then
  444.                        lmatno=imatno(/1)+1000
  445.                        segadj imatno
  446.                      endif
  447.                      imatno(kmatno)=xmatri
  448.                      kcoerno=kcoerno+1
  449.                      if (kcoerno.gt.coerno(/1))then
  450.                        lcoerno=coerno(/1)+1000
  451.                        segadj coerno
  452.                      endif
  453.                      coerno(kcoerno)=ri4.coerig(irig)
  454.                      dupli=.true.
  455.                      imatt(i)=1
  456.                   endif
  457. *  on elimine dans la matrice dupliquee
  458.                   irig1=irigd(icpr(ip),incd)
  459.                   i1=nelrigd(icpr(ip),incd)
  460.                   ipt1=ri1.irigel(1,irig1)
  461.                   des1=ri1.irigel(3,irig1)
  462.                   xmatr1=ri1.irigel(4,irig1)
  463.                   segact ipt1,des1,xmatr1
  464. *          nligrd1 =des1.lisdua(/2)
  465. *   nligrd1=1 est l'inconnue à éliminer (de la matrice de dependance)
  466.                   nligrd1=1
  467.                   nligrp1=des1.lisinc(/2)
  468. *          write (6,*) ' scnd2  nligrp1  nligrd1 i1 ',nligrp1,nligrd1,i1,
  469. *     > ip,incp
  470. *  nlig1 : inconnue primale à remettre dans la matrice (comme primale)
  471.                   ipvm=0
  472.                   do 1050 nlig1=1,nligrp1
  473.                      xmatmp = xmatr1.re(1,nlig1,i1)
  474.                      ip1=ipt1.num(des1.noelep(nlig1),i1)
  475.                      co1=des1.lisinc(nlig1)
  476.                      if (ipass.eq.2.or.descar) then
  477.                         do  incd=1,nbincd
  478.                            if (compd(incd).eq.co1) then
  479.                               co1d=compdd(incd)
  480.                               goto 1055
  481.                            endif
  482.                         enddo
  483.                         write (6,*) ' scnd2 2 pas de correspondance '
  484.      $                       ,co1
  485.                         call erreur(5)
  486.                         return
  487.  1055                   continue
  488.                      endif
  489. *           write (6,*) ' ip1 co1 nligrpi ',ip1,co1,nligrpi
  490. *  le ddl est il deja dans la raideur
  491.                      nbno=0
  492.                      if (ipass.eq.1) then
  493.                         do 1060 nligr=1,kligrp
  494.                            if (noelep(nligr).eq.0) goto 1060
  495.                            if (num(noelep(nligr),1).ne.ip1) goto 1060
  496. * on note le numéro local du noeud pour le reutiliser
  497.                            nbno=noelep(nligr)
  498.                            if (lisinc(nligr).ne.co1) goto 1060
  499.                            goto 1070
  500.  1060                   continue
  501.                      else
  502.                         do 2060 nligr=1,kligrd
  503.                            if (noeled(nligr).eq.0) goto 2060
  504.                            if (num(noeled(nligr),1).ne.ip1) goto 2060
  505. * on note le numéro local du noeud pour le reutiliser
  506.                            nbno=noeled(nligr)
  507.                            if (lisdua(nligr).ne.co1d) goto 2060
  508.                            goto 1070
  509.  2060                   continue
  510.                      endif
  511. *  on n'a pas trouve le ddl, on le rajoute
  512. *            write (6,*) ' ajout de ddl'
  513.                      if (descar) then
  514.                         nligr =kligrp+1
  515.                         kligrp=kligrp+1
  516.                         kligrd=kligrd+1
  517.                      else
  518.                         if (ipass.eq.1) then
  519.                            nligr =kligrp+1
  520.                            kligrp=kligrp+1
  521.                            kligrd=kligrd
  522.                            xmatri.symre=2
  523.                         else
  524.                            nligr =kligrd+1
  525.                            kligrp=kligrp
  526.                            kligrd=kligrd+1
  527.                            xmatri.symre=2
  528.                         endif
  529.                      endif
  530.                      nligrp=lisinc(/2)
  531.                      nligrd=lisdua(/2)
  532.                      if (kligrp.gt.nligrp.or.kligrd.gt.nligrd) then
  533.                        nligrd=kligrd+max(1,nligrp1)
  534.                        nligrp=kligrp+max(1,nligrp1)
  535. *                      write(6,*) ' premier adj xmatri '
  536.                        segadj xmatri
  537.                        segadj descr
  538.                      endif
  539.                      if (nbno.eq.0) then
  540.                         nbsous=0
  541.                         nbelem=1
  542.                         nbnn=num(/1)+1
  543.                         nbref=0
  544.                         segadj meleme
  545.                         nbno=nbnn
  546.                      endif
  547.                      num(nbno,1)=ip1
  548.                      if (ipass.eq.1.or.descar) then
  549.                         noelep(kligrp)=nbno
  550.                         lisinc(kligrp)=co1
  551.                      endif
  552.                      if (ipass.eq.2.or.descar) then
  553.                         lisdua(kligrd)=co1d
  554.                         noeled(kligrd)=nbno
  555.                      endif
  556.  1070                continue
  557. *  on met a jour la raideur
  558.                      nligrp=kligrp
  559.                      nligrd=kligrd
  560. *            segact xmatr1
  561.                         if(nligr.gt.nligpv) then
  562.                           nligpv=nligr*2
  563.                           segadj trav5
  564. **                        write(6,*) 'agrandissement trav5 ',nligpv
  565.                         endif
  566.                      if (ipass.eq.2) then
  567.                         ipvm=ipvm+1
  568.                         imatpv(ipvm)=nligr
  569.                         xmatmpv(ipvm)=xmatmp
  570.                      endif
  571.                      if (ipass.eq.1) then
  572. *  unrolling sur 2
  573.                         do j=1,nligrd-1,2
  574.                            retmp1=re(j,nligr,1)+re(j,nligri,1)
  575.      $                          *xmatmp
  576.                            re(j,nligr,1)=retmp1
  577.                            retmp2=re(j+1,nligr,1)+re(j+1,nligri,1)
  578.      $                          *xmatmp
  579.                            re(j+1,nligr,1)=retmp2
  580.                         enddo
  581.                         if (j.eq.nligrd) then
  582.                            retmp1=re(j,nligr,1)+re(j,nligri,1)
  583.      $                          *xmatmp
  584.                            re(j,nligr,1)=retmp1
  585.                         endif
  586.                      endif
  587. *    la mise a jour est sorti de la boucle pour ppouvoir l'effectuer dans l'ordre efficace pour le cache
  588. ***                  if (ipass.eq.2) then
  589. ***                     do j=1,nligrp
  590. ***                        re(nligr,j,1)=re(nligr,j,1)+re(nligri,j,1)
  591. ***  $                          *xmatmp
  592. ***                   enddo
  593. ***                endif
  594. *              segdes xmatr1,des1
  595.  1050             continue
  596. *  mise a jour transferee ici
  597.                  if(ipass.eq.2) then
  598. *  unrolling sur 2
  599.                   do j=1,nligrp
  600.                    retmp=re(nligri,j,1)
  601.                    do ipv=1,ipvm-1,2    
  602.                     nligr1=imatpv(ipv)
  603.                     re(nligr1,j,1)=re(nligr1,j,1)+retmp*xmatmpv(ipv)
  604.                     nligr2=imatpv(ipv+1)
  605.                     re(nligr2,j,1)=re(nligr2,j,1)+retmp*xmatmpv(ipv+1)
  606.                    enddo
  607.                    if (ipv.eq.ipvm) then
  608.                        nligr=imatpv(ipv)
  609.                        re(nligr,j,1)=re(nligr,j,1)+retmp*xmatmpv(ipv)
  610.                    endif
  611.                   enddo
  612.                  endif
  613. *  on note les inconnues à supprimer (primale et duale)
  614. *           if (ipass.eq.2.or.descar)
  615. *    >      write (6,*) ' il faut eliminer ',nligri,lisdua(nligri)
  616. *      write (6,*) ' lisinc contient ',(lisinc(iou),iou=1,lisinc(/2))
  617.                   if (ipass.eq.1.or.descar) noelep(nligri)=0
  618.                   if (ipass.eq.2.or.descar) noeled(nligri)=0
  619.  1010          continue
  620. *
  621. *            write(ioimp,*) 'Avant suppression des inconnues flaguees'
  622. *          WRITE(IOIMP,194)(NOELEP(J),LISINC(J),J=1,LISINC(/2))
  623. *          WRITE(IOIMP,194)(NOELED(J),LISDUA(J),J=1,LISDUA(/2))
  624. *  194     FORMAT( I6,9X,A4)
  625. *          WRITE(IOIMP,195) ((RE(L,K,1),K=1,RE(/2)),L=1,RE(/1))
  626. * 195      FORMAT(1X,6E12.5)
  627. *  on supprime les inconnues flaguées
  628.  1045       continue
  629.             if (dupli) then
  630.                nligrd=kligrd
  631.                nligrp=kligrp
  632.                idecp=0
  633.                do 1080 nligr=1,nligrp
  634.                   if (noelep(nligr).eq.0) then
  635.                      idecp=idecp+1
  636.                   elseif (idecp.ne.0) then
  637.                      lisinc(nligr-idecp)=lisinc(nligr)
  638.                      noelep(nligr-idecp)=noelep(nligr)
  639.                      do nligd=1,nligrd
  640.                         if (noeled(nligd).ne.0)
  641.      >                  re(nligd,nligr-idecp,1)=re(nligd,nligr,1)
  642.                      enddo
  643.                   endif
  644.  1080          continue
  645. *pv  on imbrique la boucle 1093 dans ce sens pour minimiser l'utilisation du cache
  646.                nligrp=nligrp-idecp
  647.              do 1093 nligr=1,nligrp
  648.                idecd=0
  649.                    do nligd=1,nligrd
  650.                 if (noeled(nligd).eq.0) then
  651.                    idecd=idecd+1
  652.                 elseif (idecd.ne.0) then
  653.                       re(nligd-idecd,nligr,1)=re(nligd,nligr,1)
  654.                 endif
  655.                    enddo
  656. 1093          continue
  657.              idecd=0
  658.              do 1092 nligd=1,nligrd
  659.                 if (noeled(nligd).eq.0) then
  660.                    idecd=idecd+1
  661.                 elseif (idecd.ne.0) then
  662.                    lisdua(nligd-idecd)=lisdua(nligd)
  663.                    noeled(nligd-idecd)=noeled(nligd)
  664.                 endif
  665. 1092         continue
  666.                nligrd=nligrd-idecd
  667. *            write (6,*) ' scnd2 avant apres   ',lisinc(/2),nligrp
  668.           if(descar.and..true.) then
  669. *  symetriser xmatri
  670.           ninc=32
  671.           do ib=1,nligrd,ninc
  672.             do jr=1,nligrd 
  673.               do ir=ib,min(jr-1,ib+ninc-1)
  674.                 re1=re(ir,jr,1)
  675.                 re2=re(jr,ir,1)
  676.                    rem=(re1+re2)/2.d0
  677.                    re(ir,jr,1)=rem
  678.                    re(jr,ir,1)=rem
  679.               enddo
  680.             enddo
  681.           enddo
  682.           endif
  683.  
  684. *  on elimine les matrices ne contenant que les multiplicateurs de Lagrange
  685. *  ou les autres termes non significatifs  
  686. *  la verif du second membre nul est faite dans vechpo
  687. *
  688. **             if(itypel.eq.49) write(6,*) 'scnd2 49 ',meleme
  689.                if (itypel.eq.49.or.itypel.eq.22) then
  690. ***          write (6,*)  (lisinc(iou),iou=1,lisinc(/2))
  691.                    iid=2
  692.                    if (itypel.eq.49) iid=3
  693. *  suppression des termes negligeables de la relation
  694.                    idec=0
  695. *  on utilise une valeur absolu car la relation initiale est supposee normalisee a 1
  696.                    do ii=iid,nligrd
  697.                      if (abs(re(ii,1,1)).lt.renorm.and.
  698.      >                   abs(re(1,ii,1)).lt.renorm) then
  699.                       idec=idec+1
  700.                      elseif(idec.ne.0) then
  701.                         re(ii-idec,1,1)=re(ii,1,1)
  702.                         if(iid.eq.3) re(ii-idec,2,1)=re(ii,2,1)
  703.                         re(1,ii-idec,1)=re(1,ii,1)
  704.                         if(iid.eq.3) re(2,ii-idec,1)=re(2,ii,1)
  705.                         lisinc(ii-idec)=lisinc(ii)
  706.                         lisdua(ii-idec)=lisdua(ii)
  707.                         noelep(ii-idec)=noelep(ii)
  708.                         noeled(ii-idec)=noeled(ii)
  709.                      endif
  710.                    enddo
  711.                    nligrp=nligrp-idec
  712.                    nligrd=nligrd-idec
  713.                    if(nligrd.lt.iid) then
  714. **           write (6,*) 'scnd2 matrice supprimee LX',
  715. **   >       nligrd,itypel,meleme
  716.                      nligrd=0
  717.                      nligrp=0
  718.                      nbdeg=nbdeg+1
  719.                    endif
  720.                endif
  721.                if (re(/1).ne.nligrd.or.re(/2).ne.nligrp) then
  722. **             write(6,*) ' ajustement final xmatri ',re(/1),re(/2),
  723. **   >         nligrd,nligrp
  724.                segadj xmatri
  725.                else
  726. **              write(6,*) ' pas d ajustement xmatri '
  727.                endif
  728.                if (nligrp.ne.nligrd.or..not.descar) then
  729.                  xmatri.symre=2
  730.                else
  731.                endif
  732.                segdes xmatri
  733.                segadj descr
  734. * la matrice est fabriquée, on elimine maintenant les noeuds
  735. * non references
  736.                lvu=num(/1)
  737.                if (lvu.gt.lrvu) then
  738.                 lrvu=lvu+100
  739.                 segadj svu
  740.                endif
  741.                do in=1,lvu
  742.                   vu(in)=.false.
  743.                enddo
  744.                do in=1,nligrp
  745.                   vu(noelep(in))=.true.
  746.                enddo
  747.                if (.not.descar) then
  748.                   do in=1,nligrd
  749.                      vu(noeled(in))=.true.
  750.                   enddo
  751.                endif
  752.                idec=0
  753.                do 1095 in=1,num(/1)
  754.                   if (.not.vu(in)) then
  755. *  noeud a eliminer
  756.                      idec=idec+1
  757.                   elseif (idec.ne.0) then
  758.                      do no=1,nligrp
  759.                         if (noelep(no).eq.in) then
  760.                            noelep(no)=noelep(no)-idec
  761.                            if (descar) noeled(no)=noeled(no)-idec
  762.                         endif
  763.                      enddo
  764.                      if (.not.descar) then
  765.                         do no=1,nligrd
  766.                            if (noeled(no).eq.in) then
  767.                               noeled(no)=noeled(no)-idec
  768.                            endif
  769.                         enddo
  770.                      endif
  771.                      num(in-idec,1)=num(in,1)
  772.                   endif
  773.  1095          continue
  774.                nbnn=num(/1)-idec
  775.                nbelem=1
  776.                nbsous=0
  777.                nbref=0
  778.                segadj meleme
  779.                if (nbnn.eq.0) imatno(kmatno)=0
  780. *  recyclage eventuel du descripteur
  781.                do 4032 ides=llisd,1,-1
  782.                des2=lisdes(ides)
  783.                segact des2
  784.          if (nligrd.ne.des2.noeled(/1)) goto 4030
  785.          if (nligrp.ne.des2.noelep(/1)) goto 4030
  786.          do id1=1,nligrd
  787.            if (noeled(id1).ne.des2.noeled(id1)) goto 4030
  788.            if (lisdua(id1).ne.des2.lisdua(id1)) goto 4030
  789.          enddo
  790.          do id1=1,nligrp
  791.            if (noelep(id1).ne.des2.noelep(id1)) goto 4030
  792.            if (lisinc(id1).ne.des2.lisinc(id1)) goto 4030
  793.          enddo
  794.          segsup descr
  795.          idsup=idsup+1
  796.          descr=des2
  797.          idesno(kdesno)=descr
  798.          goto 4031
  799.  4030    continue
  800.  4032    continue
  801.                 llisd=llisd+1
  802.                 if (llisd.gt.klisd) then
  803.                   klisd=llisd+100
  804.                   segadj trav4
  805.                 endif
  806.                 lisdes(llisd)=descr
  807. *               write(6,*) ' nouveau descripteur ',nligrp,nligrd
  808. **              segdes des2
  809.  4031    continue
  810. **             endif
  811.                if (nligrp.ne.noelep(/1).or.nligrd.ne.noeled(/1))
  812.      >             segadj descr
  813.             endif
  814.  1040    continue
  815. *  compression de ipt4 et imatr4
  816.          idec=0
  817.          ipt5=ipt4
  818.          xmatr5=xmatr4
  819.          nligrp=xmatr5.re(/2)
  820.          nligrd=xmatr5.re(/1)
  821.          nelrig=xmatr5.re(/3)
  822.          nbnn=ipt5.num(/1)
  823.          nbelem=ipt5.num(/2)
  824.          idec=0
  825.          do i=1,xmatr4.re(/3)
  826.            if (imatt(i).ne.0) idec=idec+1
  827.          enddo
  828. *  on duplique systematiquement pour pouvoir detruire dans resouc
  829.          if (idec.ne.0.or..true.) then
  830.                  nbelem=nbelem-idec
  831.                  nbsous=0
  832.                  nbref=ipt5.lisref(/1)
  833.                  segini,ipt4
  834.                  ipt4.itypel=ipt5.itypel
  835.                  do ir=1,ipt4.lisref(/1)
  836.                   ipt4.lisref(ir)=ipt5.lisref(ir)
  837.                  enddo
  838.                  ri4.irigel(1,irig)=ipt4
  839.                  nelrig=nelrig-idec
  840.                  segini,xmatr4
  841.                  ri4.irigel(4,irig)=Xmatr4
  842.                  xmatr4.symre=xmatr5.symre
  843.                  idec=0
  844.          endif
  845.          do 2000 i=1,xmatr5.re(/3)
  846.             if(imatt(i).ne.0) then
  847.                idec=idec+1
  848.             elseif (xmatr4.ne.xmatr5) then
  849. **             do io=1,nligrp
  850. **                do iu=1,nligrd
  851. **                   xmatr4.re(iu,io,i-idec)=xmatr5.re(iu,io,i)
  852. **                enddo
  853. **             enddo
  854.                      call scndc(nligrp*nligrd,xmatr5.re(1,1,i),
  855.      >               xmatr4.re(1,1,i-idec))
  856.                do in=1,nbnn
  857.                    ipt4.num(in,i-idec)=ipt5.num(in,i)
  858.                enddo
  859.             endif
  860.  2000    continue
  861.          segdes xmatr4
  862. *        write (6,*) ' nelrig dans scnd2 ',nelrig
  863.          if (nelrig.eq.0) ri4.irigel(4,irig)=0
  864.          segsup trav3
  865.          segsup trav5
  866. **       segdes descr
  867.  1000 continue
  868.  
  869. *
  870. *  compression des xmatri suivant les descripteurs
  871. *
  872.        isup=0
  873.        id=0
  874.        is=1
  875.        nid=1
  876.        nrigf=0
  877.  4100  id=nid
  878.        if (id.gt.kdesno) goto 4120
  879.        xmatri=imatno(id)
  880.        if (xmatri.eq.0) then
  881. *        write(6,*) ' xmatri 0',id
  882.          isup=isup+1
  883.          nid=id+1
  884.          goto 4100
  885.        endif
  886.        meleme=igeomno(id)
  887.        segact meleme
  888.        des4=idesno(id)
  889.        nid=0
  890.        nelrig=num(/2) 
  891. *      if (num(/2).ne.1) write(6,*) 'nelrig-1',nelrig
  892.        nelrii=nelrig
  893.        nrigf=nrigf+1
  894.        lis=id
  895.        do is=id+1,kdesno
  896.          if (idesno(is).ne.des4) goto 4110
  897.          if (imatno(is).eq.0) goto 4110
  898.          if (irg2no(is).ne.irg2no(id)) goto 4110
  899.          if (irg5no(is).ne.irg5no(id)) goto 4110
  900.          if (irg6no(is).ne.irg6no(id)) goto 4110
  901.          if (irg7no(is).ne.irg7no(id)) goto 4110
  902.          lis=is
  903.          meleme=igeomno(is)
  904.          segact meleme
  905.          nelrig=nelrig+num(/2)
  906. *      if (num(/2).ne.1) write(6,*) 'nelrig-2',num(/2)
  907.          goto 4112
  908.  4110   continue
  909.         if (nid.eq.0) nid=is
  910.  4112   continue
  911.        enddo
  912.        if (nid.eq.0) nid=is
  913.        if (nelrig.eq.nelrii) goto 4100
  914.        xmatri=imatno(id)
  915.        segact xmatri
  916.        symr=symre
  917.        nligrd=re(/1)
  918.        nligrp=re(/2)
  919.        segini xmatri
  920.        symre=symr
  921.        meleme=igeomno(id)
  922.        segact meleme
  923.        nbnn=num(/1)
  924.        nbelem=nelrig
  925.        nbsous=0
  926.        nbref=0
  927.        itype=itypel
  928.        segini meleme
  929.        itypel=itype
  930.        nelrig=0
  931.        do is=id,lis
  932.         coeref=1.d0
  933.         if (is.ne.id) then
  934.          if (idesno(is).ne.des4) goto 4111
  935.          if (imatno(is).eq.0) goto 4111
  936.          if (irg2no(is).ne.irg2no(id)) goto 4111
  937.          if (irg5no(is).ne.irg5no(id)) goto 4111
  938.          if (irg6no(is).ne.irg6no(id)) goto 4111
  939.          if (irg7no(is).ne.irg7no(id)) goto 4111
  940.          coeref=coerno(is)/coerno(id)
  941.         endif
  942.          xmatr1=imatno(is)
  943.          ipt1=igeomno(is)
  944.          segact xmatr1,ipt1
  945.          symre=max(symre,xmatr1.symre)
  946.          do nel=1,xmatr1.re(/3)
  947.           nelrig=nelrig+1
  948.           do irp=1,nligrp
  949.           do ird=1,nligrd
  950.            re(ird,irp,nelrig)=xmatr1.re(ird,irp,nel)*coeref
  951.           enddo
  952.           enddo
  953.           do ip=1,nbnn
  954.            num(ip,nelrig)=ipt1.num(ip,nel)
  955.           enddo
  956.           imatno(is)=0
  957.          enddo
  958.         segsup xmatr1,ipt1
  959.         isup=isup+1
  960.  4111   continue
  961.        enddo
  962.        imatno(id)=xmatri
  963.        igeomno(id)=meleme
  964.        segdes xmatri
  965.  
  966.        goto 4100
  967.  4120  continue
  968. *      write(6,*) ' suppression de raideurs et descr ',isup,idsup
  969. *      write(6,*) ' nb desc          ',llisd
  970. **     write(6,*) ' lisdes ',(lisdes(il),il=1,lisdes(/1))
  971.          segsup trav4
  972. *  compression de irigel
  973.       idec=0
  974.  
  975.       do 2100 irig=1,ri4.irigel(/2)
  976.          if (ri4.irigel(4,irig).eq.0.or.abs(coerig(irig)).lt.xpetit) 
  977.      >      then
  978.             idec=idec+1
  979.          elseif (idec.ne.0) then
  980.             do in=1,ri4.irigel(/1)
  981.                ri4.irigel(in,irig-idec)=ri4.irigel(in,irig)
  982.             enddo
  983.             ri4.coerig(irig-idec)=ri4.coerig(irig)
  984.          endif
  985.  2100 continue
  986.       nrigeo=ri4.irigel(/2)-idec
  987. *  il ne reste plus qu'a ranger dans ri4 les nouvelles raideurs engendrées
  988.       nrigel=nrigf+nrigeo
  989.       segadj ri4
  990.       ir=nrigeo
  991.       do 3000 irig=1,kmatno
  992.          if (imatno(irig).eq.0) goto 3000
  993.          if (abs(coerno(irig)).lt.xpetit) goto 3000
  994.          ir=ir+1
  995.          ri4.coerig(ir)=coerno(irig)
  996.          ri4.irigel(1,ir)=igeomno(irig)
  997.          ri4.irigel(3,ir)=idesno(irig)
  998.          ri4.irigel(7,ir)=irg7no(irig)
  999.          ri4.irigel(5,ir)=irg5no(irig)
  1000.          ri4.irigel(2,ir)=irg2no(irig)
  1001.          ri4.irigel(6,ir)=irg6no(irig)
  1002.          ri4.irigel(4,ir)=imatno(irig)
  1003.  3000 continue
  1004. *     write(6,*) ' ir nrigel ',ir,nrigel,nrigeo
  1005. *  il peut y avoir des imatri nul a enlever
  1006.       if (ir.lt.nrigel) then
  1007.        nrigel=ir
  1008.        segadj ri4
  1009.       endif
  1010.       if (ir.ne.nrigel) call erreur(5)
  1011. *     write (6,*) ' raideur resultante '
  1012. *       call prrigi(ri4,0)
  1013.       segsup trav1,trav2,imatno,idesno,igeomno,coerno,irg7no,icpr
  1014.      $     ,irg5no,irg2no,irg6no,svu
  1015. *  desactivation generale
  1016.       segact mrigid
  1017.       do irig=1,irigel(/2)
  1018.          xmatri=irigel(4,irig)
  1019.          meleme=irigel(1,irig)
  1020.          descr=irigel(3,irig)
  1021.          segdes xmatri,descr
  1022.       enddo
  1023.       segact ri4
  1024.       do irig=1,ri4.irigel(/2)
  1025.          xmatri=ri4.irigel(4,irig)
  1026.          meleme=ri4.irigel(1,irig)
  1027.          descr=ri4.irigel(3,irig)
  1028.          segdes xmatri,descr
  1029.       enddo
  1030.       segact ri1
  1031.       do irig=1,ri1.irigel(/2)
  1032.          xmatri=ri1.irigel(4,irig)
  1033.          meleme=ri1.irigel(1,irig)
  1034.          descr=ri1.irigel(3,irig)
  1035.          segdes xmatri,descr
  1036.       enddo
  1037.       if(iimpi.ne.0.and.nbdeg.ne.0)
  1038.      >  write (6,*) ' nombre de relations degenerees'//
  1039.      >  ' apres elimination ',nbdeg
  1040.       end
  1041.  
  1042.  
  1043.  
  1044.  
  1045.  
  1046.  
  1047.  
  1048.  
  1049.  
  1050.  
  1051.  
  1052.  
  1053.  
  1054.  
  1055.  
  1056.  
  1057.  
  1058.  
  1059.  
  1060.  
  1061.  
  1062.  
  1063.  
  1064.  
  1065.  
  1066.  
  1067.  
  1068.  
  1069.  
  1070.  
  1071.  
  1072.  
  1073.  
  1074.  
  1075.  

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