Télécharger la plaquette Cast3M


fr en

Télécharger cmct5.eso

Retour à la liste

Numérotation des lignes : 

cmct5

  1. C CMCT5     SOURCE    GOUNAND   25/06/11    21:15:05     12278          
  2.       SUBROUTINE CMCT5(IRI1,INFMEL,IRI3)
  3. *_______________________________________________________________________
  4. c
  5. c      opérateur KPRE (cmct option ELEM)
  6. c
  7. c     Réduction de la rigidité IRI1 sur le maillage IMEL
  8. c     ICPR,INODE,JELNUM,IZONE sont des objets de
  9. c     preconditionnement
  10. C     associes a IMEL
  11. C
  12. C
  13. C     En entree/sortie : IRI3 est definie sur IMEL et son nombre d'elements
  14. C                        est le meme que IMEL
  15. C                 Si IRI3=0 en entrée, on cree son descripteur DES3
  16. C                    par assemblage des descripteurs DES1 de IRI1
  17. C                    si necessaire
  18. C                 SI IRI3 different de 0 en entree, on utilise les descripteurs
  19. C                    DES3 deja presents
  20. C
  21. C Gérer le cas où la matrice est deja réduite comme il faut ???
  22. C
  23. c
  24. *_______________________________________________________________________
  25.  
  26.       IMPLICIT INTEGER(I-N)
  27.       IMPLICIT REAL*8(A-H,O-Z)
  28. *
  29.  
  30. -INC PPARAM
  31. -INC CCOPTIO
  32. -INC CCREEL
  33. -INC SMRIGID
  34. -INC SMELEME
  35. -INC CCHAMP
  36. -INC SMCOORD
  37. -INC SMCHPOI
  38.       segment icpr(nbpts)
  39.       segment inode(ino)
  40.       segment jelnum(imaxel,ino)
  41.       segment izone(imaxel,ino)
  42.       segment INFMEL
  43.       integer jmel,jcpr,jnode,kelnum,jzone
  44.       endsegment
  45. *
  46.       segment ilmasq(nele1)
  47.       segment iddlp(nligrp)
  48.       segment iddld(nligrd)
  49.       segment INFRIG
  50.          INTEGER jzon(NRIGEL)
  51.          INTEGER jlmasq(NRIGEL)
  52.          INTEGER jddlp(NRIGEL)
  53.          INTEGER jddld(NRIGEL)
  54.       endsegment
  55.       logical ldescr,ldbg
  56.       ldbg=.false.
  57. *
  58.       MELEME=jmel
  59. *      segprt,meleme
  60.       nbsou=lisous(/1)
  61.       if (jcpr.eq.0) then
  62. *
  63. *   creation d'une numerotation locale
  64. *
  65.          segini icpr
  66.          ino=0
  67.          do  1 i =1, nbsou
  68.             ipt1=lisous(i)
  69. *            write(ioimp,*) 'cmct5 i=',i,' ipt1=',ipt1
  70. *            segprt,ipt1
  71.             do 2 j=1,ipt1.num(/2)
  72.                do 22 k=1,ipt1.num(/1)
  73.                   ia= ipt1.num(k,j)
  74.                   if(icpr(ia).eq.0) then
  75.                      ino=ino+1
  76.                      icpr(ia)=ino
  77.                   endif
  78.  22            continue
  79.  2          continue
  80.  1       continue
  81. *         write(ioimp,*) 'ino=',ino
  82. *         segprt,icpr
  83. *
  84. * on compte combien d'elements touche un noeud
  85. *
  86.          segini inode
  87.          do 3 i=1,nbsou
  88.             ipt1=lisous(i)
  89.             do 4 j=1,ipt1.num(/2)
  90.                do 42 k=1,ipt1.num(/1)
  91.                   ia=ipt1.num(k,j)
  92.                   ib= icpr(ia)
  93. *               write(ioimp,*) 'ia,ib=',ia,ib
  94.                   inode(ib)=inode(ib)+1
  95.  42            continue
  96.  4          continue
  97.  3       continue
  98.          imaxel=0
  99.          do i=1,ino
  100.             imaxel=max(imaxel,inode(i))
  101.             inode(i)=0
  102.          enddo
  103.          segini jelnum,izone
  104. *
  105. * jelnum(i,k) dira le Ieme element qui touche le noeud k
  106. * izone (i,k) dira dans quel lisous le ieme element se trouve
  107. * attention le noeud K est le plus petit de l'élément
  108. *
  109.          do 5 i=1,nbsou
  110.             ipt1=lisous(i)
  111.             do 6 j=1,ipt1.num(/2)
  112.                ipeti=igrand
  113.                do  k=1,ipt1.num(/1)
  114.                   ipeti=min(ipt1.num(k,j),ipeti)
  115.                enddo
  116. *           write(ioimp,*) 'ipeti=',ipeti
  117.                ib= icpr(ipeti)
  118.                inode(ib)=inode(ib)+1
  119.                ic=inode(ib)
  120.                jelnum(ic,ib)=j
  121.                izone(ic,ib)=i
  122.  6          continue
  123.  5       continue
  124. *         segprt,izone
  125. *         segprt,jelnum
  126. * Sauvegarde des objets créés dans INFMEL
  127.          jcpr=icpr
  128.          jnode=inode
  129.          kelnum=jelnum
  130.          jzone=izone
  131.       endif
  132.       icpr=jcpr
  133.       inode=jnode
  134.       jelnum=kelnum
  135.       izone=jzone
  136. *
  137. * Matrice réduite (sortie)
  138. *
  139.       RI1=IRI1
  140.       SEGACT MELEME
  141.       IF (IRI3.EQ.0) THEN
  142.          ldescr=.true.
  143.          NRIG3=LISOUS(/1)
  144.          NRIGEL=NRIG3
  145.          SEGINI RI3
  146.          IRI3=RI3
  147.          RI3.MTYMAT='CMCT5'
  148.          RI3.IFORIG=IFOUR
  149.          RI3.MCRCNF=MCOORD
  150.          NLIGRP=0
  151.          NLIGRD=0
  152.          DO irig3=1,nrig3
  153.             SEGINI,DES3
  154.             RI3.IRIGEL(3,irig3)=DES3
  155.          ENDDO
  156.       ELSE
  157.          ldescr=.false.
  158.          RI3=IRI3
  159.          SEGACT RI3*MOD
  160.          NRIG3=RI3.COERIG(/1)
  161.          IF (NRIG3.NE.LISOUS(/1)) THEN
  162.             write(ioimp,*) 'Erreur grave 1 dans cmct5'
  163.             call erreur(5)
  164.             return
  165.          endif
  166.          DO irig3=1,nrig3
  167.             DES3=RI3.IRIGEL(3,irig3)
  168.             SEGACT,DES3
  169.          ENDDO
  170.       ENDIF
  171.       RI3.MTYMAT='CMCT5'
  172.       RI3.IFORIG=IFOUR
  173.       RI3.MCRCNF=MCOORD
  174.       DO irig3=1,nrig3
  175.          RI3.COERIG(irig3)=1.D0
  176.          ipt3=LISOUS(irig3)
  177.          RI3.IRIGEL(1,irig3)=ipt3
  178.       ENDDO
  179. *
  180. * Traitement de la matrice d'entree
  181. *
  182.       SEGACT RI1
  183.       NRIG1=RI1.COERIG(/1)
  184.       NRIGEL=NRIG1
  185.       SEGINI INFRIG
  186. *
  187.       DO IRIG1=1,nrig1
  188.          IPT1=RI1.IRIGEL(1,IRIG1)
  189. * Y a-t-il des elements a considerer
  190.          SEGACT IPT1
  191. *         write(ioimp,*) 'IRIG1=',irig1, ' ipt1=',ipt1
  192. *         segprt,ipt1
  193.          DES1=RI1.IRIGEL(3,IRIG1)
  194.          SEGACT DES1
  195.          nele1=ipt1.num(/2)
  196.          segini ilmasq
  197.          nnod1=ipt1.num(/1)
  198.          izon =0
  199.          ielzo=0
  200.          DO 11 iele1=1,nele1
  201.             ipeti=IGRAND
  202.             DO inod1=1,nnod1
  203.                ipeti=min(ipeti,ipt1.num(inod1,iele1))
  204.             enddo
  205. * on regarde s'il existe un element de ipt1 ayant ce noeud en plus
  206. * petite position sinon on passe a l'élément suivant
  207.             ib=icpr(ipeti)
  208. *            write(ioimp,*) 'iele1=',iele1,' ipeti=',ipeti,' ib=',ib
  209.             if(ib.eq.0) go to 11
  210.             if(inode(ib).eq.0) go to 11
  211.             do 13 mm=1,inode(ib)
  212.                kzon=izone(mm,ib)
  213.                ipt3 = lisous(kzon)
  214.                if (ipt3.eq.ipt1) then
  215.                   segsup ilmasq
  216.                   ilmasq=-1
  217.                   izon=kzon
  218.                   ielzo=ipt3.num(/2)
  219.                   goto 14
  220.                endif
  221.                if (ipt3.num(/1).ne.ipt1.num(/1)) go to 13
  222.                iel=jelnum(mm,ib)
  223.                do  in=1,ipt1.num(/1)
  224.                   if( ipt1.num(in,iele1).ne.ipt3.num(in,iel))go to 13
  225.                enddo
  226. * on a trouve un element a conserver
  227.                ielzo=ielzo+1
  228.                ilmasq(iele1)=iel
  229.                if (izon.eq.0) then
  230.                   izon=kzon
  231.                else
  232.                   if (izon.ne.kzon) then
  233.                      write(ioimp,*)
  234.      $                    'irig1,ipt1,nele1,nnod1,iele1,ipeti=', irig1
  235.      $                    ,ipt1,nele1,nnod1,iele1,ipeti
  236.                      write(ioimp,*) 'izon,kzon=',izon,kzon
  237.                      write(ioimp,*) 'Erreur grave 2 dans cmct5'
  238.                      call erreur(5)
  239.                      return
  240.                   endif
  241.                endif
  242.                goto 11
  243.  13         continue
  244.  11      continue
  245.  14      continue
  246.          jzon(irig1)=izon
  247.          if (ldbg) write(ioimp,*) 'irig1=',irig1,' ielzo=',ielzo
  248.      $        ,'/ nele1,nnod1=',nele1,nnod1
  249.          if (izon.eq.0.or.ielzo.eq.0) then
  250.             segsup ilmasq
  251.             ilmasq=0
  252.          endif
  253.          jlmasq(irig1)=ilmasq
  254.          if (ilmasq.ne.0) then
  255. *            if (ilmasq.gt.0) then
  256. *               segprt,ilmasq
  257. *            else
  258. *               write(ioimp,*) 'ilmasq=',ilmasq
  259. *            endif
  260. *     tableau de correspondance des DESCR
  261.             DES3=RI3.IRIGEL(3,IZON)
  262.             DES1=RI1.IRIGEL(3,IRIG1)
  263.             SEGACT DES1
  264. *       pour la primale
  265.             NLIGP1=DES1.NOELEP(/1)
  266.             NLIGRP=NLIGP1
  267.             SEGINI IDDLP
  268.             DO ILIGP1=1,NLIGP1
  269.                nno1=des1.noelep(iligp1)
  270.                ilig3=0
  271.                nligp3=des3.noelep(/1)
  272.                do 113 iligp3=1,nligp3
  273.                   if (des3.noelep(iligp3).ne.nno1) goto 113
  274.                   if (des3.lisinc(iligp3).ne.des1.lisinc(iligp1))
  275.      $                 goto 113
  276.                   ilig3=iligp3
  277.                   goto 114
  278.  113           continue
  279.  114           continue
  280.                IF (ilig3.EQ.0) THEN
  281.                   if (ldescr) then
  282.                      NLIGP3=NLIGP3+1
  283.                      NLIGRP=NLIGP3
  284.                      NLIGRD=des3.noeled(/1)
  285.                      SEGADJ DES3
  286.                      des3.noelep(nligp3)=des1.noelep(iligp1)
  287.                      des3.lisinc(nligp3)=des1.lisinc(iligp1)
  288.                      IDDLP(ILIGP1)=nligp3
  289.                   else
  290.                      write(ioimp,*) 'Inconnue primale de nom : '
  291.      $                    ,des1.lisinc(iligp1),' noeud : '
  292.      $                    ,des1.noelep(iligp1)
  293.                      write(ioimp,*) ' presente dans la rigidite : ',ri1
  294.                      write(ioimp,*) ' non presente dans la rigidite : '
  295.      $                    ,ri3
  296.                      write(ioimp,*) 'des1'
  297.                      segprt,des1
  298.                      write(ioimp,*) 'des3'
  299.                      segprt,des3
  300.                      call erreur(21)
  301.                      return
  302.                   endif
  303.                ELSE
  304.                   IDDLP(ILIGP1)=ilig3
  305.                ENDIF
  306.             ENDDO
  307.             JDDLP(irig1)=IDDLP
  308. *            segprt,iddlp
  309. *       pour la duale
  310.             NLIGD1=DES1.NOELED(/1)
  311.             NLIGRD=NLIGD1
  312.             SEGINI IDDLD
  313.             DO ILIGD1=1,NLIGD1
  314.                nno1=des1.noeled(iligd1)
  315.                ilig3=0
  316.                nligd3=des3.noeled(/1)
  317.                do 213 iligd3=1,nligd3
  318.                   if (des3.noeled(iligd3).ne.nno1) goto 213
  319.                   if (des3.lisdua(iligd3).ne.des1.lisdua(iligd1))
  320.      $                 goto 213
  321.                   ilig3=iligd3
  322.                   goto 214
  323.  213           continue
  324.  214           continue
  325.                IF (ilig3.EQ.0) THEN
  326.                   if (ldescr) then
  327.                      NLIGD3=NLIGD3+1
  328.                      NLIGRD=NLIGD3
  329.                      NLIGRP=des3.noelep(/1)
  330.                      SEGADJ DES3
  331.                      des3.noeled(nligd3)=des1.noeled(iligd1)
  332.                      des3.lisdua(nligd3)=des1.lisdua(iligd1)
  333.                      IDDLD(iligd1)=NLIGD3
  334.                   else
  335.                      write(ioimp,*) 'Inconnue duale de nom : '
  336.      $                    ,des1.lisdua(iligd1),' noeud : '
  337.      $                    ,des1.noeled(iligd1)
  338.                      write(ioimp,*) ' presente dans la rigidite : ',ri1
  339.                      write(ioimp,*) ' non presente dans la rigidite : '
  340.      $                    ,ri3
  341.                      write(ioimp,*) 'des1'
  342.                      segprt,des1
  343.                      write(ioimp,*) 'des3'
  344.                      segprt,des3
  345.                      call erreur(21)
  346.                      return
  347.                   endif
  348.                ELSE
  349.                   IDDLD(iligd1)=ILIG3
  350.                ENDIF
  351.             ENDDO
  352.             JDDLD(irig1)=IDDLD
  353. *            segprt,iddld
  354.          endif
  355.       ENDDO
  356. *      do irig3=1,nrig3
  357. *         des3=ri3.irigel(3,irig3)
  358. *         write(ioimp,*) 'irig3=',irig3
  359. *         segprt,des3
  360. *      enddo
  361. *
  362. * Transfert des xmatr selon le descripteur final des3
  363. *
  364.       DO IRIG1=1,nrig1
  365.          ilmasq=jlmasq(irig1)
  366. *         if (irig1.EQ.1) write(ioimp,*)  'ilmasq=',ilmasq
  367. *         if (irig1.EQ.1) segprt,ilmasq
  368.          if (ilmasq.ne.0) then
  369.             izon=jzon(irig1)
  370. *            if (irig1.EQ.1) write(ioimp,*)  'izon=',izon
  371.             iddlp=jddlp(irig1)
  372. *            if (irig1.EQ.1) segprt,iddlp
  373.             iddld=jddld(irig1)
  374. *            if (irig1.EQ.1) segprt,iddld
  375.             xmatr3=ri3.irigel(4,izon)
  376.             xc1=RI1.COERIG(IRIG1)
  377. *            if (irig1.EQ.1) write(ioimp,*) 'xc1=',xc1
  378.             XMATR1=RI1.IRIGEL(4,IRIG1)
  379.             SEGACT XMATR1
  380.             if (xmatr3.eq.0) then
  381.                ipt3=ri3.irigel(1,izon)
  382.                des3=ri3.irigel(3,izon)
  383.                NLIGRD=des3.NOELED(/1)
  384.                NLIGRP=des3.NOELEP(/1)
  385.                NELRIG=ipt3.num(/2)
  386.                SEGINI XMATR3
  387.                ri3.irigel(4,izon)=xmatr3
  388.                ri3.irigel(7,izon)=max(ri3.irigel(7,izon),ri1.irigel(7
  389.      $              ,irig1))
  390.                xmatr3.symre=max(xmatr3.symre,xmatr1.symre)
  391.             endif
  392.             NLIGD1=XMATR1.RE(/1)
  393.             NLIGP1=XMATR1.RE(/2)
  394.             NELE1 =XMATR1.RE(/3)
  395.             if (ilmasq.gt.0) then
  396.                DO iele1=1,nele1
  397.                   iele3=ilmasq(iele1)
  398.                   IF (iele3.ne.0) THEN
  399.                      DO ILIGP1=1,NLIGP1
  400.                         ILIGP3=IDDLP(ILIGP1)
  401.                         DO ILIGD1=1,NLIGD1
  402.                            ILIGD3=IDDLD(ILIGD1)
  403.                            XMATR3.RE(ILIGD3,ILIGP3,IELE3)=
  404.      $                          XMATR3.RE(ILIGD3,ILIGP3,IELE3)
  405.      $                          +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1
  406.                         ENDDO
  407.                      ENDDO
  408.                   ENDIF
  409.                ENDDO
  410.             else
  411.                DO iele1=1,nele1
  412.                   DO ILIGP1=1,NLIGP1
  413.                      ILIGP3=IDDLP(ILIGP1)
  414.                      DO ILIGD1=1,NLIGD1
  415.                         ILIGD3=IDDLD(ILIGD1)
  416.                         XMATR3.RE(ILIGD3,ILIGP3,IELE1)=
  417.      $                       XMATR3.RE(ILIGD3,ILIGP3,IELE1)
  418.      $                       +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1
  419.                      ENDDO
  420.                   ENDDO
  421.                ENDDO
  422.             endif
  423.          endif
  424.       ENDDO
  425. *      do irig3=1,nrig3
  426. *         xmatr3=ri3.irigel(4,irig3)
  427. *         write(ioimp,*) 'irig3=',irig3
  428. *         if (xmatr3.gt.0) then
  429. *            segprt,xmatr3
  430. *         else
  431. *            write(ioimp,*) 'xmatr3=',xmatr3
  432. *         endif
  433. *      enddo
  434. *  Menage
  435.       DO IRIG1=1,nrig1
  436.          ilmasq=jlmasq(irig1)
  437.          if (ilmasq.gt.0) segsup ilmasq
  438.          iddlp=jddlp(irig1)
  439.          if (iddlp.gt.0) segsup iddlp
  440.          iddld=jddld(irig1)
  441.          if (iddld.gt.0) segsup iddld
  442.       ENDDO
  443.       SEGSUP INFRIG
  444.       RETURN
  445.       END
  446.  
  447.  

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