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/05/15    21:15:05     12266          
  2.       SUBROUTINE CMCT5(IRI1,INFMEL,IRI3)
  3. *_______________________________________________________________________
  4. c
  5. c      opérateur 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
  56. *
  57. * Recuperation IMEL et objets de preconditionnements eventuels
  58. *
  59.       MELEME=jmel
  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.             do 2 j=1,ipt1.num(/2)
  70.                do 22 k=1,ipt1.num(/1)
  71.                   ia= ipt1.num(k,j)
  72.                   if(icpr(ia).eq.0) then
  73.                      ino=ino+1
  74.                      icpr(ia)=ino
  75.                   endif
  76.  22            continue
  77.  2          continue
  78.  1       continue
  79. *     write(ioimp,*) 'ino=',ino
  80. *     segprt,icpr
  81. *
  82. * on compte combien d'elements touche un noeud
  83. *
  84.          segini inode
  85.          do 3 i=1,nbsou
  86.             ipt1=lisous(i)
  87.             do 4 j=1,ipt1.num(/2)
  88.                do 42 k=1,ipt1.num(/1)
  89.                   ia=ipt1.num(k,j)
  90.                   ib= icpr(ia)
  91. *               write(ioimp,*) 'ia,ib=',ia,ib
  92.                   inode(ib)=inode(ib)+1
  93.  42            continue
  94.  4          continue
  95.  3       continue
  96.          imaxel=0
  97.          do i=1,ino
  98.             imaxel=max(imaxel,inode(i))
  99.             inode(i)=0
  100.          enddo
  101.          segini jelnum,izone
  102. *
  103. * jelnum(i,k) dira le Ieme element qui touche le noeud k
  104. * izone (i,k) dira dans quel lisous le ieme element se trouve
  105. * attention le noeud K est le plus petit de l'élément
  106. *
  107.          do 5 i=1,nbsou
  108.             ipt1=lisous(i)
  109.             do 6 j=1,ipt1.num(/2)
  110.                ipeti=igrand
  111.                do  k=1,ipt1.num(/1)
  112.                   ipeti=min(ipt1.num(k,j),ipeti)
  113.                enddo
  114. *           write(ioimp,*) 'ipeti=',ipeti
  115.                ib= icpr(ipeti)
  116.                inode(ib)=inode(ib)+1
  117.                ic=inode(ib)
  118.                jelnum(ic,ib)=j
  119.                izone(ic,ib)=i
  120.  6          continue
  121.  5       continue
  122. *      segprt,izone
  123. *     segprt,jelnum
  124. * Sauvegarde des objets créés dans INFMEL
  125.          jcpr=icpr
  126.          jnode=inode
  127.          kelnum=jelnum
  128.          jzone=izone
  129.       endif
  130.       icpr=jcpr
  131.       inode=jnode
  132.       jelnum=kelnum
  133.       izone=jzone
  134. *
  135. * Matrice réduite (sortie)
  136. *
  137.       RI1=IRI1
  138.       SEGACT MELEME
  139.       IF (IRI3.EQ.0) THEN
  140.          ldescr=.true.
  141.          NRIG3=LISOUS(/1)
  142.          NRIGEL=NRIG3
  143.          SEGINI RI3
  144.          IRI3=RI3
  145.          RI3.MTYMAT='CMCT5'
  146.          RI3.IFORIG=IFOUR
  147.          RI3.MCRCNF=MCOORD
  148.          NLIGRP=0
  149.          NLIGRD=0
  150.          DO irig3=1,nrig3
  151.             SEGINI,DES3
  152.             RI3.IRIGEL(3,irig3)=DES3
  153.          ENDDO
  154.       ELSE
  155.          ldescr=.false.
  156.          RI3=IRI3
  157.          SEGACT RI3*MOD
  158.          NRIG3=RI3.COERIG(/1)
  159.          IF (NRIG3.NE.LISOUS(/1)) THEN
  160.             write(ioimp,*) 'Erreur grave 1 dans cmct5'
  161.             call erreur(5)
  162.             return
  163.          endif
  164.          DO irig3=1,nrig3
  165.             DES3=RI3.IRIGEL(3,irig3)
  166.             SEGACT,DES3
  167.          ENDDO
  168.       ENDIF
  169.       RI3.MTYMAT='CMCT5'
  170.       RI3.IFORIG=IFOUR
  171.       RI3.MCRCNF=MCOORD
  172.       DO irig3=1,nrig3
  173.          RI3.COERIG(irig3)=1.D0
  174.          ipt3=LISOUS(irig3)
  175.          RI3.IRIGEL(1,irig3)=ipt3
  176.       ENDDO
  177. *
  178. * Traitement de la matrice d'entree
  179. *
  180.       SEGACT RI1
  181.       NRIG1=RI1.COERIG(/1)
  182.       NRIGEL=NRIG1
  183.       SEGINI INFRIG
  184. *
  185.       DO IRIG1=1,nrig1
  186.          IPT1=RI1.IRIGEL(1,IRIG1)
  187. * Y a-t-il des elements a considerer
  188.          SEGACT IPT1
  189.          nele1=ipt1.num(/2)
  190.          segini ilmasq
  191.          nnod1=ipt1.num(/1)
  192.          izon =0
  193.          ielzo=0
  194.          DO 11 iele1=1,nele1
  195.             ipeti=IGRAND
  196.             DO inod1=1,nnod1
  197.                ipeti=min(ipeti,ipt1.num(inod1,iele1))
  198.             enddo
  199. * on regarde s'il existe un element de ipt1 ayant ce noeud en plus
  200. * petit position si non on passe a l'élément suivant
  201.             ib=icpr(ipeti)
  202.             if(ib.eq.0) go to 11
  203.             if(inode(ib).eq.0) go to 11
  204.             do 13 mm=1,inode(ib)
  205.                kzon=izone(mm,ib)
  206.                ipt3 = lisous(kzon)
  207.                if (ipt3.eq.ipt1) then
  208.                   segsup ilmasq
  209.                   ilmasq=-1
  210.                   izon=kzon
  211.                   ielzo=ipt3.num(/2)
  212.                   goto 14
  213.                endif
  214.                if (ipt3.num(/1).ne.ipt1.num(/1)) go to 13
  215.                iel=jelnum(mm,ib)
  216.                do  in=1,ipt1.num(/1)
  217.                   if( ipt1.num(in,iele1).ne.ipt3.num(in,iel))go to 13
  218.                enddo
  219. * on a trouve un element a conserver
  220.                ielzo=ielzo+1
  221.                ilmasq(iele1)=iel
  222.                if (izon.eq.0) then
  223.                   izon=kzon
  224.                else
  225.                   if (izon.ne.kzon) then
  226.                      write(ioimp,*) 'Erreur grave 2 dans cmct5'
  227.                      call erreur(5)
  228.                      return
  229.                   endif
  230.                endif
  231.                goto 11
  232.  13         continue
  233.  11      continue
  234.  14      continue
  235.          jzon(irig1)=izon
  236. *         write(ioimp,*) 'irig1=',irig1,' ielzo=',ielzo
  237.          if (izon.eq.0.or.ielzo.eq.0) then
  238.             segsup ilmasq
  239.             ilmasq=0
  240.          endif
  241.          jlmasq(irig1)=ilmasq
  242.          if (ilmasq.ne.0) then
  243. *            if (ilmasq.gt.0) then
  244. *               segprt,ilmasq
  245. *            else
  246. *               write(ioimp,*) 'ilmasq=',ilmasq
  247. *            endif
  248. *     tableau de correspondance des DESCR
  249.             DES3=RI3.IRIGEL(3,IZON)
  250.             DES1=RI1.IRIGEL(3,IRIG1)
  251.             SEGACT DES1
  252. *       pour la primale
  253.             NLIGP1=DES1.NOELEP(/1)
  254.             NLIGRP=NLIGP1
  255.             SEGINI IDDLP
  256.             DO ILIGP1=1,NLIGP1
  257.                nno1=des1.noelep(iligp1)
  258.                ilig3=0
  259.                nligp3=des3.noelep(/1)
  260.                do 113 iligp3=1,nligp3
  261.                   if (des3.noelep(iligp3).ne.nno1) goto 113
  262.                   if (des3.lisinc(iligp3).ne.des1.lisinc(iligp1))
  263.      $                 goto 113
  264.                   ilig3=iligp3
  265.                   goto 114
  266.  113           continue
  267.  114           continue
  268.                IF (ilig3.EQ.0) THEN
  269.                   if (ldescr) then
  270.                      NLIGP3=NLIGP3+1
  271.                      NLIGRP=NLIGP3
  272.                      NLIGRD=des3.noeled(/1)
  273.                      SEGADJ DES3
  274.                      des3.noelep(nligp3)=des1.noelep(iligp1)
  275.                      des3.lisinc(nligp3)=des1.lisinc(iligp1)
  276.                      IDDLP(ILIGP1)=nligp3
  277.                   else
  278.                      write(ioimp,*) 'Inconnue de nom : '
  279.      $                    ,des1.lisinc(iligp1),' noeud : '
  280.      $                    ,des1.noelep(iligp1)
  281.                      write(ioimp,*) ' presente dans la rigidite : ',ri1
  282.                      write(ioimp,*) ' non presente dans la rigidite : '
  283.      $                    ,ri3
  284.                      call erreur(21)
  285.                      return
  286.                   endif
  287.                ELSE
  288.                   IDDLP(ILIGP1)=ilig3
  289.                ENDIF
  290.             ENDDO
  291.             JDDLP(irig1)=IDDLP
  292. *            segprt,iddlp
  293. *       pour la duale
  294.             NLIGD1=DES1.NOELED(/1)
  295.             NLIGRD=NLIGD1
  296.             SEGINI IDDLD
  297.             DO ILIGD1=1,NLIGD1
  298.                nno1=des1.noeled(iligd1)
  299.                ilig3=0
  300.                nligd3=des3.noeled(/1)
  301.                do 213 iligd3=1,nligd3
  302.                   if (des3.noeled(iligd3).ne.nno1) goto 213
  303.                   if (des3.lisdua(iligd3).ne.des1.lisdua(iligd1))
  304.      $                 goto 213
  305.                   ilig3=iligd3
  306.                   goto 214
  307.  213           continue
  308.  214           continue
  309.                IF (ilig3.EQ.0) THEN
  310.                   if (ldescr) then
  311.                      NLIGD3=NLIGD3+1
  312.                      NLIGRD=NLIGD3
  313.                      NLIGRP=des3.noelep(/1)
  314.                      SEGADJ DES3
  315.                      des3.noeled(nligd3)=des1.noeled(iligd1)
  316.                      des3.lisdua(nligd3)=des1.lisdua(iligd1)
  317.                      IDDLD(iligd1)=NLIGD3
  318.                   else
  319.                      write(ioimp,*) 'Inconnue de nom : '
  320.      $                    ,des1.lisdua(iligd1),' noeud : '
  321.      $                    ,des1.noeled(iligd1)
  322.                      write(ioimp,*) ' presente dans la rigidite : ',ri1
  323.                      write(ioimp,*) ' non presente dans la rigidite : '
  324.      $                    ,ri3
  325.                      call erreur(21)
  326.                      return
  327.                   endif
  328.                ELSE
  329.                   IDDLD(iligd1)=ILIG3
  330.                ENDIF
  331.             ENDDO
  332.             JDDLD(irig1)=IDDLD
  333. *            segprt,iddld
  334.          endif
  335.       ENDDO
  336. *      do irig3=1,nrig3
  337. *         des3=ri3.irigel(3,irig3)
  338. *         write(ioimp,*) 'irig3=',irig3
  339. *         segprt,des3
  340. *      enddo
  341. *
  342. * Transfert des xmatr selon le descripteur final des3
  343. *
  344.       DO IRIG1=1,nrig1
  345.          ilmasq=jlmasq(irig1)
  346. *         if (irig1.EQ.1) write(ioimp,*)  'ilmasq=',ilmasq
  347. *         if (irig1.EQ.1) segprt,ilmasq
  348.          if (ilmasq.ne.0) then
  349.             izon=jzon(irig1)
  350. *            if (irig1.EQ.1) write(ioimp,*)  'izon=',izon
  351.             iddlp=jddlp(irig1)
  352. *            if (irig1.EQ.1) segprt,iddlp
  353.             iddld=jddld(irig1)
  354. *            if (irig1.EQ.1) segprt,iddld
  355.             xmatr3=ri3.irigel(4,izon)
  356.             xc1=RI1.COERIG(IRIG1)
  357. *            if (irig1.EQ.1) write(ioimp,*) 'xc1=',xc1
  358.             XMATR1=RI1.IRIGEL(4,IRIG1)
  359.             SEGACT XMATR1
  360.             if (xmatr3.eq.0) then
  361.                ipt3=ri3.irigel(1,izon)
  362.                des3=ri3.irigel(3,izon)
  363.                NLIGRD=des3.NOELED(/1)
  364.                NLIGRP=des3.NOELEP(/1)
  365.                NELRIG=ipt3.num(/2)
  366.                SEGINI XMATR3
  367.                ri3.irigel(4,izon)=xmatr3
  368.                ri3.irigel(7,izon)=max(ri3.irigel(7,izon),ri1.irigel(7
  369.      $              ,irig1))
  370.                xmatr3.symre=max(xmatr3.symre,xmatr1.symre)
  371.             endif
  372.             NLIGD1=XMATR1.RE(/1)
  373.             NLIGP1=XMATR1.RE(/2)
  374.             NELE1 =XMATR1.RE(/3)
  375.             if (ilmasq.gt.0) then
  376.                DO iele1=1,nele1
  377.                   iele3=ilmasq(iele1)
  378.                   IF (iele3.ne.0) THEN
  379.                      DO ILIGP1=1,NLIGP1
  380.                         ILIGP3=IDDLP(ILIGP1)
  381.                         DO ILIGD1=1,NLIGD1
  382.                            ILIGD3=IDDLD(ILIGD1)
  383.                            XMATR3.RE(ILIGD3,ILIGP3,IELE3)=
  384.      $                          XMATR3.RE(ILIGD3,ILIGP3,IELE3)
  385.      $                          +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1
  386.                         ENDDO
  387.                      ENDDO
  388.                   ENDIF
  389.                ENDDO
  390.             else
  391.                DO iele1=1,nele1
  392.                   DO ILIGP1=1,NLIGP1
  393.                      ILIGP3=IDDLP(ILIGP1)
  394.                      DO ILIGD1=1,NLIGD1
  395.                         ILIGD3=IDDLD(ILIGD1)
  396.                         XMATR3.RE(ILIGD3,ILIGP3,IELE1)=
  397.      $                       XMATR3.RE(ILIGD3,ILIGP3,IELE1)
  398.      $                       +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1
  399.                      ENDDO
  400.                   ENDDO
  401.                ENDDO
  402.             endif
  403.          endif
  404.       ENDDO
  405. *      do irig3=1,nrig3
  406. *         xmatr3=ri3.irigel(4,irig3)
  407. *         write(ioimp,*) 'irig3=',irig3
  408. *         if (xmatr3.gt.0) then
  409. *            segprt,xmatr3
  410. *         else
  411. *            write(ioimp,*) 'xmatr3=',xmatr3
  412. *         endif
  413. *      enddo
  414. *  Menage
  415.       DO IRIG1=1,nrig1
  416.          ilmasq=jlmasq(irig1)
  417.          if (ilmasq.gt.0) segsup ilmasq
  418.          iddlp=jddlp(irig1)
  419.          if (iddlp.gt.0) segsup iddlp
  420.          iddld=jddld(irig1)
  421.          if (iddld.gt.0) segsup iddld
  422.       ENDDO
  423.       SEGSUP INFRIG
  424.       RETURN
  425.       END
  426.  
  427.  

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