Télécharger la plaquette Cast3M


fr en

Télécharger exdiar.eso

Retour à la liste

Numérotation des lignes : 

exdiar

  1. C EXDIAR    SOURCE    CB215821  22/07/20    15:39:40     11411          
  2.       SUBROUTINE EXDIAR(MRIGID,MCHPOI)
  3.       IMPLICIT REAL*8 (A-H,O-Z)
  4.       IMPLICIT INTEGER (I-N)
  5. C*************************************************************************
  6. C     Operateur EXDIAR
  7. C
  8. C     OBJET   : Extrait la diagonale d'une matrice.
  9. C
  10. C***********************************************************************
  11. C HISTORIQUE : 06/06/1 : Première version
  12. C
  13. C***********************************************************************
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17. -INC CCHAMP
  18. -INC SMCHPOI
  19. -INC SMRIGID
  20. -INC SMELEME
  21. -INC SMCOORD
  22. -INC TMTRAV
  23. *
  24.       segment trav1
  25.       character*4 compd(nbincd),compdd(nbincd)
  26.       endsegment
  27.       segment trav2
  28.       integer cor1p(nligrp),cor1d(nligrd)
  29.       endsegment
  30. *
  31. * Pour chaque segment DESCR, on a un segment CORES2 où
  32. *   CORDIA donne les couples (NLIGRP,NLIGRD) qui sont des
  33. *   termes diagonaux avec les cas particuliers suivant :
  34. *   CORES2=0 : il n'y a pas de termes diagonaux dans ces matrices
  35. *      élémentaires
  36. *   CORES2=-1 : les termes diagonaux de la matrices élémentaire (carrée)
  37. *      sont des termes diagonaux de la matrice globale
  38. *   on a aussi a segment CORES3 où COMTRA indique
  39. *   le nom d'inconnue primal+numéro d'harmonique suivant
  40. *   le repérage induit par MTRA1 et MTRA2
  41. *
  42.       SEGMENT CORES1
  43.       INTEGER ITRAV2(NRIGEL)
  44.       INTEGER IPCOR2(NRIGEL)
  45.       INTEGER IPCOR3(NRIGEL)
  46.       ENDSEGMENT
  47.       SEGMENT CORES2
  48.       INTEGER CORDIA(2,NLIGP)
  49.       ENDSEGMENT
  50.       SEGMENT CORES3
  51.       INTEGER COMTRA(NLIGP)
  52.       ENDSEGMENT
  53. *
  54.       SEGMENT ICPR(nbpts)
  55.       SEGMENT MTRA1
  56. *      CHARACTER*4 ICOMP(0)
  57.       INTEGER ICOMP(0)
  58.       ENDSEGMENT
  59.       SEGMENT MTRA2
  60.       INTEGER MHAR(0)
  61.       ENDSEGMENT
  62.       logical descar
  63. *
  64. * Construisons d'abord l'ensemble des noms d'inconnues
  65. * primales, puis les duales associées (segment trav1)
  66. *
  67. C      WRITE(IOIMP,*) 'Coucou exdiar'
  68.       SEGACT MRIGID
  69.       nbincd=0
  70.       segini trav1
  71.       NRIGEL=COERIG(/1)
  72.       DO IRIG=1,NRIGEL
  73.          DESCR=IRIGEL(3,IRIG)
  74.          SEGACT DESCR
  75.          do 45 iligrp=1,lisinc(/2)
  76.             do 55 inc=1,nbincd
  77.                if (compd(inc).eq.lisinc(iligrp)) goto 45
  78.  55         continue
  79.             nbincd=nbincd+1
  80.             segadj trav1
  81.             compd(nbincd)=lisinc(iligrp)
  82.  45      continue
  83.          SEGDES DESCR
  84.       ENDDO
  85.       do iincd=1,nbincd
  86.          CALL PLACE(NOMDD,LNOMDD,idx,compd(iincd))
  87.          IF (idx.EQ.0) THEN
  88.             compdd(iincd)=compd(iincd)
  89.          ELSE
  90.             compdd(iincd)=NOMDU(idx)
  91.          ENDIF
  92.       ENDDO
  93. C         WRITE(IOIMP,*) '  1) compd'
  94. C         WRITE (IOIMP,2019) (compd(I),I=1,compd(/2))
  95. C         WRITE(IOIMP,*) '  2) compdd'
  96. C         WRITE (IOIMP,2019) (compdd(I),I=1,compdd(/2))
  97. *
  98. * Maintenant, on construit les CORES1 et CORES2 qui nous
  99. * permettent de savoir quels sont les termes diagonaux
  100. * dans les sous-matrices
  101. *
  102.       SEGINI CORES1
  103.       DO IRIG=1,NRIGEL
  104.          DESCR=IRIGEL(3,IRIG)
  105.          SEGACT DESCR
  106. *  Tout d'abord le segment trav2 qui permet d'avoir
  107. *  la correspondance sur les noms d'inconnues
  108.          NLIGRP=LISINC(/2)
  109.          NLIGRD=LISDUA(/2)
  110.          segini trav2
  111.          do 410 iligrp=1,nligrp
  112.             do 420 incd=1,nbincd
  113.                if (compd(incd).eq.lisinc(iligrp)) then
  114.                   cor1p(iligrp)=incd
  115.                   goto 410
  116.                endif
  117.  420        continue
  118.  410     continue
  119.          do 510 iligrd=1,nligrd
  120.             do 520 incd=1,nbincd
  121.                if (compdd(incd).eq.lisdua(iligrd)) then
  122.                   cor1d(iligrd)=incd
  123.                   goto 510
  124.                endif
  125.  520        continue
  126. * Si on a pas trouvé dans la liste des duales, on regarde
  127. * dans la liste des primales.
  128. * Ceci est pratique en mécaflu où il arrive qu'on nomme les duales
  129. * pareil que les primales (ex. 'UX' 'UX')
  130.             do 525 incd=1,nbincd
  131.                if (compd(incd).eq.lisdua(iligrd)) then
  132.                   cor1d(iligrd)=incd
  133.                   goto 510
  134.                endif
  135.  525        continue
  136.  510     continue
  137. *   Voyons si le descripteur est carré, ce qui est un
  138. *   cas courant qui simplifie notablement les choses
  139.          descar=.false.
  140.          if (nligrp.ne.nligrd) goto 600
  141.          do iligr=1,nligrp
  142.             if (cor1p(iligr).ne.cor1d(iligr)) goto 600
  143.             if (noelep(iligr).ne.noeled(iligr)) goto 600
  144.          enddo
  145.          descar=.true.
  146.  600     continue
  147. *   Si oui, tout est simple, sinon il faut rechercher
  148. *   une correspondance primale-duale
  149.          if (descar) then
  150.             ipcor2(irig)=-1
  151.          else
  152.             nligp=nligrp
  153.             segini cores2
  154.             nligp=0
  155.             do iligrp=1,nligrp
  156.                do iligrd=1,nligrd
  157.                   if (cor1p(iligrp).eq.cor1d(iligrd).and.
  158.      $                 noelep(iligrp).eq.noeled(iligrd)) then
  159.                      nligp=nligp+1
  160.                      cordia(1,nligp)=iligrp
  161.                      cordia(2,nligp)=iligrd
  162.                   endif
  163.                enddo
  164.             enddo
  165.             if (nligp.eq.0) then
  166.                segsup cores2
  167.                ipcor2(irig)=0
  168.             else
  169.                segadj cores2
  170.                ipcor2(irig)=cores2
  171.             endif
  172.          endif
  173.          itrav2(irig)=trav2
  174. C         Write(ioimp,*) 'irig=',irig
  175. C         Write(ioimp,*) 'descar=',descar
  176. C         if (ipcor2.gt.0) then
  177. C            write(ioimp,*) 'cordia(1, )'
  178. C         WRITE (IOIMP,2020) (cordia(1,I),I=1,cordia(/2))
  179. C            write(ioimp,*) 'cordia(2, )'
  180. C         WRITE (IOIMP,2020) (cordia(2,I),I=1,cordia(/2))
  181. C         else
  182. C            WRITE(IOIMP,*) 'ipcor2=',ipcor2
  183. C         endif
  184.          SEGDES DESCR
  185.       ENDDO
  186. *
  187. * Maintenant, on peut construire MTRA1, MTRA2
  188. * (repérage des couples noms d'inconnues primales,
  189. *  harmoniques) ainsi que CORES3
  190. * et ICPR (repérage global->local des noeuds)
  191. * ce qui permettra dans une deuxième passe de remplir
  192. * le segment TMTRAV qui sera transformé en CHPOINT
  193. *
  194.       SEGINI ICPR
  195.       NNNOE=0
  196.       SEGINI MTRA1,MTRA2
  197.       NNIN=0
  198.       DO 100 IRIG=1,NRIGEL
  199.          DESCR=IRIGEL(3,IRIG)
  200.          trav2=itrav2(irig)
  201.          MELEME=IRIGEL(1,IRIG)
  202.          cores2=ipcor2(irig)
  203.          if (cores2.eq.0) goto 99
  204.          SEGACT DESCR
  205.          SEGACT MELEME
  206.          if (cores2.eq.-1) then
  207.             nligp=noelep(/1)
  208.          else
  209.             nligp=cordia(/2)
  210.          endif
  211.          segini cores3
  212.          do 110 iligp=1,nligp
  213.             if (cores2.eq.-1) then
  214.                iligrp=iligp
  215.             else
  216.                iligrp=cordia(1,iligp)
  217.             endif
  218.             do ic=1,icomp(/1)
  219.                if (icomp(ic).eq.cor1p(iligrp).and.
  220.      $              mhar(ic).eq.irigel(5,irig)) then
  221.                   comtra(iligp)=ic
  222.                   goto 110
  223.                endif
  224.             enddo
  225.             icomp(**)=cor1p(iligrp)
  226.             mhar(**)=irigel(5,irig)
  227.             nnin=nnin+1
  228.             comtra(iligp)=nnin
  229.  110     continue
  230. C         WRITE(IOIMP,*) '  3) comtra'
  231. C         WRITE (IOIMP,2020) (comtra(I),I=1,comtra(/1))
  232.          ipcor3(irig)=cores3
  233. * Dans le cas cores2=-1
  234. *   tous les noeuds du meleme sont référencés
  235. *   sinon, on regarde ceux qui le sont
  236.          if (cores2.eq.-1) then
  237.             nno=num(/1)
  238.          else
  239.             nno=nligp
  240.          endif
  241.          do iel=1,num(/2)
  242.             do ino=1,nno
  243.                if (cores2.eq.-1) then
  244.                   ilno=ino
  245.                else
  246.                   ilno=noelep(cordia(1,ino))
  247.                endif
  248.                ipt=num(ilno,iel)
  249.                if (icpr(ipt).eq.0) then
  250.                   nnnoe=nnnoe+1
  251.                   icpr(ipt)=nnnoe
  252.                endif
  253.             enddo
  254.          enddo
  255.  99      continue
  256.          SEGDES MELEME
  257.          SEGDES DESCR
  258.          segsup trav2
  259.  100  CONTINUE
  260. C      WRITE(IOIMP,*) 'nnin=',nnin
  261. C      WRITE(IOIMP,*) 'nnin2=',icomp(/1)
  262. C      WRITE(IOIMP,*) 'nnin3=',mhar(/1)
  263. C      WRITE(IOIMP,*) '  1) icomp'
  264. C      WRITE (IOIMP,2020) (icomp(I),I=1,icomp(/1))
  265. C      WRITE(IOIMP,*) '  2) mhar'
  266. C      WRITE (IOIMP,2020) (mhar(I),I=1,mhar(/1))
  267. C      WRITE(IOIMP,*) 'nnnoe=',nnnoe
  268. C      WRITE(IOIMP,*) '  4) icpr'
  269. C      WRITE (IOIMP,2020) (icpr(I),I=1,icpr(/1))
  270. *
  271. * Maintenant, on peut remplir MTRAV
  272. *
  273.       SEGINI MTRAV
  274.       DO ININ=1,NNIN
  275.          INCO(ININ)=compd(ICOMP(ININ))
  276.          NHAR(ININ)=MHAR(ININ)
  277.       ENDDO
  278.       SEGSUP MTRA1,MTRA2
  279.       segsup trav1
  280. C      xt=0.D0
  281.       DO 1100 IRIG=1,NRIGEL
  282.          MELEME=IRIGEL(1,IRIG)
  283.          DESCR=IRIGEL(3,IRIG)
  284.          XMATRI=IRIGEL(4,IRIG)
  285.          cores2=ipcor2(irig)
  286.          cores3=ipcor3(irig)
  287.          if (cores2.eq.0) goto 1099
  288.          SEGACT MELEME
  289.          SEGACT DESCR
  290.          SEGACT XMATRI
  291.          nel=num(/2)
  292.          if (cores2.eq.-1) then
  293.             nligp=noelep(/1)
  294.          else
  295.             nligp=cordia(/2)
  296.          endif
  297.          do iel=1,nel
  298.             do iligp=1,nligp
  299.                if (cores2.eq.-1) then
  300.                   iligrp=iligp
  301.                   iligrd=iligrp
  302.                else
  303.                   iligrp=cordia(1,iligp)
  304.                   iligrd=cordia(2,iligp)
  305.                endif
  306.                ININ=comtra(iligp)
  307.                ino=noelep(iligrp)
  308.                ipt=num(ino,iel)
  309.                INNOE=icpr(ipt)
  310. C               if (ipt.eq.2.and.lisinc(iligrp).eq.'UX  ') then
  311. C                  write(ioimp,*) 'irig=',irig
  312. C                  write(ioimp,*) 'iel=',iel
  313. C                  write(ioimp,*) 'iligrp=',iligrp
  314. C                  write(ioimp,*) 'iligrd=',iligrd
  315. C                  xt=xt+RE(iligrd,iligrp,iel)
  316. C                  WRITE(IOIMP,*) 'RE=',RE(iligrd,iligrp,iel)
  317. C                  WRITE(IOIMP,*) 'xt=',xt
  318. C                  WRITE(IOIMP,*) 'inin=',inin
  319. C                  WRITE(IOIMP,*) 'innoe=',innoe
  320. C               endif
  321.                IGEO(INNOE)=ipt
  322.                IBIN(ININ,INNOE)=1
  323.                BB(ININ,INNOE)=BB(ININ,INNOE)+
  324.      $              RE(iligrd,iligrp,iel)*COERIG(IRIG)
  325.             enddo
  326.          enddo
  327.          if (cores2.ne.-1) segsup cores2
  328.          segsup cores3
  329.          SEGDES XMATRI
  330.          SEGDES DESCR
  331.          SEGDES MELEME
  332.  1099    continue
  333.  1100 CONTINUE
  334.       SEGSUP ICPR
  335.       SEGSUP CORES1
  336.       SEGDES MRIGID
  337. *
  338. * Création du CHPOINT résultat
  339. *
  340.       CALL CRECHP(MTRAV,MCHPOI)
  341.      $
  342.       SEGSUP MTRAV
  343.       RETURN
  344.  2019 FORMAT (20(2X,A4) )
  345.  2020 FORMAT (20(2X,I4) )
  346.  2021 FORMAT (20(2X,L4) )
  347.  2022 FORMAT(10(1X,1PG12.5))
  348. *
  349. * End of EXDIAG
  350. *
  351.       END
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  

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