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

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