Télécharger la plaquette Cast3M


fr en

Télécharger cmct4.eso

Retour à la liste

Numérotation des lignes : 

cmct4

  1. C CMCT4     SOURCE    GOUNAND   25/06/11    21:15:04     12278          
  2.       SUBROUTINE CMCT4(IRIGC,IRIGD,IRIGB,IRIGE,IMEL,IRIG2)
  3. *_______________________________________________________________________
  4. c
  5. c      opérateur KPRE (cmct option ELEM)
  6. c
  7. c entrée
  8. c   IRIGB : rigidité B
  9. c   IRIGC : rigidité C
  10. c   IRIGD : rigidité D
  11. c   IMEL  : maillage
  12. c
  13. c sortie
  14. c     IRIG2 : rigidité contenant la matrice C* D*-1 B*t construite par
  15. c     elements sur les elements de MELEME
  16. c     avec C*, D*, B* les matrices reduites sur les elements de MELEME
  17. c
  18. *_______________________________________________________________________
  19.  
  20.       IMPLICIT INTEGER(I-N)
  21.       IMPLICIT REAL*8(A-H,O-Z)
  22. *
  23.  
  24. -INC PPARAM
  25. -INC CCOPTIO
  26. -INC CCREEL
  27. -INC SMRIGID
  28.       POINTEUR IRIGC.MRIGID,IRIGD.MRIGID,IRIGB.MRIGID
  29.       POINTEUR IRIGC2.MRIGID,IRIGD2.MRIGID,IRIGB2.MRIGID
  30.       POINTEUR IRIGC3.MRIGID,IRIGB3.MRIGID
  31.       POINTEUR IRIGE.MRIGID,IRIGE2.MRIGID
  32.       POINTEUR IRIG2.MRIGID
  33.       POINTEUR XMATRC.XMATRI,XMATRD.XMATRI,XMATRB.XMATRI
  34.       POINTEUR XMATRE.XMATRI
  35.       POINTEUR XMAT.XMATRI
  36.       POINTEUR XMATV.XMATRI
  37.       SEGMENT vvec(ntot)
  38.       SEGMENT amat(ntot,ntot)
  39.       SEGMENT amats(ntot,ntot)
  40.       SEGMENT amatv(ntot,ntot)
  41. -INC SMELEME
  42. -INC CCHAMP
  43. -INC SMCOORD
  44. -INC SMCHPOI
  45. *
  46.       CHARACTER*(LOCHPO) NOMINC,NOMDUA
  47. *
  48.       segment icpr(nbpts)
  49.       segment inode(ino)
  50.       segment jelnum(imaxel,ino)
  51.       segment izone(imaxel,ino)
  52.       segment INFMEL
  53.       integer jmel,jcpr,jnode,kelnum,jzone
  54.       endsegment
  55.  
  56.  
  57.       logical lvdec,ldbg
  58. *
  59.       ldbg=.false.
  60.       if (ldbg )write(ioimp,*) 'irigc,irigd,irigb,irige=',irigc,irigd
  61.      $     ,irigb,irige
  62.       lvdec=.false.
  63.       nprint=40
  64.       iprint=0
  65.       npivhg=0
  66.       npivbd=0
  67.       nelpnu=0
  68.       xzpiv=xzprec**0.75d0
  69.       xzfac=xzprec**0.5d0
  70.       IF (IMEL.NE.0) THEN
  71.          CALL reduri(irigc,imel,irigc2)
  72.          if (ierr.ne.0) return
  73.          irigc=irigc2
  74.       ENDIF
  75. *
  76.       segact irigc
  77.       IFOC=IRIGC.IFORIG
  78.       IFO2=IFOC
  79.       nbsou=irigc.irigel(/2)
  80. C   Cas d'un MELEME vide => RIGIDITE vide
  81.       IF (nbsou.eq.0) THEN
  82.          NRIGEL = 0
  83.          SEGINI,IRIG2
  84.          IRIG2.MTYMAT = 'KPRE'
  85.          IRIG2.IFORIG = IFO2
  86.          RETURN
  87.       ENDIF
  88. C
  89. C     On commence par régulariser irigc au sens ou tous les points de
  90. C     son maillage sont atteints
  91. C
  92.       if (ldbg) write(ioimp,*) 'IRIGC->MELEME'
  93.       call rdscr2(irigc,irigc2)
  94.       IF (IERR.NE.0) RETURN
  95.       call melrig(irigc2,ipt5)
  96.       IF (IERR.NE.0) RETURN
  97.       call regmai(ipt5,meleme)
  98.       IF (IERR.NE.0) RETURN
  99. C
  100.       segact meleme
  101.       nbsou=lisous(/1)
  102.       SEGINI INFMEL
  103.       INFMEL.JMEL=MELEME
  104. *
  105. * Réduction de C sur le maillage MELEME
  106. *
  107.       if (ldbg) write(ioimp,*) ' IRIGC '
  108. *      call ecrobj('RIGIDITE',irigc)
  109. *      call ecrcha('RESU')
  110. *      call prlist
  111.       IRIGC3=0
  112.       CALL CMCT5(IRIGC2,INFMEL,IRIGC3)
  113.       IF (IERR.NE.0) RETURN
  114.       segsup irigc2
  115.       irigc2=irigc3
  116. *      write(ioimp,*) ' IRIGC2 '
  117. *      segprt,irigc2
  118. *      call ecrobj('RIGIDITE',irigc2)
  119. *      call ecrcha('RESU')
  120. *      call prlist
  121. *      return
  122. *
  123. * Réduction eventuelle de B sur le maillage MELEME et sur le DESCR de C2
  124. *
  125.       IF (IRIGB.NE.IRIGC) THEN
  126.          SEGACT IRIGB
  127.          IFOB=IRIGB.IFORIG
  128.          IF (IFOB.NE.IFOC) THEN
  129.             moterr(1:8)='RIGIDITE'
  130.             interr(1)=IFOB
  131.             interr(2)=IFOC
  132.             interr(3)=IFOUR
  133.             call erreur(1132)
  134.             IFO2 = IFOUR
  135.          ENDIF
  136.          call rdscr2(irigb,irigb2)
  137.          IF (IERR.NE.0) RETURN
  138. * Amélioration possible : transposition auto de B si necessaire
  139.          if (ldbg) write(ioimp,*) ' IRIGB .NE. IRIGC'
  140.          segini,irigb3=irigc2
  141. * On enleve les xmatri
  142.          nrig=irigb3.coerig(/1)
  143.          do irig=1,nrig
  144.             irigb3.irigel(4,irig)=0
  145.          enddo
  146.          CALL CMCT5(IRIGB2,INFMEL,IRIGB3)
  147.          IF (IERR.NE.0) RETURN
  148.          segsup irigb2
  149.          irigb2=irigb3
  150. *         write(ioimp,*) ' IRIGB2 '
  151. *     segprt,irigb2
  152. *         call ecrobj('RIGIDITE',irigb2)
  153. *         call ecrcha('RESU')
  154. *         call prlist
  155.       ELSE
  156.          if (ldbg) write(ioimp,*) ' IRIGB = IRIGC'
  157.       ENDIF
  158. *
  159. * Reduction de D sur le maillage MELEME et sur le DESCR primal de C2
  160. *
  161.       SEGACT IRIGD
  162.       IFOD=IRIGD.IFORIG
  163.       IF (IFOD.NE.IFOC) THEN
  164.          moterr(1:8)='RIGIDITE'
  165.          interr(1)=IFOD
  166.          interr(2)=IFOC
  167.          interr(3)=IFOUR
  168.          call erreur(1132)
  169.          IFO2 = IFOUR
  170.       ENDIF
  171.       NRIGEL=nbsou
  172.       SEGINI,IRIGD2
  173.       SEGACT,IRIGC2
  174.       DO irig=1,NRIGEL
  175.          ipt1=irigc2.irigel(1,irig)
  176.          des1=irigc2.irigel(3,irig)
  177.          call descar(des1,1,des3)
  178.          IF (IERR.NE.0) RETURN
  179.          irigd2.irigel(1,irig)=ipt1
  180.          irigd2.irigel(3,irig)=des3
  181.       ENDDO
  182.       if (ldbg) write(ioimp,*) ' IRIGD '
  183.       CALL CMCT6(IRIGD,IRIGD2)
  184.       IF (IERR.NE.0) RETURN
  185. *      call ecrobj('RIGIDITE',irigd2)
  186. *      call ecrcha('RESU')
  187. *      call prlist
  188. *      return
  189. *
  190. * Reduction (eventuelle) de E sur le maillage MELEME et sur le DESCR dual de C2
  191. *
  192.       IF (IRIGE.NE.0) THEN
  193.          SEGACT IRIGD
  194.          IFOE=IRIGE.IFORIG
  195.          IF (IFOE.NE.IFOC) THEN
  196.             moterr(1:8)='RIGIDITE'
  197.             interr(1)=IFOE
  198.             interr(2)=IFOC
  199.             interr(3)=IFOUR
  200.             call erreur(1132)
  201.             IFO2 = IFOUR
  202.          ENDIF
  203.          NRIGEL=nbsou
  204.          SEGINI,IRIGE2
  205.          SEGACT,IRIGC2
  206.          DO irig=1,NRIGEL
  207.             ipt1=irigc2.irigel(1,irig)
  208.             des1=irigc2.irigel(3,irig)
  209.             call descar(des1,2,des3)
  210.             IF (IERR.NE.0) RETURN
  211.             irige2.irigel(1,irig)=ipt1
  212.             irige2.irigel(3,irig)=des3
  213.          ENDDO
  214.          if (ldbg) write(ioimp,*) ' IRIGE '
  215.          CALL CMCT6(IRIGE,IRIGE2)
  216.          IF (IERR.NE.0) RETURN
  217.       ELSE
  218.          IRIGE2=0
  219.       ENDIF
  220. *
  221. * On est pret maintenant à calculer E - C D-1 Bt
  222. *     On envisage trois algorithmes :
  223. *     1) Force Brute : inversion de D puis produit
  224. *     2) Factorisation de la matrice  (D Bt) = (Lu   0 ) ( Uu  Upu)
  225. *                                     (C 0 )   (Lup  Lp) ( 0   Up )
  226. *      Si je comprends bien, Lp * Up = C D-1 Bt Oui !
  227. *     3) Brute brute : d'après la formulation d'inversion, on inverse
  228. *     tout, puis on inverse le bloc pp de l'inverse. Ca nous redonne le
  229. *     complément de Schur
  230. *
  231. *
  232. *     Algorithme factorisation LDMt: on pourrait pivoter symétriquement ?
  233. *     Repris de Golub et Van Loan !
  234. *
  235.       NRIGEL=nbsou
  236.       SEGINI IRIG2
  237.       IRIG2.MTYMAT = 'KPRE'
  238.       IRIG2.IFORIG = IFO2
  239.       SEGACT IRIGC2
  240.       DO irig=1,NRIGEL
  241.          IRIG2.COERIG(irig)=1.D0
  242.          ipt1=irigc2.irigel(1,irig)
  243.          des1=irigc2.irigel(3,irig)
  244.          nligd=des1.noeled(/1)
  245.          nligp=des1.noelep(/1)
  246.          call descar(des1,2,des3)
  247.          IF (IERR.NE.0) RETURN
  248.          nlrig=ipt1.num(/2)
  249.          nelpnu=nelpnu+nlrig
  250. * Matrice de stockage totale
  251.          isym=2
  252.          IF (IRIGB.EQ.IRIGC) THEN
  253.             IF (irigd2.irigel(7,irig).eq.0) THEN
  254.                if (irige2.ne.0) then
  255.                   isym=irige2.irigel(7,irig)
  256.                else
  257.                   isym=0
  258.                endif
  259.             endif
  260.          endif
  261.          call rdscr1(ipt1,des3,isym,ipt2,des2)
  262.          IF (IERR.NE.0) RETURN
  263.          irig2.irigel(1,irig)=ipt2
  264.          irig2.irigel(3,irig)=des2
  265.          IRIG2.IRIGEL(7,irig)=isym
  266. *
  267. * Remplissage des matrices élémentaires
  268. *
  269.          nelrig=nlrig
  270.          nligrp=des3.noelep(/1)
  271.          nligrd=des3.noeled(/1)
  272.          segini xmatr3
  273.          xmatr3.symre=isym
  274.          xmatrd=irigd2.irigel(4,irig)
  275. *         segprt,xmatrd
  276.          xmatrc=irigc2.irigel(4,irig)
  277.          if (irigb.ne.irigc) then
  278.             xmatrb=irigb2.irigel(4,irig)
  279.          else
  280.             xmatrb=xmatrc
  281.          endif
  282.          if (irige.ne.0) then
  283.             xmatre=irige2.irigel(4,irig)
  284.          else
  285.             xmatre=0
  286.          endif
  287.          ntot=nligp+nligd
  288.          segini amat
  289.          if (lvdec) then
  290.             segini amats
  291.             segini amatv
  292.          endif
  293.          segini vvec
  294.          do ilrig=1,nlrig
  295. *     Copies
  296.             do iligp=1,nligp
  297.                do iligd=1,nligp
  298. *                  segprt,xmatrd
  299. *                  write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd
  300.                   amat(iligd,iligp)=xmatrd.re(iligd,iligp,ilrig)
  301. *                  write(ioimp,*) 'amat(iligd,iligp)=',amat(iligp,iligd)
  302.                enddo
  303.             enddo
  304.             do iligp=1,nligp
  305.                do iligd=1,nligd
  306. *                  write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd
  307.                   amat(iligd+nligp,iligp)=
  308.      $                 xmatrc.re(iligd,iligp,ilrig)
  309. *                  write(ioimp,*) 'amat(iligd+nligp,iligp)=',amat(iligd
  310. *     $                 +nligp,iligp)
  311.                enddo
  312.             enddo
  313.             do iligp=1,nligp
  314.                do iligd=1,nligd
  315. *                  write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd
  316.                   amat(iligp,iligd+nligp)=
  317.      $                 xmatrb.re(iligd,iligp,ilrig)
  318. *                  write(ioimp,*) 'amat(iligp,iligd+nligp)=',amat(iligp
  319. *     $                 ,iligd+nligp)
  320.                enddo
  321.             enddo
  322.             if (xmatre.ne.0) then
  323.                do iligp=1,nligd
  324.                   do iligd=1,nligd
  325.                      amat(iligp+nligp,iligd+nligp)=
  326.      $                    xmatre.re(iligd,iligp,ilrig)
  327.                   enddo
  328.                enddo
  329.             else
  330.                do iligp=1,nligd
  331.                   do iligd=1,nligd
  332.                      amat(iligp+nligp,iligd+nligp)=0.D0
  333.                   enddo
  334.                enddo
  335.             endif
  336. * Sauvegarde de A dans As
  337.             if (lvdec) then
  338.                do j=1,ntot
  339.                   do i=1,ntot
  340.                      amats(i,j)=amat(i,j)
  341.                   enddo
  342.                enddo
  343. *               segprt,amats
  344.             endif
  345. *            if (ilrig.eq.1) segprt,amats
  346. *  Echelle par la partie haut gauche
  347.             xscahg=0.d0
  348.             do j=1,nligp
  349.                xscahg=max(xscahg,amat(j,j))
  350.             enddo
  351. * xpethg positi
  352.             xpethg=xscahg*xzpiv
  353.             if (xpethg.EQ.0.d0) then
  354.                write(ioimp,*) 'D is a null matrix'
  355.                segprt,des3
  356.                write(ioimp,*) 'irig,nrigel=',irig,nrigel
  357.                write(ioimp,*) 'ilrig,nlrig=',ilrig,nlrig
  358.                write(ioimp,*) 'nligp,nligd=',nligp,nligd
  359.                segprt,amat
  360.                call erreur(185)
  361.                return
  362.             endif
  363. *
  364.             xscabd=0.d0
  365.             do j=1,nligp
  366.                xest=0.d0
  367.                xpiv=max(xpethg,amat(j,j))
  368.                do i=1,nligd
  369.                   xest=xest-amat(i+nligp,j)*amat(j,i+nligp)/xpiv
  370.                enddo
  371.                xscabd=min(xscabd,xest)
  372.             enddo
  373. * xpetbd negatif
  374.             xpetbd=xscabd*xzpiv
  375.             if (xpetbd.EQ.0.d0) then
  376.                write(ioimp,*) 'B or C are null matrices'
  377.                segprt,des3
  378.                write(ioimp,*) 'irig,nrigel=',irig,nrigel
  379.                write(ioimp,*) 'ilrig,nlrig=',ilrig,nlrig
  380.                write(ioimp,*) 'nligp,nligd=',nligp,nligd
  381.                segprt,amat
  382.                call erreur(185)
  383.                return
  384.             endif
  385.             if (lvdec) then
  386.                xtol=max(xscahg,abs(xscabd))*xzfac
  387.             endif
  388. *
  389. *     Factorisation de amat (LDMt ou LDLt d'apres Golub et Van Loan 3eme
  390. *     edition pp. 137-139)
  391. * Amélioration possible : prendre en compte la symétrie
  392.             do j=1,ntot
  393. *     Résoudre L(1:1j,1:j) v(1:j) = A(1:j,j)
  394.                do k=1,j
  395.                   vvec(k)=amat(k,j)
  396.                enddo
  397.                do k=1,j-1
  398.                   do i=k+1,j
  399.                      vvec(i)=vvec(i)-vvec(k)*amat(i,k)
  400.                   enddo
  401.                enddo
  402. *     Calculer M(j,1:j-1) et stocker dans A (1:j-1,j)
  403.                do i=1,j-1
  404.                   amat(i,j)=vvec(i)/amat(i,i)
  405.                enddo
  406. *     Stocker d(j) dans A(j,j)
  407.                amat(j,j)=vvec(j)
  408.                if (j.le.nligp) then
  409.                   if (abs(amat(j,j)).le.xpethg) then
  410.                      if (ldbg) then
  411.                         write(ioimp,*) 'Pivot petit detecte ilrig,j='
  412.      $                       ,ilrig,j
  413.                         write(ioimp,*) '| ',amat(j,j),' | < ',xpethg
  414.                         write(ioimp,*) 'Noeud ',ipt1.num(des1.noelep(j)
  415.      $                       ,ilrig), ' Inconnue ',des1.lisinc(j)
  416.                         segprt,amats
  417.                         write(ioimp,*) 'Diag hg :'
  418.                         write(ioimp,*) (amats(jj,jj),jj=1,nligp)
  419.                         write(ioimp,*) 'Diag bd :'
  420.                         write(ioimp,*) (amats(jj,jj),jj=nligp+1,ntot)
  421.                         write(ioimp,*) 'Pivots hg :'
  422.                         write(ioimp,*) (amat(jj,jj),jj=1,nligp)
  423.                         write(ioimp,*) 'Pivots bd :'
  424.                         write(ioimp,*) (amat(jj,jj),jj=nligp+1,ntot)
  425.                      endif
  426.                      amat(j,j)=xpethg
  427.                      npivhg=npivhg+1
  428.                   endif
  429.                else
  430.                   if (abs(amat(j,j)).le.-xpetbd) then
  431.                      if (ldbg) then
  432.                         write(ioimp,*) 'Pivot petit detecte ilrig,j='
  433.      $                       ,ilrig,j
  434.                         write(ioimp,*) '| ',amat(j,j),' | < -1* ',xpetbd
  435.                         write(ioimp,*) 'Noeud ',ipt1.num(des3.noelep(j
  436.      $                       -nligp),ilrig), ' Inconnue ',des3.lisinc(j
  437.      $                       -nligp)
  438.                         segprt,amats
  439.                         write(ioimp,*) 'Diag hg :'
  440.                         write(ioimp,*) (amats(jj,jj),jj=1,nligp)
  441.                         write(ioimp,*) 'Diag bd :'
  442.                         write(ioimp,*) (amats(jj,jj),jj=nligp+1,ntot)
  443.                         write(ioimp,*) 'Pivots hg :'
  444.                         write(ioimp,*) (amat(jj,jj),jj=1,nligp)
  445.                         write(ioimp,*) 'Pivots bd :'
  446.                         write(ioimp,*) (amat(jj,jj),jj=nligp+1,ntot)
  447.                      endif
  448.                      amat(j,j)=xpetbd
  449.                      npivbd=npivbd+1
  450.                   endif
  451.                endif
  452. *     Calculer L(j+1:n,j) et stocker dans A (j+1:n,j)
  453.                do k=1,j-1
  454.                   do i=j+1,ntot
  455.                      amat(i,j)=amat(i,j)-vvec(k)*amat(i,k)
  456.                   enddo
  457.                enddo
  458. *               write(ioimp,*) 'ilrig,j,a(j,j),xpetbd,xpethg=',ilrig,j
  459. *     $              ,amat(j,j),xpetbd,xpethg
  460.                do i=j+1,ntot
  461.                   amat(i,j)=amat(i,j)/amat(j,j)
  462.                enddo
  463.             enddo
  464. *     Vérification de la factorisation
  465.             if (lvdec) then
  466. *               segprt,amat
  467.                do k=1,ntot
  468.                   do i=1,ntot
  469.                      amatv(i,k)=0.D0
  470.                   enddo
  471.                enddo
  472.                do k=1,ntot
  473.                   do i=1,ntot
  474.                      do j=1,min(i,k)
  475.                         xajj=amat(j,j)
  476.                         if (k.eq.j) then
  477.                            xukj=1.d0
  478.                         else
  479.                            xukj=amat(j,k)
  480.                         endif
  481.                         if (i.eq.j) then
  482.                            xlij=1.d0
  483.                         else
  484.                            xlij=amat(i,j)
  485.                         endif
  486.                         amatv(i,k)=amatv(i,k)+xlij*xajj*xukj
  487.                      enddo
  488.                   enddo
  489.                enddo
  490. *               segprt,amatv
  491.                do k=1,ntot
  492.                   do i=1,ntot
  493.                      xda=abs(amatv(i,k)-amats(i,k))
  494.                      if (xda.gt.xtol) then
  495.                         iprint=iprint+1
  496.                         write(ioimp,*)
  497.      $                       'Erreur factorisation detectee ilrig='
  498.      $                       ,ilrig
  499.                         write(ioimp,*) 'aik=',amats(i,k),'aikv=',amatv(i
  500.      $                       ,k)
  501.                         write(ioimp,*) 'i=',i,' k=',k
  502.                         write(ioimp,*) 'xpethg,xpetbd,xtol=',xpethg
  503.      $                       ,xpetbd,xtol
  504.                      write(ioimp,*) 'Diag hg orig :'
  505.                      write(ioimp,*) (amats(jj,jj),jj=1,nligp)
  506.                      write(ioimp,*) 'Diag bd orig :'
  507.                      write(ioimp,*) (amats(jj,jj),jj=nligp+1,ntot)
  508.                      write(ioimp,*) 'Diag hg :'
  509.                      write(ioimp,*) (amatv(jj,jj),jj=1,nligp)
  510.                      write(ioimp,*) 'Diag bd :'
  511.                      write(ioimp,*) (amatv(jj,jj),jj=nligp+1,ntot)
  512.                      segprt,amats
  513.                      segprt,amatv
  514.                      ierr=1
  515.                      return
  516. *                        write(ioimp,*) 'xda = ',xda,' > ',xtol
  517. *      write(ioimp,*) 'Noeud dual i ',ipt1.num(des3.noeled(i)
  518. *     $                    ,ilrig), ' Inconnue ',des3.lisdua(i)
  519. *      write(ioimp,*) 'Noeud primal k',ipt1.num(des3.noelep(k)
  520. *     $                    ,ilrig), ' Inconnue ',des3.lisinc(k)
  521.                         if (iprint.gt.nprint) then
  522.                            call erreur(185)
  523.                            return
  524.                         endif
  525.                      endif
  526.                   enddo
  527.                enddo
  528.             endif
  529. * Produit LDMt avec le bloc en bas à droite et stockage dans RE
  530. * Amélioration possible : prendre en compte la symétrie
  531.             do k=1,nligd
  532.                k2 = nligp+k
  533.                do i=1,nligd
  534.                   i2 = nligp+i
  535.                   do j=1,min(i,k)
  536.                      j2 = nligp+j
  537.                      xajj=amat(j2,j2)
  538.                      if (k.eq.j) then
  539.                         xukj=1.d0
  540.                      else
  541.                         xukj=amat(j2,k2)
  542.                      endif
  543.                      if (i.eq.j) then
  544.                         xlij=1.d0
  545.                      else
  546.                         xlij=amat(i2,j2)
  547.                      endif
  548.                      xmatr3.re(i,k,ilrig)=xmatr3.re(i,k,ilrig)
  549.      $                    +xlij*xajj*xukj
  550.                   enddo
  551.                enddo
  552.             enddo
  553.          ENDDO
  554.          segsup vvec
  555.          if (lvdec) then
  556.             segsup amats
  557.             segsup amatv
  558.          endif
  559.          segsup amat
  560.          IRIG2.IRIGEL(4,irig)=xmatr3
  561.          IRIG2.IRIGEL(5,irig)=IRIGC2.IRIGEL(5,irig)
  562.       enddo
  563. *     Impression pivots nuls
  564.       npivnu=npivhg+npivbd
  565.       if (npivnu.gt.0.and.iimpi.ne.0) then
  566. *      if (npivnu.gt.0) then
  567.          xpivel=float(npivnu)/float(nelpnu)
  568.          xpilhg=float(npivhg)/float(nelpnu)
  569.          xpilbd=float(npivbd)/float(nelpnu)
  570.          write(ioimp,*) npivnu,' pivots nuls detectes / ',nelpnu,
  571.      $        ' elements = ',xpivel,' pivots nuls/element'
  572.          write(ioimp,*) '  dont : ', xpilhg
  573.      $        ,' pivots nuls bloc(1,1)/element'
  574.          write(ioimp,*) '         ', xpilbd
  575.      $        ,' pivots nuls bloc(2,2)/element'
  576.       endif
  577. *      segprt,irig2
  578. *
  579. * Menage
  580. *
  581.       DO irig=1,nbsou
  582. *         descr=irigc2.irigel(3,irig)
  583. *         segsup descr
  584.          xmatri=irigc2.irigel(4,irig)
  585.          segsup xmatri
  586.          IF (IRIGB.NE.IRIGC) THEN
  587. *            descr=irigb2.irigel(3,irig)
  588. *            segsup descr
  589.             xmatri=irigb2.irigel(4,irig)
  590.             segsup xmatri
  591.          ENDIF
  592. *         descr=irigd2.irigel(3,irig)
  593. *        segsup descr
  594.          xmatri=irigd2.irigel(4,irig)
  595.          segsup xmatri
  596.          IF (IRIGE.NE.0) THEN
  597. *            descr=irige2.irigel(3,irig)
  598. *            segsup descr
  599.             xmatri=irige2.irigel(4,irig)
  600.             segsup xmatri
  601.          ENDIF
  602.       ENDDO
  603.       segsup irigc2
  604.       IF (IRIGB.NE.IRIGC) segsup irigb2
  605.       segsup irigd2
  606.       IF (IRIGE.NE.0) segsup irige2
  607. *
  608.       icpr=jcpr
  609.       inode=jnode
  610.       jelnum=kelnum
  611.       izone=jzone
  612.       segsup izone
  613.       segsup jelnum
  614.       segsup inode
  615.       segsup icpr
  616.       SEGSUP INFMEL
  617. *
  618.       segsup meleme
  619.       RETURN
  620.       END
  621.  
  622.  

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