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    PV090527  26/04/30    21:15:19     12529          
  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 IRIGE
  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.          rigrel=0
  273.          segini xmatr3
  274.          xmatr3.symre=isym
  275.          xmatrd=irigd2.irigel(4,irig)
  276. *         segprt,xmatrd
  277.          xmatrc=irigc2.irigel(4,irig)
  278.          if (irigb.ne.irigc) then
  279.             xmatrb=irigb2.irigel(4,irig)
  280.          else
  281.             xmatrb=xmatrc
  282.          endif
  283.          if (irige.ne.0) then
  284.             xmatre=irige2.irigel(4,irig)
  285.          else
  286.             xmatre=0
  287.          endif
  288.          ntot=nligp+nligd
  289.          segini amat
  290.          if (lvdec) then
  291.             segini amats
  292.             segini amatv
  293.          endif
  294.          segini vvec
  295.          do ilrig=1,nlrig
  296. *     Copies
  297.             do iligp=1,nligp
  298.                do iligd=1,nligp
  299. *                  segprt,xmatrd
  300. *                  write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd
  301.                   amat(iligd,iligp)=xmatrd.re(iligd,iligp,ilrig)
  302. *                  write(ioimp,*) 'amat(iligd,iligp)=',amat(iligp,iligd)
  303.                enddo
  304.             enddo
  305.             do iligp=1,nligp
  306.                do iligd=1,nligd
  307. *                  write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd
  308.                   amat(iligd+nligp,iligp)=
  309.      $                 xmatrc.re(iligd,iligp,ilrig)
  310. *                  write(ioimp,*) 'amat(iligd+nligp,iligp)=',amat(iligd
  311. *     $                 +nligp,iligp)
  312.                enddo
  313.             enddo
  314.             do iligp=1,nligp
  315.                do iligd=1,nligd
  316. *                  write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd
  317.                   amat(iligp,iligd+nligp)=
  318.      $                 xmatrb.re(iligd,iligp,ilrig)
  319. *                  write(ioimp,*) 'amat(iligp,iligd+nligp)=',amat(iligp
  320. *     $                 ,iligd+nligp)
  321.                enddo
  322.             enddo
  323.             if (xmatre.ne.0) then
  324.                do iligp=1,nligd
  325.                   do iligd=1,nligd
  326.                      amat(iligp+nligp,iligd+nligp)=
  327.      $                    xmatre.re(iligd,iligp,ilrig)
  328.                   enddo
  329.                enddo
  330.             else
  331.                do iligp=1,nligd
  332.                   do iligd=1,nligd
  333.                      amat(iligp+nligp,iligd+nligp)=0.D0
  334.                   enddo
  335.                enddo
  336.             endif
  337. * Sauvegarde de A dans As
  338.             if (lvdec) then
  339.                do j=1,ntot
  340.                   do i=1,ntot
  341.                      amats(i,j)=amat(i,j)
  342.                   enddo
  343.                enddo
  344. *               segprt,amats
  345.             endif
  346. *            if (ilrig.eq.1) segprt,amats
  347. *  Echelle par la partie haut gauche
  348.             xscahg=0.d0
  349.             do j=1,nligp
  350.                xscahg=max(xscahg,amat(j,j))
  351.             enddo
  352. * xpethg positi
  353.             xpethg=xscahg*xzpiv
  354.             if (xpethg.EQ.0.d0) then
  355.                write(ioimp,*) 'D is a null matrix'
  356.                segprt,des3
  357.                write(ioimp,*) 'irig,nrigel=',irig,nrigel
  358.                write(ioimp,*) 'ilrig,nlrig=',ilrig,nlrig
  359.                write(ioimp,*) 'nligp,nligd=',nligp,nligd
  360.                segprt,amat
  361.                call erreur(185)
  362.                return
  363.             endif
  364. *
  365.             xscabd=0.d0
  366.             do j=1,nligp
  367.                xest=0.d0
  368.                xpiv=max(xpethg,amat(j,j))
  369.                do i=1,nligd
  370.                   xest=xest-amat(i+nligp,j)*amat(j,i+nligp)/xpiv
  371.                enddo
  372.                xscabd=min(xscabd,xest)
  373.             enddo
  374. * xpetbd negatif
  375.             xpetbd=xscabd*xzpiv
  376.             if (xpetbd.EQ.0.d0) then
  377.                write(ioimp,*) 'B or C are null matrices'
  378.                segprt,des3
  379.                write(ioimp,*) 'irig,nrigel=',irig,nrigel
  380.                write(ioimp,*) 'ilrig,nlrig=',ilrig,nlrig
  381.                write(ioimp,*) 'nligp,nligd=',nligp,nligd
  382.                segprt,amat
  383.                call erreur(185)
  384.                return
  385.             endif
  386.             if (lvdec) then
  387.                xtol=max(xscahg,abs(xscabd))*xzfac
  388.             endif
  389. *
  390. *     Factorisation de amat (LDMt ou LDLt d'apres Golub et Van Loan 3eme
  391. *     edition pp. 137-139)
  392. * Amélioration possible : prendre en compte la symétrie
  393.             do j=1,ntot
  394. *     Résoudre L(1:1j,1:j) v(1:j) = A(1:j,j)
  395.                do k=1,j
  396.                   vvec(k)=amat(k,j)
  397.                enddo
  398.                do k=1,j-1
  399.                   do i=k+1,j
  400.                      vvec(i)=vvec(i)-vvec(k)*amat(i,k)
  401.                   enddo
  402.                enddo
  403. *     Calculer M(j,1:j-1) et stocker dans A (1:j-1,j)
  404.                do i=1,j-1
  405.                   amat(i,j)=vvec(i)/amat(i,i)
  406.                enddo
  407. *     Stocker d(j) dans A(j,j)
  408.                amat(j,j)=vvec(j)
  409.                if (j.le.nligp) then
  410.                   if (abs(amat(j,j)).le.xpethg) then
  411.                      if (ldbg) then
  412.                         write(ioimp,*) 'Pivot petit detecte ilrig,j='
  413.      $                       ,ilrig,j
  414.                         write(ioimp,*) '| ',amat(j,j),' | < ',xpethg
  415.                         write(ioimp,*) 'Noeud ',ipt1.num(des1.noelep(j)
  416.      $                       ,ilrig), ' Inconnue ',des1.lisinc(j)
  417.                         segprt,amats
  418.                         write(ioimp,*) 'Diag hg :'
  419.                         write(ioimp,*) (amats(jj,jj),jj=1,nligp)
  420.                         write(ioimp,*) 'Diag bd :'
  421.                         write(ioimp,*) (amats(jj,jj),jj=nligp+1,ntot)
  422.                         write(ioimp,*) 'Pivots hg :'
  423.                         write(ioimp,*) (amat(jj,jj),jj=1,nligp)
  424.                         write(ioimp,*) 'Pivots bd :'
  425.                         write(ioimp,*) (amat(jj,jj),jj=nligp+1,ntot)
  426.                      endif
  427.                      amat(j,j)=xpethg
  428.                      npivhg=npivhg+1
  429.                   endif
  430.                else
  431.                   if (abs(amat(j,j)).le.-xpetbd) then
  432.                      if (ldbg) then
  433.                         write(ioimp,*) 'Pivot petit detecte ilrig,j='
  434.      $                       ,ilrig,j
  435.                         write(ioimp,*) '| ',amat(j,j),' | < -1* ',xpetbd
  436.                         write(ioimp,*) 'Noeud ',ipt1.num(des3.noelep(j
  437.      $                       -nligp),ilrig), ' Inconnue ',des3.lisinc(j
  438.      $                       -nligp)
  439.                         segprt,amats
  440.                         write(ioimp,*) 'Diag hg :'
  441.                         write(ioimp,*) (amats(jj,jj),jj=1,nligp)
  442.                         write(ioimp,*) 'Diag bd :'
  443.                         write(ioimp,*) (amats(jj,jj),jj=nligp+1,ntot)
  444.                         write(ioimp,*) 'Pivots hg :'
  445.                         write(ioimp,*) (amat(jj,jj),jj=1,nligp)
  446.                         write(ioimp,*) 'Pivots bd :'
  447.                         write(ioimp,*) (amat(jj,jj),jj=nligp+1,ntot)
  448.                      endif
  449.                      amat(j,j)=xpetbd
  450.                      npivbd=npivbd+1
  451.                   endif
  452.                endif
  453. *     Calculer L(j+1:n,j) et stocker dans A (j+1:n,j)
  454.                do k=1,j-1
  455.                   do i=j+1,ntot
  456.                      amat(i,j)=amat(i,j)-vvec(k)*amat(i,k)
  457.                   enddo
  458.                enddo
  459. *               write(ioimp,*) 'ilrig,j,a(j,j),xpetbd,xpethg=',ilrig,j
  460. *     $              ,amat(j,j),xpetbd,xpethg
  461.                do i=j+1,ntot
  462.                   amat(i,j)=amat(i,j)/amat(j,j)
  463.                enddo
  464.             enddo
  465. *     Vérification de la factorisation
  466.             if (lvdec) then
  467. *               segprt,amat
  468.                do k=1,ntot
  469.                   do i=1,ntot
  470.                      amatv(i,k)=0.D0
  471.                   enddo
  472.                enddo
  473.                do k=1,ntot
  474.                   do i=1,ntot
  475.                      do j=1,min(i,k)
  476.                         xajj=amat(j,j)
  477.                         if (k.eq.j) then
  478.                            xukj=1.d0
  479.                         else
  480.                            xukj=amat(j,k)
  481.                         endif
  482.                         if (i.eq.j) then
  483.                            xlij=1.d0
  484.                         else
  485.                            xlij=amat(i,j)
  486.                         endif
  487.                         amatv(i,k)=amatv(i,k)+xlij*xajj*xukj
  488.                      enddo
  489.                   enddo
  490.                enddo
  491. *               segprt,amatv
  492.                do k=1,ntot
  493.                   do i=1,ntot
  494.                      xda=abs(amatv(i,k)-amats(i,k))
  495.                      if (xda.gt.xtol) then
  496.                         iprint=iprint+1
  497.                         write(ioimp,*)
  498.      $                       'Erreur factorisation detectee ilrig='
  499.      $                       ,ilrig
  500.                         write(ioimp,*) 'aik=',amats(i,k),'aikv=',amatv(i
  501.      $                       ,k)
  502.                         write(ioimp,*) 'i=',i,' k=',k
  503.                         write(ioimp,*) 'xpethg,xpetbd,xtol=',xpethg
  504.      $                       ,xpetbd,xtol
  505.                      write(ioimp,*) 'Diag hg orig :'
  506.                      write(ioimp,*) (amats(jj,jj),jj=1,nligp)
  507.                      write(ioimp,*) 'Diag bd orig :'
  508.                      write(ioimp,*) (amats(jj,jj),jj=nligp+1,ntot)
  509.                      write(ioimp,*) 'Diag hg :'
  510.                      write(ioimp,*) (amatv(jj,jj),jj=1,nligp)
  511.                      write(ioimp,*) 'Diag bd :'
  512.                      write(ioimp,*) (amatv(jj,jj),jj=nligp+1,ntot)
  513.                      segprt,amats
  514.                      segprt,amatv
  515.                      ierr=1
  516.                      return
  517. *                        write(ioimp,*) 'xda = ',xda,' > ',xtol
  518. *      write(ioimp,*) 'Noeud dual i ',ipt1.num(des3.noeled(i)
  519. *     $                    ,ilrig), ' Inconnue ',des3.lisdua(i)
  520. *      write(ioimp,*) 'Noeud primal k',ipt1.num(des3.noelep(k)
  521. *     $                    ,ilrig), ' Inconnue ',des3.lisinc(k)
  522.                         if (iprint.gt.nprint) then
  523.                            call erreur(185)
  524.                            return
  525.                         endif
  526.                      endif
  527.                   enddo
  528.                enddo
  529.             endif
  530. * Produit LDMt avec le bloc en bas à droite et stockage dans RE
  531. * Amélioration possible : prendre en compte la symétrie
  532.             do k=1,nligd
  533.                k2 = nligp+k
  534.                do i=1,nligd
  535.                   i2 = nligp+i
  536.                   do j=1,min(i,k)
  537.                      j2 = nligp+j
  538.                      xajj=amat(j2,j2)
  539.                      if (k.eq.j) then
  540.                         xukj=1.d0
  541.                      else
  542.                         xukj=amat(j2,k2)
  543.                      endif
  544.                      if (i.eq.j) then
  545.                         xlij=1.d0
  546.                      else
  547.                         xlij=amat(i2,j2)
  548.                      endif
  549.                      xmatr3.re(i,k,ilrig)=xmatr3.re(i,k,ilrig)
  550.      $                    +xlij*xajj*xukj
  551.                   enddo
  552.                enddo
  553.             enddo
  554.          ENDDO
  555.          segsup vvec
  556.          if (lvdec) then
  557.             segsup amats
  558.             segsup amatv
  559.          endif
  560.          segsup amat
  561.          IRIG2.IRIGEL(4,irig)=xmatr3
  562.          IRIG2.IRIGEL(5,irig)=IRIGC2.IRIGEL(5,irig)
  563.       enddo
  564. *     Impression pivots nuls
  565.       npivnu=npivhg+npivbd
  566.       if (npivnu.gt.0.and.iimpi.ne.0) then
  567. *      if (npivnu.gt.0) then
  568.          xpivel=float(npivnu)/float(nelpnu)
  569.          xpilhg=float(npivhg)/float(nelpnu)
  570.          xpilbd=float(npivbd)/float(nelpnu)
  571.          write(ioimp,*) npivnu,' pivots nuls detectes / ',nelpnu,
  572.      $        ' elements = ',xpivel,' pivots nuls/element'
  573.          write(ioimp,*) '  dont : ', xpilhg
  574.      $        ,' pivots nuls bloc(1,1)/element'
  575.          write(ioimp,*) '         ', xpilbd
  576.      $        ,' pivots nuls bloc(2,2)/element'
  577.       endif
  578. *      segprt,irig2
  579. *
  580. * Menage
  581. *
  582.       DO irig=1,nbsou
  583. *         descr=irigc2.irigel(3,irig)
  584. *         segsup descr
  585.          xmatri=irigc2.irigel(4,irig)
  586.          segsup xmatri
  587.          IF (IRIGB.NE.IRIGC) THEN
  588. *            descr=irigb2.irigel(3,irig)
  589. *            segsup descr
  590.             xmatri=irigb2.irigel(4,irig)
  591.             segsup xmatri
  592.          ENDIF
  593. *         descr=irigd2.irigel(3,irig)
  594. *        segsup descr
  595.          xmatri=irigd2.irigel(4,irig)
  596.          segsup xmatri
  597.          IF (IRIGE.NE.0) THEN
  598. *            descr=irige2.irigel(3,irig)
  599. *            segsup descr
  600.             xmatri=irige2.irigel(4,irig)
  601.             segsup xmatri
  602.          ENDIF
  603.       ENDDO
  604.       segsup irigc2
  605.       IF (IRIGB.NE.IRIGC) segsup irigb2
  606.       segsup irigd2
  607.       IF (IRIGE.NE.0) segsup irige2
  608. *
  609.       icpr=jcpr
  610.       inode=jnode
  611.       jelnum=kelnum
  612.       izone=jzone
  613.       segsup izone
  614.       segsup jelnum
  615.       segsup inode
  616.       segsup icpr
  617.       SEGSUP INFMEL
  618. *
  619.       segsup meleme
  620.       RETURN
  621.       END
  622.  
  623.  
  624.  
  625.  

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