Télécharger la plaquette Cast3M


fr en

Télécharger cmct6.eso

Retour à la liste

Numérotation des lignes : 

cmct6

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

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