Télécharger la plaquette Cast3M


fr en

Télécharger cmct3b.eso

Retour à la liste

Numérotation des lignes : 

cmct3b

  1. C CMCT3B    SOURCE    GOUNAND   25/06/11    21:15:03     12278          
  2.       SUBROUTINE CMCT3B(DESCOB,LSINCO,LTINCO,LSINCD,IRIG2)
  3. *_______________________________________________________________________
  4. c
  5. c      opérateur cmct
  6. c
  7. c entrée
  8. c   ICHP  : champ par point qui stocke la masse inversée M-1
  9. c   IRIGB : rigidité B
  10. c   IRIGC : rigidité C
  11. c
  12. c sortie
  13. c   IRIG2 : rigidité contenant la matrice condensée C M-1 Bt
  14. c
  15. *_______________________________________________________________________
  16.  
  17.       IMPLICIT INTEGER(I-N)
  18.       IMPLICIT REAL*8(A-H,O-Z)
  19. *
  20.  
  21. -INC PPARAM
  22. -INC CCOPTIO
  23. -INC SMRIGID
  24. -INC SMELEME
  25. -INC CCHAMP
  26. -INC SMCOORD
  27. -INC SMCHPOI
  28. *
  29. * Description des objets à traiter
  30. * PNTOB : pointeur de l'objet
  31. * TYPOB : type de l'objet
  32. * PNTCOB : pointeur vers une version locale de l'objet
  33. *          pour la RIGIDITE IMAT
  34. *                                PNTCOB(1,IMAT) pointe vers un CORES1
  35. *                                PNTCOB(2,IMAT) pointe vers un LSINCO
  36. *                                PNTCOB(3,IMAT) pointe vers un MCOEF
  37.  
  38.       SEGMENT DESCOB
  39.         INTEGER PNTOB(NMAT)
  40.         INTEGER PNTCOB(3,NMAT)
  41.       ENDSEGMENT
  42. *  stockage des noms des composantes duales.
  43.       SEGMENT LSINCD
  44.         CHARACTER*(LOCOMP) LISIND(NLIGD)
  45.       ENDSEGMENT
  46. *
  47. * Noms d'inconnues primales correspondant à LISIND
  48. *
  49.       SEGMENT CORIDP
  50.         CHARACTER*(LOCOMP) LISIDP(NLIGD)
  51.       ENDSEGMENT
  52. *
  53. * tableau pour pointer vers MCOEF à partir du nombre d'inconnues
  54. *
  55.       SEGMENT LSINCO
  56.         INTEGER LESINC(NINC+1,2,NMAT)
  57.       ENDSEGMENT
  58.       SEGMENT LTINCO
  59.         REAL*8  XMAS(NINC)
  60.       ENDSEGMENT
  61. *
  62. *  tableau des coefficient de la matrice C
  63. *  ordonné dans l'ordre des inconnues
  64.       SEGMENT MCOEF
  65. *        numero du noeud support du multiplicateur ligne 1
  66.       INTEGER ICOEF(2,NCOEF)
  67. *         valeur des coefficients
  68.         REAL*8 XCOEF(NCOEF)
  69.       ENDSEGMENT
  70. *
  71.       SEGMENT WORK1
  72.         REAL*8 XDUM(NBNN)
  73.       ENDSEGMENT
  74. *
  75.       segment inddl(ninc)
  76.       segment icat(0)
  77.       segment ivcat(jg)
  78.       segment incat(ncat)
  79.       segment jcat(ncat+1)
  80.       segment kcat(ninc)
  81.       segment khash(ninc)
  82. * dans ce segment, valeurs positives => iinc uniques dans knuniq
  83. *     valeurs negatives => iinc en doublon valeur absolue pointe sur le
  84. *     unique dans knuniq
  85.       segment kuniq(ninc)
  86. * nombre de iinc uniques dans chaque categorie
  87.       segment inuniq(ncat)
  88. * stokage des iinc
  89.       segment knuniq(nuniq)
  90. * nombre de iinc partageant un meme unique
  91.       segment nnuniq(nuniq)
  92. * un segment de correspondance
  93.       segment jnuniq(nuniq)
  94.  
  95.       LOGICAL LDBG,LFOUND
  96. *
  97.       LDBG=.FALSE.
  98. *_______________________________________________________________________
  99. *
  100. *   il ne reste plus qu' a creer les matrices élémentaires
  101. *
  102. * Construction des noms d'inconnues primales correspondant à LISIND
  103.       NLIGD=LISIND(/2)
  104.       SEGINI CORIDP
  105.       DO ILIGD=1,NLIGD
  106.          CALL PLACE(NOMDU,LNOMDU,idx,LISIND(ILIGD))
  107.          IF (idx.EQ.0) THEN
  108.             LISIDP(ILIGD)=LISIND(ILIGD)
  109.          ELSE
  110.             LISIDP(ILIGD)=NOMDD(idx)
  111.          ENDIF
  112.       ENDDO
  113. *
  114. *     Pour diminuer le nombre de NRIGEL a creer, on va chercher les
  115. *     redondances dans le segment LESINC(I,1,IMAT), d'abord par nombre
  116. *     de ddls (LESINC(I,2,IMAT)) puis par egalite des suites de ddls
  117. *
  118.       NINC=LESINC(/1)-1
  119.       NMAT=LESINC(/3)
  120.       MAXDDL=0
  121.       segini inddl
  122.       do IINC=1,NINC
  123.          do IMAT=NMAT,1,-1
  124.             inddl(iinc)=inddl(iinc)+LESINC(IINC,2,IMAT)
  125.          enddo
  126.          MAXDDL=max(maxddl,inddl(iinc))
  127.       enddo
  128. *
  129.       segini icat
  130.       jg=maxddl
  131.       segini ivcat
  132.       ncat=maxddl
  133.       segini incat
  134.       do IINC=1,NINC
  135.          NDDL=inddl(iinc)
  136.          if (incat(nddl).eq.0) then
  137.             icat(**)=nddl
  138.             ivcat(nddl)=icat(/1)
  139.          endif
  140.          incat(nddl)=incat(nddl)+1
  141.       enddo
  142.       IF (LDBG) THEN
  143.          write(ioimp,*) 'Classement par categorie (nombre de ddls)'
  144.          write(ioimp,*) '  Categories'
  145.          segprt,icat
  146.          segprt,ivcat
  147.          write(ioimp,*) '  Nombre par categorie'
  148.          segprt,incat
  149.       ENDIF
  150.       ncat=icat(/1)
  151.       segini jcat
  152.       jcat(1)=1
  153.       do ic=1,ncat
  154.          jcat(ic+1)=jcat(ic)+incat(icat(ic))
  155.          incat(icat(ic))=0
  156.       enddo
  157.       segsup icat
  158.       segadj incat
  159. * verif
  160.       if (jcat(ncat+1)-1.ne.ninc) then
  161.          write(ioimp,*) 'ncat,ninc=',ncat,ninc
  162.          segprt,jcat
  163.          call erreur(5)
  164.          return
  165.       endif
  166.       segini kcat
  167.       do IINC=1,NINC
  168.          nucat=ivcat(inddl(iinc))
  169.          idx=jcat(nucat)+incat(nucat)
  170.          kcat(idx)=iinc
  171.          incat(nucat)=incat(nucat)+1
  172.       enddo
  173.       segsup incat
  174.       segsup ivcat
  175.       segsup inddl
  176. *
  177. * Classement par categorie (nombre de ddls )
  178. *
  179.       IF (LDBG) THEN
  180.          write(ioimp,*) '  Classement des inconnues duales'
  181.          segprt,jcat
  182.          segprt,kcat
  183.       ENDIF
  184. * Un hash pour trier plus vite les non-identiques
  185.       segini khash
  186.       do IINC=1,NINC
  187.          IHASH=0
  188.          DO IMAT=NMAT,1,-1
  189.             MCOEF=PNTCOB(3,IMAT)
  190.             IDEB=LESINC(IINC,1,IMAT)
  191.             IFIN=IDEB+LESINC(IINC,2,IMAT)-1
  192.             DO IDX=IDEB,IFIN
  193.                IHASH=IHASH*(NLIGD+1)+ICOEF(2,IDX)
  194.             ENDDO
  195.          ENDDO
  196.          khash(IINC)=IHASH
  197.       enddo
  198.       IF (LDBG) THEN
  199.          write(ioimp,*) 'Hash pour les inconnues duales'
  200.          segprt,khash
  201.       ENDIF
  202.       segini kuniq
  203.       segini inuniq
  204.       nuniq=ninc
  205.       segini knuniq
  206.       segini nnuniq
  207.       do ic=1,ncat
  208.          idxdeb=jcat(ic)
  209.          idxfin=jcat(ic+1)-1
  210.          do j=idxdeb,idxfin
  211. * Le premier de la passe est forcement unique
  212.             iinc=kcat(j)
  213.             lfound=.false.
  214.             if (j.ne.idxdeb) then
  215.                ihash=khash(iinc)
  216.                do 32 k=1,inuniq(ic)
  217.                   juniq=idxdeb+k-1
  218.                   iinck=knuniq(juniq)
  219.                   ihashk=khash(iinck)
  220.                   if (ihashk.ne.ihash) goto 32
  221.                   do imat=nmat,1,-1
  222.                      mcoef=pntcob(3,imat)
  223.                      kdeb=lesinc(iinc,1,imat)
  224.                      knb =lesinc(iinc,2,imat)
  225.                      kdebk=lesinc(iinck,1,imat)
  226.                      do 33 inb=1,knb
  227.                         icoefu=icoef(2,kdeb+inb-1)
  228.                         icoefk=icoef(2,kdebk+inb-1)
  229.                         if (icoefu.ne.icoefk) goto 32
  230.  33                  continue
  231.                   enddo
  232.                   lfound=.true.
  233.                   goto 34
  234.  32            continue
  235.  34            continue
  236.             endif
  237.             if (.not.lfound) then
  238.                juniq=idxdeb+inuniq(ic)
  239.                inuniq(ic)=inuniq(ic)+1
  240.                knuniq(juniq)=iinc
  241.                nnuniq(juniq)=1
  242.                kuniq(j)=juniq
  243. *               kuniq(iinc)=juniq
  244.             else
  245.                kuniq(j)=-juniq
  246. *               kuniq(iinc)=-juniq
  247.                nnuniq(juniq)=nnuniq(juniq)+1
  248.             endif
  249.          enddo
  250.       enddo
  251.       segsup khash
  252.       IF (LDBG) THEN
  253.          write(ioimp,*) 'Primales unique par duale'
  254.          write(ioimp,*) '  Nombre de groupes de primales/duale'
  255.          segprt,inuniq
  256.          write(ioimp,*) '  Stockage des duales uniques'
  257.          segprt,knuniq
  258.          write(ioimp,*) '  Combien partagent un meme uniq'
  259.          segprt,nnuniq
  260.          write(ioimp,*) '  Pointeur pour chaque duale vers l''unique'
  261.          segprt,kuniq
  262.       ENDIF
  263.       NRIGEL=0
  264.       do i=1,ncat
  265.          nrigel=nrigel+inuniq(i)
  266.       enddo
  267.       IF (LDBG) THEN
  268.          write(ioimp,2030) 'Nombre de sous-zones =',NRIGEL
  269.       ENDIF
  270. *
  271. *     Compacter knuniq, nnuniq et mettre à jour kuniq
  272. *     juniq varie au debut entre 1 et ninc puis apres entre 1 et nrigel
  273. *     On utilise jnuniq comme tableau de conversion juniq_avant juniq_apres
  274. *
  275. *      if (.false.) then
  276.       if (nrigel.lt.ninc) then
  277.          segini jnuniq
  278.          juniqp=0
  279.          do juniqv=1,ninc
  280.             if (nnuniq(juniqv).ne.0) then
  281.                juniqp=juniqp+1
  282.                nnuniq(juniqp)=nnuniq(juniqv)
  283.                knuniq(juniqp)=knuniq(juniqv)
  284.                jnuniq(juniqv)=juniqp
  285.             endif
  286.          enddo
  287.          do iinc=1,ninc
  288.             if (kuniq(iinc).ge.0) then
  289.                kuniq(iinc)=jnuniq(kuniq(iinc))
  290.             else
  291.                ixxx=-kuniq(iinc)
  292.                kuniq(iinc)=-jnuniq(ixxx)
  293.             endif
  294.          enddo
  295.          if (juniqp.ne.nrigel) then
  296.             write(ioimp,*) 'juniqp.ne.nrigel',juniqp,nrigel
  297.             call erreur(5)
  298.             return
  299.          endif
  300.          nuniq=nrigel
  301.          segadj,knuniq
  302.          segadj,nnuniq
  303.          segsup jnuniq
  304.          IF (LDBG) THEN
  305.             write(ioimp,*) 'Apres compactage'
  306.             write(ioimp,*) '  Stockage des duales uniques'
  307.             segprt,knuniq
  308.             write(ioimp,*) '  Combien partagent un meme uniq'
  309.             segprt,nnuniq
  310.             write(ioimp,*) '  Pointeur pour chaque duale vers l''unique'
  311.             segprt,kuniq
  312.          ENDIF
  313.       endif
  314. *
  315.       SEGINI MRIGID
  316.       IRIG2 = MRIGID
  317.       MTYMAT = 'CMCT    '
  318.       IFORIG = IFOUR
  319. *     Pour le mode de calcul on devrait s'appuyer sur celui des matrices/champs en entree ?
  320. *      VRAI ! A faire
  321.  
  322. *
  323. *    boucle sur les categories
  324. *
  325.       NBNN=maxddl
  326.       SEGINI WORK1
  327. *!      segini jnuniq
  328.       iricnt=0
  329.  2030 FORMAT (A,20(2X,I5))
  330.  
  331.       do ic=1,ncat
  332. *!         idxdeb=jcat(ic)
  333.          nrigun=inuniq(ic)
  334.          do irigun=1,nrigun
  335.             iricnt=iricnt+1
  336.             coerig(iricnt)=1.D0
  337. *!            juniq=idxdeb+irigun-1
  338.             juniq=iricnt
  339.             nbl=nnuniq(juniq)
  340.             iinc=knuniq(juniq)
  341.             if (LDBG) THEN
  342. *!               write(ioimp,2030)
  343. *!     $              'ic,idxdeb,irigun,iricnt,juniq,nbl,iinc=',ic
  344. *!     $              ,idxdeb,irigun,iricnt,juniq,nbl,iinc
  345.                write(ioimp,2030)
  346.      $              'ic,irigun,iricnt,juniq,nbl,iinc=',ic
  347.      $              ,irigun,iricnt,juniq,nbl,iinc
  348.             endif
  349.             NBNN=0
  350.             IDXDUA=0
  351. *
  352. * S'il y a deux rigidités à multiplier, on stocke la primale
  353. * puis la duale dans le MELEME et dans XDUM. A ce moment-là,
  354. * IDXDUA+1 est l'index de départ des noeuds et des valeurs de la duale.
  355. * On parcourt donc DESCOB à l'envers car on avait d'abord mis
  356. * la duale C avant la primale B
  357. *
  358.             DO IMAT=NMAT,1,-1
  359.                NBNN=NBNN+LESINC(iinc,2,IMAT)
  360.                IF (IDXDUA.EQ.0.AND.NMAT.EQ.2) IDXDUA=NBNN
  361.             ENDDO
  362. * creation du maillage et du vecteur des coefficients
  363.             NBELEM = nbl
  364.             NBSOUS = 0
  365.             NBREF = 0
  366.             SEGINI MELEME
  367.             ITYPEL = 32
  368.             IRIGEL(1,iricnt) = MELEME
  369. *
  370. *    segment descripteur DESCR
  371. *
  372.             NLIGRP=NBNN
  373.             IF (NMAT.EQ.1) THEN
  374.                NLIGRP = NBNN
  375.                NLIGRD = NBNN
  376.             ELSE
  377.                NLIGRP = IDXDUA
  378.                NLIGRD = NBNN-IDXDUA
  379.             ENDIF
  380.             SEGINI DESCR
  381.             IRIGEL(3,iricnt) = DESCR
  382. *
  383.             NELRIG =nbl
  384.             SEGINI XMATRI
  385.             IRIGEL(4,iricnt)=XMATRI
  386. *     S'il y a deux matrices en entrée, le résultat
  387. *     n'est sans doute pas symétrique
  388.             IF (NMAT.EQ.2) then
  389.                IRIGEL(7,iricnt)=2
  390.                xmatri.symre=2
  391.             endif
  392. *! On recycle nnuniq pour stocker iricnt
  393. *!            nnuniq(juniq)=iricnt
  394. * On recycle nnuniq comme compteur
  395.             nnuniq(juniq)=0
  396.          enddo
  397. *
  398. *     Remplissage des objets
  399. *
  400.          idxdeb=jcat(ic)
  401.          idxfin=jcat(ic+1)-1
  402.          do jc=idxdeb,idxfin
  403.             iinc=kcat(jc)
  404. *     lfirst=kuniq(iinc).gt.0
  405.             juniq=abs(kuniq(jc))
  406.             iricou=juniq
  407. *!            iricou=nnuniq(juniq)
  408.             meleme=irigel(1,iricou)
  409.             descr=irigel(3,iricou)
  410.             xmatri=irigel(4,iricou)
  411.             iincu=knuniq(juniq)
  412. *!            jnuniq(juniq)=jnuniq(juniq)+1
  413. *!            ibl=jnuniq(juniq)
  414.             nnuniq(juniq)=nnuniq(juniq)+1
  415.             ibl=nnuniq(juniq)
  416.             if (LDBG) THEN
  417.                write(ioimp,2030)
  418.      $              'ic,jc,iinc,juniq,iincu,ibl,iricou='
  419.      $              ,ic,jc,iinc,juniq,iincu,ibl,iricou
  420.             endif
  421. *     remplissage MELEME
  422.             INOEU = 0
  423.             DO IMAT=NMAT,1,-1
  424.                MCOEF=PNTCOB(3,IMAT)
  425.                DO 1200 J=0,LESINC(iinc,2,IMAT)-1
  426.                   INOEU = INOEU + 1
  427.                   NUM(INOEU,ibl) = ICOEF(1,J+LESINC(iinc,1,IMAT))
  428.                   XDUM(INOEU) = XCOEF(J+LESINC(iinc,1,IMAT))
  429.  1200          CONTINUE
  430.             ENDDO
  431.  
  432. *     remplissage DESCR (une seule fois, on a tout fait pour)
  433.             if (ibl.eq.1) then
  434.                IPASS=0
  435.                DO IMAT=NMAT,1,-1
  436.                   MCOEF=PNTCOB(3,IMAT)
  437.                   IPASS=IPASS+1
  438.                   INOEU=0
  439.                   DO 1300 J=0,LESINC(iincu,2,IMAT)-1
  440.                      INOEU = INOEU + 1
  441.                      IF (IPASS.EQ.1) THEN
  442.                         LISINC(INOEU) = LISIDP(ICOEF(2,J+LESINC(iincu,1
  443.      $                       ,IMAT)))
  444.                         NOELEP(INOEU)=INOEU
  445.                         IF (NMAT.EQ.1) THEN
  446.                            LISDUA(INOEU) = LISIND(ICOEF(2,J+LESINC(iincu
  447.      $                          ,1,IMAT)))
  448.                            NOELED(INOEU)=INOEU
  449.                         ENDIF
  450.                      ENDIF
  451.                      IF (IPASS.EQ.2) THEN
  452.                         LISDUA(INOEU)
  453.      $                       = LISIND(ICOEF(2,J+LESINC(iincu,1,IMAT)))
  454.                         NOELED(INOEU)=IDXDUA+INOEU
  455.                      ENDIF
  456.  1300             CONTINUE
  457.                ENDDO
  458.             endif
  459. *     remplissage XMATRI
  460.             IF (LTINCO.NE.0) THEN
  461.                XCOEF=XMAS(iinc)
  462.             ELSE
  463.                XCOEF=1.D0
  464.             ENDIF
  465.             nligrd=re(/1)
  466.             nligrp=re(/2)
  467.             if (nmat.eq.2) then
  468.                idxdua=nligrp
  469.             else
  470.                idxdua=0
  471.             endif
  472.             DO 600 J=1,NLIGRP
  473.                DO 500 K=1,NLIGRD
  474.                   RE(K,J,ibl)=XDUM(IDXDUA+K)*XDUM(J)*XCOEF
  475.  500           CONTINUE
  476.  600        CONTINUE
  477. * recyclage XDUM
  478.             INOEU = 0
  479.             DO IMAT=NMAT,1,-1
  480.                MCOEF=PNTCOB(3,IMAT)
  481.                DO 2200 J=0,LESINC(iinc,2,IMAT)-1
  482.                   INOEU = INOEU + 1
  483.                   XDUM(INOEU) = 0.D0
  484.  2200          CONTINUE
  485.             ENDDO
  486.          enddo
  487.       enddo
  488. *!      segsup jnuniq
  489.       segsup nnuniq
  490.       segsup knuniq
  491.       segsup inuniq
  492.       segsup kuniq
  493.       segsup jcat
  494.       segsup kcat
  495.       SEGSUP WORK1
  496.       SEGSUP CORIDP
  497. *
  498. *_______________________________________________________________________
  499.       RETURN
  500.       END
  501.  
  502.  

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