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/05/15    21:15:04     12266          
  2.       SUBROUTINE CMCT4(IRIGC,IRIGD,IRIGB,IRIGE,IMEL,IRIG2)
  3. *_______________________________________________________________________
  4. c
  5. c      opérateur 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 IRIGE.MRIGID,IRIGE2.MRIGID
  31.       POINTEUR IRIG2.MRIGID
  32.       POINTEUR XMATRC.XMATRI,XMATRD.XMATRI,XMATRB.XMATRI
  33.       POINTEUR XMATRE.XMATRI
  34.       POINTEUR XMAT.XMATRI
  35.       POINTEUR XMATV.XMATRI
  36.       SEGMENT vvec(ntot)
  37.       SEGMENT amat(ntot,ntot)
  38.       SEGMENT amats(ntot,ntot)
  39.       SEGMENT amatv(ntot,ntot)
  40. -INC SMELEME
  41. -INC CCHAMP
  42. -INC SMCOORD
  43. -INC SMCHPOI
  44. *
  45.       CHARACTER*(LOCHPO) NOMINC,NOMDUA
  46. *
  47.       segment icpr(nbpts)
  48.       segment inode(ino)
  49.       segment jelnum(imaxel,ino)
  50.       segment izone(imaxel,ino)
  51.       segment INFMEL
  52.       integer jmel,jcpr,jnode,kelnum,jzone
  53.       endsegment
  54.  
  55.  
  56.       logical lvdec
  57. *
  58.       lvdec=.false.
  59.       nprint=40
  60.       iprint=0
  61.       npivnu=0
  62.       nelpnu=0
  63.       IF (IMEL.NE.0) THEN
  64.          IPT1=IMEL
  65.       ELSE
  66.          CALL MELRIG(IRIGD,IPT1)
  67.       ENDIF
  68.       CALL REGMAI(IPT1,MELEME)
  69.       IF (IERR.NE.0) RETURN
  70. *
  71. *  On suppose que l'objet GEOMETRIE de IRIG2 est le MELEME donne en
  72. *     entree
  73. *
  74. * Etape I : reduire les matrices d'entree sur ce MELEME
  75. *           puis construire les descripteurs globaux associes
  76. *
  77.       CALL ACTOBJ('MAILLAGE',MELEME,1)
  78.       nbsou=lisous(/1)
  79. C   Cas d'un MELEME vide => RIGIDITE vide
  80.       IF (nbsou.eq.0) THEN
  81.          NRIGEL = 0
  82.          SEGINI,IRIG2
  83.          IRIG2.MTYMAT = 'KPRE'
  84.          IRIG2.IFORIG = IFOUR
  85.          SEGDES,IRIG2
  86.          RETURN
  87.       ENDIF
  88.       SEGINI INFMEL
  89.       INFMEL.JMEL=MELEME
  90. *
  91. * Réduction de C sur le maillage MELEME
  92. *
  93. *      write(ioimp,*) ' IRIGC '
  94. *      call ecrobj('RIGIDITE',irigc)
  95. *      call ecrcha('RESU')
  96. *      call prlist
  97.       IRIGC2=0
  98.       CALL CMCT5(IRIGC,INFMEL,IRIGC2)
  99.       IF (IERR.NE.0) RETURN
  100. *      write(ioimp,*) ' IRIGC2 '
  101. *      call ecrobj('RIGIDITE',irigc2)
  102. *      call ecrcha('RESU')
  103. *      call prlist
  104. *      return
  105. *
  106. * Réduction eventuelle de B sur le maillage MELEME et sur le DESCR de C2
  107. *
  108.       IF (IRIGB.NE.IRIGC) THEN
  109. * Amélioration possible : transposition auto de B si necessaire
  110. *         write(ioimp,*) ' IRIGB .NE. IRIGC'
  111.          segini,irigb2=irigc2
  112. * On enleve les xmatri
  113.          nrig=irigb2.coerig(/1)
  114.          do irig=1,nrig
  115.             irigb2.irigel(4,irig)=0
  116.          enddo
  117.          CALL CMCT5(IRIGB,INFMEL,IRIGB2)
  118.          IF (IERR.NE.0) RETURN
  119.       ELSE
  120. *         write(ioimp,*) ' IRIGB = IRIGC'
  121.       ENDIF
  122. *
  123. * Reduction de D sur le maillage MELEME et sur le DESCR primal de C2
  124. *
  125.       NRIGEL=nbsou
  126.       SEGINI,IRIGD2
  127.       SEGACT,IRIGC2
  128.       DO irig=1,NRIGEL
  129.          des1=irigc2.irigel(3,irig)
  130.          nligrp=des1.noelep(/1)
  131.          nligrd=nligrp
  132.          segini des3
  133.          do iligrp=1,nligrp
  134.             nno=des1.noelep(iligrp)
  135.             des3.noelep(iligrp)=nno
  136.             des3.noeled(iligrp)=nno
  137.             nominc=des1.lisinc(iligrp)
  138.             des3.lisinc(iligrp)=nominc
  139.             CALL PLACE(NOMDD,LNOMDD,IPLA,NOMINC)
  140.             IF (IPLA.EQ.0) THEN
  141.                des3.lisdua(iligrp)=nominc
  142.             ELSE
  143.                des3.lisdua(iligrp)=nomdu(ipla)
  144.             ENDIF
  145.          enddo
  146.          irigd2.irigel(3,irig)=des3
  147.       ENDDO
  148. *      write(ioimp,*) ' IRIGD '
  149.       CALL CMCT5(IRIGD,INFMEL,IRIGD2)
  150.       IF (IERR.NE.0) RETURN
  151. *      call ecrobj('RIGIDITE',irigd2)
  152. *      call ecrcha('RESU')
  153. *      call prlist
  154. *      return
  155. *
  156. * Reduction (eventuelle) de E sur le maillage MELEME et sur le DESCR dual de C2
  157. *
  158.       IF (IRIGE.NE.0) THEN
  159.          NRIGEL=nbsou
  160.          SEGINI,IRIGE2
  161.          SEGACT,IRIGC2
  162.          DO irig=1,NRIGEL
  163.             des1=irigc2.irigel(3,irig)
  164.             nligrd=des1.noeled(/1)
  165.             nligrp=nligrd
  166.             segini des3
  167.             do iligrd=1,nligrd
  168.                nno=des1.noeled(iligrd)
  169.                des3.noeled(iligrd)=nno
  170.                des3.noelep(iligrd)=nno
  171.                nomdua=des1.lisdua(iligrd)
  172.                des3.lisdua(iligrd)=nomdua
  173.                CALL PLACE(NOMDU,LNOMDU,IPLA,NOMDUA)
  174.                IF (IPLA.EQ.0) THEN
  175.                   des3.lisinc(iligrd)=nomdua
  176.                ELSE
  177.                   des3.lisinc(iligrd)=nomdd(ipla)
  178.                ENDIF
  179.             enddo
  180.             irige2.irigel(3,irig)=des3
  181.          ENDDO
  182. *         write(ioimp,*) ' IRIGE '
  183.          CALL CMCT5(IRIGE,INFMEL,IRIGE2)
  184.          IF (IERR.NE.0) RETURN
  185.       ELSE
  186.          IRIGE2=0
  187.       ENDIF
  188. *
  189. * On est pret maintenant à calculer E - C D-1 Bt
  190. *     On envisage trois algorithmes :
  191. *     1) Force Brute : inversion de D puis produit
  192. *     2) Factorisation de la matrice  (D Bt) = (Lu   0 ) ( Uu  Upu)
  193. *                                     (C 0 )   (Lup  Lp) ( 0   Up )
  194. *      Si je comprends bien, Lp * Up = C D-1 Bt Oui !
  195. *     3) Brute brute : d'après la formulation d'inversion, on inverse
  196. *     tout, puis on inverse le bloc pp de l'inverse. Ca nous redonne le
  197. *     complément de Schur
  198. *
  199. *
  200. *     Algorithme factorisation LDMt: on pourrait pivoter symétriquement ?
  201. *     Repris de Golub et Van Loan !
  202. *
  203.       NRIGEL=nbsou
  204.       SEGINI IRIG2
  205.       SEGACT IRIGC2
  206.       DO irig=1,NRIGEL
  207.          IRIG2.MTYMAT = 'KPRE'
  208.          IRIG2.IFORIG = IFOUR
  209.          IRIG2.COERIG(irig)=1.D0
  210.          ipt1=LISOUS(irig)
  211.          IRIG2.IRIGEL(1,irig)=ipt1
  212.          des1=irigc2.irigel(3,irig)
  213.          nligd=des1.noeled(/1)
  214.          nligp=des1.noelep(/1)
  215.          nlrig=ipt1.num(/2)
  216.          nelpnu=nelpnu+nlrig
  217. * Matrice de stockage totale
  218.          nligrd=nligd
  219.          nligrp=nligd
  220.          segini des3
  221.          do iligrd=1,nligrd
  222.             nno=des1.noeled(iligrd)
  223.             des3.noeled(iligrd)=nno
  224.             des3.noelep(iligrd)=nno
  225.             nomdua=des1.lisdua(iligrd)
  226.             des3.lisdua(iligrd)=nomdua
  227.             CALL PLACE(NOMDU,LNOMDU,IPLA,NOMDUA)
  228.             IF (IPLA.EQ.0) THEN
  229.                des3.lisinc(iligrd)=nomdua
  230.             ELSE
  231.                des3.lisinc(iligrd)=nomdd(ipla)
  232.             ENDIF
  233.          enddo
  234.          IRIG2.IRIGEL(3,irig)=des3
  235.          isym=2
  236.          IF (IRIGB.EQ.IRIGC) THEN
  237.             IF (irigd2.irigel(7,irig).eq.0) THEN
  238.                if (irige2.ne.0) then
  239.                   isym=irige2.irigel(7,irig)
  240.                else
  241.                   isym=0
  242.                endif
  243.             endif
  244.          endif
  245.          IRIG2.IRIGEL(7,irig)=isym
  246. *
  247. * Remplissage des matrices élémentaires
  248. *
  249.          nelrig=nlrig
  250.          segini xmatr3
  251.          xmatrd=irigd2.irigel(4,irig)
  252. *         segprt,xmatrd
  253.          xmatrc=irigc2.irigel(4,irig)
  254.          if (irigb.ne.irigc) then
  255.             xmatrb=irigb2.irigel(4,irig)
  256.          else
  257.             xmatrb=xmatrc
  258.          endif
  259.          if (irige.ne.irige) then
  260.             xmatre=irige2.irigel(4,irig)
  261.          else
  262.             xmatre=0
  263.          endif
  264.          ntot=nligp+nligd
  265.          segini amat
  266.          if (lvdec) then
  267.             segini amats
  268.             segini amatv
  269.          endif
  270.          segini vvec
  271.          do ilrig=1,nlrig
  272. *     Copies
  273.             do iligp=1,nligp
  274.                do iligd=1,nligp
  275. *                  segprt,xmatrd
  276. *                  write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd
  277.                   amat(iligd,iligp)=xmatrd.re(iligd,iligp,ilrig)
  278. *                  write(ioimp,*) 'amat(iligd,iligp)=',amat(iligp,iligd)
  279.                enddo
  280.             enddo
  281.             do iligp=1,nligp
  282.                do iligd=1,nligd
  283.                   amat(iligd+nligp,iligp)=
  284.      $                 xmatrc.re(iligd,iligp,ilrig)
  285.                enddo
  286.             enddo
  287.             do iligp=1,nligp
  288.                do iligd=1,nligd
  289.                   amat(iligp,iligd+nligp)=
  290.      $                 xmatrb.re(iligd,iligp,ilrig)
  291.                enddo
  292.             enddo
  293.             if (xmatre.ne.0) then
  294.                do iligp=1,nligd
  295.                   do iligd=1,nligd
  296.                      amat(iligp+nligp,iligd+nligp)=
  297.      $                    xmatre.re(iligd,iligp,ilrig)
  298.                   enddo
  299.                enddo
  300.             else
  301.                do iligp=1,nligd
  302.                   do iligd=1,nligd
  303.                      amat(iligp+nligp,iligd+nligp)=0.D0
  304.                   enddo
  305.                enddo
  306.             endif
  307. * Sauvegarde de A dans As
  308.             if (lvdec) then
  309.                do j=1,ntot
  310.                   do i=1,ntot
  311.                      amats(i,j)=amat(i,j)
  312.                   enddo
  313.                enddo
  314. *               segprt,amats
  315.             endif
  316. *  Echelle par la partie haut gauche
  317.             xscahg=0.d0
  318.             do j=1,nligp
  319.                xscahg=max(xscahg,amat(j,j))
  320.             enddo
  321. * xpethg positi
  322.             xpethg=xscahg*xzprec*1000.d0
  323.             if (xpethg.EQ.0.d0) then
  324.                segprt,des3
  325.                write(ioimp,*) 'irig,nrigel=',irig,nrigel
  326.                write(ioimp,*) 'ilrig,nlrig=',ilrig,nlrig
  327.                write(ioimp,*) 'nligp,nligd=',nligp,nligd
  328.                segprt,amat
  329.                call erreur(185)
  330.                return
  331.             endif
  332. *
  333.             xscabd=0.d0
  334.             do j=1,nligp
  335.                xest=0.d0
  336.                xpiv=max(xpethg,amat(j,j))
  337.                do i=1,nligd
  338.                   xest=xest-amat(i+nligp,j)*amat(j,i+nligp)/xpiv
  339.                enddo
  340.                xscabd=min(xscabd,xest)
  341.             enddo
  342. * xpetbd negatif
  343.             xpetbd=xscabd*xzprec*1000.d0
  344.             if (lvdec) then
  345.                xtol=max(xpethg,abs(xpetbd))*100.d0
  346.             endif
  347. *
  348. *     Factorisation de amat (LDMt ou LDLt d'apres Golub et Van Loan 3eme
  349. *     edition pp. 137-139)
  350. * Amélioration possible : prendre en compte la symétrie
  351.             do j=1,ntot
  352. *     Résoudre L(1:1j,1:j) v(1:j) = A(1:j,j)
  353.                do k=1,j
  354.                   vvec(k)=amat(k,j)
  355.                enddo
  356.                do k=1,j-1
  357.                   do i=k+1,j
  358.                      vvec(i)=vvec(i)-vvec(k)*amat(i,k)
  359.                   enddo
  360.                enddo
  361. *     Calculer M(j,1:j-1) et stocker dans A (1:j-1,j)
  362.                do i=1,j-1
  363.                   amat(i,j)=vvec(i)/amat(i,i)
  364.                enddo
  365. *     Stocker d(j) dans A(j,j)
  366.                amat(j,j)=vvec(j)
  367.                if (j.le.nligp) then
  368.                   if (abs(amat(j,j)).lt.xpethg) then
  369. c$$$                     write(ioimp,*) 'Pivot petit detecte | ',amat(j
  370. c$$$     $                    ,j),' | < ',xpethg
  371. c$$$                     write(ioimp,*) 'Noeud ',ipt1.num(des1.noelep(j)
  372. c$$$     $                    ,ilrig), ' Inconnue ',des1.lisinc(j)
  373.                      amat(j,j)=xpethg
  374.                      npivnu=npivnu+1
  375.                   endif
  376.                else
  377.                   if (abs(amat(j,j)).lt.-xpetbd) then
  378. c$$$                     write(ioimp,*) 'Pivot petit detecte | ',amat(j
  379. c$$$     $                    ,j),' | < -1* ',xpetbd
  380. c$$$                     write(ioimp,*) 'Noeud ',ipt1.num(des1.noeled(j
  381. c$$$     $                    -nligp),ilrig), ' Inconnue ',des1.lisdua(j
  382. c$$$     $                    -nligp)
  383.                      amat(j,j)=xpetbd
  384.                      npivnu=npivnu+1
  385.                   endif
  386.                endif
  387. *     Calculer L(j+1:n,j) et stocker dans A (j+1:n,j)
  388.                do k=1,j-1
  389.                   do i=j+1,ntot
  390.                      amat(i,j)=amat(i,j)-vvec(k)*amat(i,k)
  391.                   enddo
  392.                enddo
  393.                do i=j+1,ntot
  394.                   amat(i,j)=amat(i,j)/amat(j,j)
  395.                enddo
  396.             enddo
  397. *     Vérification de la factorisation
  398.             if (lvdec) then
  399. *               segprt,amat
  400.                do k=1,ntot
  401.                   do i=1,ntot
  402.                      amatv(i,k)=0.D0
  403.                   enddo
  404.                enddo
  405.                do k=1,ntot
  406.                   do i=1,ntot
  407.                      do j=1,min(i,k)
  408.                         xajj=amat(j,j)
  409.                         if (k.eq.j) then
  410.                            xukj=1.d0
  411.                         else
  412.                            xukj=amat(j,k)
  413.                         endif
  414.                         if (i.eq.j) then
  415.                            xlij=1.d0
  416.                         else
  417.                            xlij=amat(i,j)
  418.                         endif
  419.                         amatv(i,k)=amatv(i,k)+xlij*xajj*xukj
  420.                      enddo
  421.                   enddo
  422.                enddo
  423. *               segprt,amatv
  424.                do k=1,ntot
  425.                   do i=1,ntot
  426.                      xda=abs(amatv(i,k)-amats(i,k))
  427.                      if (xda.gt.xtol) then
  428.                         iprint=iprint+1
  429.                         write(ioimp,*) 'Erreur factorisation detectee '
  430.      $                       ,'aik=',amats(i,k),'aikv=',amatv(i,k)
  431.                         write(ioimp,*) 'i=',i,' k=',k
  432.                         write(ioimp,*) 'xpethg,xpetbd,xtol=',xpethg
  433.      $                       ,xpetbd,xtol
  434.                         write(ioimp,*) 'xda = ',xda,' > ',xtol
  435. *     write(ioimp,*) 'Noeud i ',ipt1.num(des1.noelep(j)
  436. *     $                    ,ilrig), ' Inconnue ',des1.lisinc(j)
  437. *     write(ioimp,*) 'Noeud ',ipt1.num(des1.noelep(j)
  438. *     $                    ,ilrig), ' Inconnue ',des1.lisinc(j)
  439.                         if (iprint.gt.nprint) then
  440.                            call erreur(185)
  441.                            return
  442.                         endif
  443.                      endif
  444.                   enddo
  445.                enddo
  446.             endif
  447. * Produit LDMt avec le bloc en bas à droite et stockage dans RE
  448. * Amélioration possible : prendre en compte la symétrie
  449.             do k=1,nligd
  450.                k2 = nligp+k
  451.                do i=1,nligd
  452.                   i2 = nligp+i
  453.                   do j=1,min(i,k)
  454.                      j2 = nligp+j
  455.                      xajj=amat(j2,j2)
  456.                      if (k.eq.j) then
  457.                         xukj=1.d0
  458.                      else
  459.                         xukj=amat(j2,k2)
  460.                      endif
  461.                      if (i.eq.j) then
  462.                         xlij=1.d0
  463.                      else
  464.                         xlij=amat(i2,j2)
  465.                      endif
  466.                      xmatr3.re(i,k,ilrig)=xmatr3.re(i,k,ilrig)
  467.      $                    +xlij*xajj*xukj
  468.                   enddo
  469.                enddo
  470.             enddo
  471.          ENDDO
  472.          segsup vvec
  473.          if (lvdec) then
  474.             segsup amats
  475.             segsup amatv
  476.          endif
  477.          segsup amat
  478.          IRIG2.IRIGEL(4,irig)=xmatr3
  479.          IRIG2.IRIGEL(5,irig)=IRIGC2.IRIGEL(5,irig)
  480.       enddo
  481. * Impression pivots nuls
  482.       if (npivnu.gt.0.and.iimpi.ne.0) then
  483.          xpivel=float(npivnu)/float(nelpnu)
  484.          write(ioimp,*) npivnu,' pivots nuls detectes / ',nelpnu,
  485.      $        ' elements = ',xpivel,' pivots nuls/element'
  486.       endif
  487. *
  488. * Menage
  489. *
  490.       DO irig=1,nbsou
  491.          descr=irigc2.irigel(3,irig)
  492.          segsup descr
  493.          xmatri=irigc2.irigel(4,irig)
  494.          segsup xmatri
  495.          IF (IRIGB.NE.IRIGC) THEN
  496.             descr=irigb2.irigel(3,irig)
  497.             segsup descr
  498.             xmatri=irigb2.irigel(4,irig)
  499.             segsup xmatri
  500.          ENDIF
  501.          descr=irigd2.irigel(3,irig)
  502.          segsup descr
  503.          xmatri=irigd2.irigel(4,irig)
  504.          segsup xmatri
  505.          IF (IRIGE.NE.0) THEN
  506.             descr=irige2.irigel(3,irig)
  507.             segsup descr
  508.             xmatri=irige2.irigel(4,irig)
  509.             segsup xmatri
  510.          ENDIF
  511.       ENDDO
  512.       segsup irigc2
  513.       IF (IRIGB.NE.IRIGC) segsup irigb2
  514.       segsup irigd2
  515.       IF (IRIGE.NE.0) segsup irige2
  516. *
  517.       icpr=jcpr
  518.       inode=jnode
  519.       jelnum=kelnum
  520.       izone=jzone
  521.       segsup izone
  522.       segsup jelnum
  523.       segsup inode
  524.       segsup icpr
  525.       SEGSUP INFMEL
  526. *
  527.       segsup meleme
  528.       RETURN
  529.       END
  530.  
  531.  

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