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/03/17    21:15:03     12207          
  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.  13         continue
  232.  11      continue
  233.  14      continue
  234.          jzon(irig1)=izon
  235.          if (izon.eq.0.or.ielzo.eq.0) then
  236.             segsup ilmasq
  237.             ilmasq=0
  238.          endif
  239.          jlmasq(irig1)=ilmasq
  240.          if (ilmasq.ne.0) then
  241. *            if (ilmasq.gt.0) then
  242. *               segprt,ilmasq
  243. *            else
  244. *               write(ioimp,*) 'ilmasq=',ilmasq
  245. *            endif
  246. *     tableau de correspondance des DESCR
  247.             DES3=RI3.IRIGEL(3,IZON)
  248.             DES1=RI1.IRIGEL(3,IRIG1)
  249.             SEGACT DES1
  250. *       pour la primale
  251.             NLIGP1=DES1.NOELEP(/1)
  252.             NLIGRP=NLIGP1
  253.             SEGINI IDDLP
  254.             DO ILIGP1=1,NLIGP1
  255.                nno1=des1.noelep(iligp1)
  256.                ilig3=0
  257.                nligp3=des3.noelep(/1)
  258.                do 113 iligp3=1,nligp3
  259.                   if (des3.noelep(iligp3).ne.nno1) goto 113
  260.                   if (des3.lisinc(iligp3).ne.des1.lisinc(iligp1))
  261.      $                 goto 113
  262.                   ilig3=iligp3
  263.                   goto 114
  264.  113           continue
  265.  114           continue
  266.                IF (ilig3.EQ.0) THEN
  267.                   if (ldescr) then
  268.                      NLIGP3=NLIGP3+1
  269.                      NLIGRP=NLIGP3
  270.                      NLIGRD=des3.noeled(/1)
  271.                      SEGADJ DES3
  272.                      des3.noelep(nligp3)=des1.noelep(iligp1)
  273.                      des3.lisinc(nligp3)=des1.lisinc(iligp1)
  274.                      IDDLP(ILIGP1)=nligp3
  275.                   else
  276.                      write(ioimp,*) 'Inconnue de nom : '
  277.      $                    ,des1.lisinc(iligp1),' noeud : '
  278.      $                    ,des1.noelep(iligp1)
  279.                      write(ioimp,*) ' presente dans la rigidite : ',ri1
  280.                      write(ioimp,*) ' non presente dans la rigidite : '
  281.      $                    ,ri3
  282.                      call erreur(21)
  283.                      return
  284.                   endif
  285.                ELSE
  286.                   IDDLP(ILIGP1)=ilig3
  287.                ENDIF
  288.             ENDDO
  289.             JDDLP(irig1)=IDDLP
  290. *            segprt,iddlp
  291. *       pour la duale
  292.             NLIGD1=DES1.NOELED(/1)
  293.             NLIGRD=NLIGD1
  294.             SEGINI IDDLD
  295.             DO ILIGD1=1,NLIGD1
  296.                nno1=des1.noeled(iligd1)
  297.                ilig3=0
  298.                nligd3=des3.noeled(/1)
  299.                do 213 iligd3=1,nligd3
  300.                   if (des3.noeled(iligd3).ne.nno1) goto 213
  301.                   if (des3.lisdua(iligd3).ne.des1.lisdua(iligd1))
  302.      $                 goto 213
  303.                   ilig3=iligd3
  304.                   goto 214
  305.  213           continue
  306.  214           continue
  307.                IF (ilig3.EQ.0) THEN
  308.                   if (ldescr) then
  309.                      NLIGD3=NLIGD3+1
  310.                      NLIGRD=NLIGD3
  311.                      NLIGRP=des3.noelep(/1)
  312.                      SEGADJ DES3
  313.                      des3.noeled(nligd3)=des1.noeled(iligd1)
  314.                      des3.lisdua(nligd3)=des1.lisdua(iligd1)
  315.                      IDDLD(iligd1)=NLIGD3
  316.                   else
  317.                      write(ioimp,*) 'Inconnue de nom : '
  318.      $                    ,des1.lisdua(iligd1),' noeud : '
  319.      $                    ,des1.noeled(iligd1)
  320.                      write(ioimp,*) ' presente dans la rigidite : ',ri1
  321.                      write(ioimp,*) ' non presente dans la rigidite : '
  322.      $                    ,ri3
  323.                      call erreur(21)
  324.                      return
  325.                   endif
  326.                ELSE
  327.                   IDDLD(iligd1)=ILIG3
  328.                ENDIF
  329.             ENDDO
  330.             JDDLD(irig1)=IDDLD
  331. *            segprt,iddld
  332.          endif
  333.       ENDDO
  334. *      do irig3=1,nrig3
  335. *         des3=ri3.irigel(3,irig3)
  336. *         write(ioimp,*) 'irig3=',irig3
  337. *         segprt,des3
  338. *      enddo
  339. *
  340. * Transfert des xmatr selon le descripteur final des3
  341. *
  342.       DO IRIG1=1,nrig1
  343.          ilmasq=jlmasq(irig1)
  344.          if (ilmasq.ne.0) then
  345.             izon=jzon(irig1)
  346.             iddlp=jddlp(irig1)
  347.             iddld=jddld(irig1)
  348.             xmatr3=ri3.irigel(4,izon)
  349.             xc1=RI1.COERIG(IRIG1)
  350.             XMATR1=RI1.IRIGEL(4,IRIG1)
  351.             SEGACT XMATR1
  352.             if (xmatr3.eq.0) then
  353.                ipt3=ri3.irigel(1,izon)
  354.                des3=ri3.irigel(3,izon)
  355.                NLIGRD=des3.NOELED(/1)
  356.                NLIGRP=des3.NOELEP(/1)
  357.                NELRIG=ipt3.num(/2)
  358.                SEGINI XMATR3
  359.                ri3.irigel(4,izon)=xmatr3
  360.                ri3.irigel(7,izon)=max(ri3.irigel(7,izon),ri1.irigel(7
  361.      $              ,irig1))
  362.                xmatr3.symre=max(xmatr3.symre,xmatr1.symre)
  363.             endif
  364.             NLIGD1=XMATR1.RE(/1)
  365.             NLIGP1=XMATR1.RE(/2)
  366.             NELE1 =XMATR1.RE(/3)
  367.             if (ilmasq.gt.0) then
  368.                DO iele1=1,nele1
  369.                   iele3=ilmasq(iele1)
  370.                   IF (iele3.ne.0) THEN
  371.                      DO ILIGP1=1,NLIGP1
  372.                         ILIGP3=IDDLP(ILIGP1)
  373.                         DO ILIGD1=1,NLIGD1
  374.                            ILIGD3=IDDLD(ILIGD1)
  375.                            XMATR3.RE(ILIGD3,ILIGP3,IELE3)=
  376.      $                          XMATR3.RE(ILIGD3,ILIGP3,IELE3)
  377.      $                          +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1
  378.                         ENDDO
  379.                      ENDDO
  380.                   ENDIF
  381.                ENDDO
  382.             else
  383.                DO iele1=1,nele1
  384.                   DO ILIGP1=1,NLIGP1
  385.                      ILIGP3=IDDLP(ILIGP1)
  386.                      DO ILIGD1=1,NLIGD1
  387.                         ILIGD3=IDDLD(ILIGD1)
  388.                         XMATR3.RE(ILIGD3,ILIGP3,IELE1)=
  389.      $                       XMATR3.RE(ILIGD3,ILIGP3,IELE1)
  390.      $                       +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1
  391.                      ENDDO
  392.                   ENDDO
  393.                ENDDO
  394.             endif
  395.          endif
  396.       ENDDO
  397. *      do irig3=1,nrig3
  398. *         xmatr3=ri3.irigel(4,irig3)
  399. *         write(ioimp,*) 'irig3=',irig3
  400. *         if (xmatr3.gt.0) then
  401. *            segprt,xmatr3
  402. *         else
  403. *            write(ioimp,*) 'xmatr3=',xmatr3
  404. *         endif
  405. *      enddo
  406. *  Menage
  407.       DO IRIG1=1,nrig1
  408.          ilmasq=jlmasq(irig1)
  409.          if (ilmasq.gt.0) segsup ilmasq
  410.          iddlp=jddlp(irig1)
  411.          if (iddlp.gt.0) segsup iddlp
  412.          iddld=jddld(irig1)
  413.          if (iddld.gt.0) segsup iddld
  414.       ENDDO
  415.       SEGSUP INFRIG
  416.       RETURN
  417.       END
  418.  
  419.  

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