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    PV090527  26/04/30    21:15:20     12529          
  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.                rigrel=0
  379.                SEGINI XMATR3
  380.                ri3.irigel(4,izon)=xmatr3
  381.                ri3.irigel(7,izon)=max(ri3.irigel(7,izon),ri1.irigel(7
  382.      $              ,irig1))
  383.                xmatr3.symre=max(xmatr3.symre,xmatr1.symre)
  384.             endif
  385.             NLIGD1=XMATR1.RE(/1)
  386.             NLIGP1=XMATR1.RE(/2)
  387.             NELE1 =XMATR1.RE(/3)
  388.             if (ielzo(izon).gt.0) then
  389.                DO iele1=1,nele1
  390.                   iele3=ilmasq(iele1,izon)
  391.                   IF (iele3.ne.0) THEN
  392.                      DO ILIGP1=1,NLIGP1
  393.                         ILIGP3=IDDLP(ILIGP1)
  394.                         DO ILIGD1=1,NLIGD1
  395.                            ILIGD3=IDDLD(ILIGD1)
  396.                            XMATR3.RE(ILIGD3,ILIGP3,IELE3)=
  397.      $                          XMATR3.RE(ILIGD3,ILIGP3,IELE3)
  398.      $                          +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1
  399.                         ENDDO
  400.                      ENDDO
  401.                   ENDIF
  402.                ENDDO
  403.             else
  404.                DO iele1=1,nele1
  405.                   DO ILIGP1=1,NLIGP1
  406.                      ILIGP3=IDDLP(ILIGP1)
  407.                      DO ILIGD1=1,NLIGD1
  408.                         ILIGD3=IDDLD(ILIGD1)
  409.                         XMATR3.RE(ILIGD3,ILIGP3,IELE1)=
  410.      $                       XMATR3.RE(ILIGD3,ILIGP3,IELE1)
  411.      $                       +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1
  412.                      ENDDO
  413.                   ENDDO
  414.                ENDDO
  415.             endif
  416.  31      continue
  417.       ENDDO
  418. *      do irig3=1,nrig3
  419. *         xmatr3=ri3.irigel(4,irig3)
  420. *         write(ioimp,*) 'irig3=',irig3
  421. *         if (xmatr3.gt.0) then
  422. *            segprt,xmatr3
  423. *         else
  424. *            write(ioimp,*) 'xmatr3=',xmatr3
  425. *         endif
  426. *      enddo
  427. *  Menage
  428.       DO IRIG1=1,nrig1
  429.          infzon=jzon(irig1)
  430.          do 41 izon=1,nbsous
  431.             if (ielzo(izon).eq.0) goto 41
  432.             iddlp=jddlp(izon)
  433.             if (iddlp.gt.0) segsup iddlp
  434.             iddld=jddld(izon)
  435.             if (iddld.gt.0) segsup iddld
  436.  41      continue
  437.          segsup infzon
  438.       ENDDO
  439.       SEGSUP INFRIG
  440.       segsup meleme
  441.       segsup jelnum,izone
  442.       segsup inode
  443.       segsup icpr
  444.       RETURN
  445.       END
  446.  
  447.  
  448.  
  449.  

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