Télécharger la plaquette Cast3M


fr en

Télécharger opto3.eso

Retour à la liste

Numérotation des lignes : 

opto3

  1. C OPTO3     SOURCE    GOUNAND   26/06/09    21:15:11     12566          
  2.       SUBROUTINE OPTO3(TRAVJ,TRAVX,JELEM,TRAVK,TRAVL,ICPRX,IDCPX,KELEMX
  3.      $     ,JTBES,JCAND)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5.       IMPLICIT INTEGER (I-N)
  6. C***********************************************************************
  7. C NOM         : OPTO3 (anciennement optt3d)
  8. C DESCRIPTION : Une implémentation de l'amélioration d'une topologie
  9. C               autour d'un élément. On reprend OPTITOPO pour le corps
  10. C     du programme. On reprend l'extraction et la topologie inverse de
  11. C     EXTO. Le point crucial sera d'implémenter la modification de la
  12. C     topologie : enlever les anciens éléments et mettre les nouveaux.
  13. C
  14. C
  15. C     Ici, on passe dans une deuxième numérotation locale à JEXTO avant
  16. C     d'interfacer à TOPV2
  17. C
  18. C     On ressemble fortement à opto1.eso
  19. C     optt3 -> optt3d : TRAVX devient un segment TRAVJ
  20. C     optt3c -> optt3d : on interface à TOPV2 au lieu de TOPVOL
  21. C          du coup on ressemble à topv1
  22. C
  23. C LANGAGE     : ESOPE
  24. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SEMT/LTA)
  25. C               mél : gounand@semt2.smts.cea.fr
  26. C***********************************************************************
  27. C APPELES          : TOPV2
  28. C APPELES (E/S)    :
  29. C APPELES (BLAS)   :
  30. C APPELES (CALCUL) :
  31. C APPELE PAR       : OPTO2
  32. C***********************************************************************
  33. C SYNTAXE GIBIANE    :
  34. C ENTREES            : JCOORD, TRAVX, JELEM
  35. C ENTREES/SORTIES    :
  36. C SORTIES            : JTBES,JCAND
  37. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  38. C***********************************************************************
  39. C VERSION    : v1, 11/11/2017, version initiale
  40. C HISTORIQUE : v1, 11/11/2017, création
  41. C HISTORIQUE :
  42. C HISTORIQUE :
  43. C***********************************************************************
  44. -INC PPARAM
  45. -INC CCOPTIO
  46. -INC TMATOP2
  47. -INC SMELEME
  48. *
  49. *     Le nombre d'éléments de JTOPO et le nombre de points de JCOORD
  50. *     vont être variables. Pour ne pas avoir à ajuster ces segments en
  51. *     permanence, on va dimensionner plus large, mais du coup, il faut
  52. *     aussi maintenir à la main le nombre de noeuds et d'éléments
  53. *     courants.
  54. *
  55. *     Le nombre d'éléments courants est NVCOU et le nombre d'éléments
  56. *     max est NVMAX. Idem pour le nombre de noeuds courants et max :
  57. *     NPCOU et NPMAX.
  58. *
  59.       POINTEUR JELEM.MELEME
  60.       POINTEUR JEXTO.MELEME,KEXTO.MELEME
  61.       POINTEUR JTBES.MELEME,KTBES.MELEME
  62.       POINTEUR JPVIRT.MELEME
  63. -INC SMCOORD
  64. * Numerotation locale (de 1 à NBPTS)
  65.       POINTEUR JCOORD.MCOORD
  66.       POINTEUR KCOORD.MCOORD
  67. -INC TMATOP1
  68. *-INC STOPINV
  69. *-INC SMETRIQ
  70.       POINTEUR JCMETR.METRIQ
  71.       POINTEUR KCMETR.METRIQ
  72. *-INC STRAVJ
  73.       POINTEUR TRAVX.TRAVJ
  74.       POINTEUR TRAVK.TRAVJ
  75. *-INC STRAVL
  76. * Passage de numerotation globale -> locale
  77. *   et locale -> globale
  78.       POINTEUR ICPRX.MLENTX
  79.       POINTEUR IDCPX.MLENTX
  80.       POINTEUR KELEMX.MELEMX
  81. *
  82.       logical lchang,lchtop, ldbg
  83.       CHARACTER*40 CFORM
  84. *
  85. * Executable statements
  86. *
  87. *      if (impr.ge.3) WRITE(IOIMP,*) 'Entrée dans opto3.eso'
  88. *
  89. * Initialisation et extension des segments JTOPO et JCOORD
  90. *
  91.       IDIMP=IDIM+1
  92. *
  93.       JCOORD=TRAVJ.COORD
  94.       JCMETR=TRAVJ.CMETR
  95.       JEXTO=TRAVX.TOPO
  96.       KEXTO=TRAVK.TOPO
  97.       KCOORD=TRAVK.COORD
  98.       KCMETR=TRAVK.CMETR
  99. * Correspondances de numérotation
  100.       JGDONN=XCOOR(/1)/(IDIM+1)
  101.       CALL mtxadj(ICPRX,JGDONN,lchang,'opto3 : ICPRX_dim')
  102.       if (ierr.ne.0) return
  103. *     Ici, on met le IDCP à la même dimension pour simplifier le
  104. *     remplissage (fait en une passe)
  105.       CALL mtxadj(IDCPX,JGDONN,lchang,'opto3 : IDCPX_dim1')
  106.       if (ierr.ne.0) return
  107. *     SEGINI ICPR
  108.  
  109.       JPVIRT=TRAVJ.PVIRT
  110.       NJPVIR=JPVIRT.NUM(/2)
  111. * On va d'abord noter les noeuds virtuels
  112.       IK=0
  113.       DO 13 IEL=1,TRAVX.NVCOU
  114.          DO 130 INO=1,JEXTO.NUM(/1)
  115.             IP=JEXTO.NUM(INO,IEL)
  116.             IF (IP.LE.NJPVIR) THEN
  117.                IF (ICPRX.LECTX(IP).EQ.0) THEN
  118.                   IK=IK+1
  119.                   ICPRX.LECTX(IP)=IK
  120.                   IDCPX.LECTX(IK)=IP
  121.                ENDIF
  122.             ENDIF
  123.  130     CONTINUE
  124.  13   CONTINUE
  125.       NKPVIR=IK
  126. * Puis les autres noeuds
  127.       DO 23 IEL=1,TRAVX.NVCOU
  128.          DO 230 INO=1,JEXTO.NUM(/1)
  129.             IP=JEXTO.NUM(INO,IEL)
  130.             IF (IP.GT.NJPVIR) THEN
  131.                IF (ICPRX.LECTX(IP).EQ.0) THEN
  132.                   IK=IK+1
  133.                   ICPRX.LECTX(IP)=IK
  134.                   IDCPX.LECTX(IK)=IP
  135.                ENDIF
  136.             ENDIF
  137.  230     CONTINUE
  138.  23   CONTINUE
  139. *
  140.       NPTINI=IK
  141.       CALL mtxadj(IDCPX,NPTINI,lchang,'opto3 : IDCPX_dim2')
  142.       if (ierr.ne.0) return
  143.       NPTBAS=TRAVJ.NPCOU
  144.       if (impr.ge.4) then
  145. *         write(ioimp,*) 'NJPVIR=',NJPVIR,'NKPVIR=',NKPVIR
  146.          write(ioimp,'(2(A,2I3,1X))')
  147.      $        'opto3 : Nb noeuds locaux, tlocaux=',NPTBAS,IK
  148.      $        ,'virtuel, tvirtuels',NJPVIR,NKPVIR
  149.          if (impr.ge.6) then
  150.             write(ioimp,*) 'IDCPX'
  151.             write(ioimp,187) (IDCPX.LECTX(I),I=1,IDCPX.JGCOU)
  152.          endif
  153.       endif
  154.  125  format ('(A,2(1X,I3),20(2X,''[''',I1,'I3,'']''))')
  155.       IF (IMPR.GE.4) THEN
  156.          nnod=jexto.num(/1)
  157.          nel=travx.nvcou
  158.          write(cform,125) nnod
  159.          write(ioimp,cform)
  160.      $        'opto3 : topologie ext. en coord locales  nnod,nel=',nnod
  161.      $        ,nel,((jexto.num(ino,iel),ino=1,nnod),iel=1,nel)
  162. *         write(ioimp,*)
  163. *     $        'opto3.eso : topologie ext. en coord locales : '
  164. *         call ecmai1(jexto,0)
  165. *         segact jexto*mod
  166.       ENDIF
  167. *
  168. * Melemes en coordonnées tlocales
  169. *
  170.       NBLEXT=TRAVX.NVCOU
  171.       CALL TOPADV(TRAVK,NBLEXT,1,lchang,'opto3 : TRAVK')
  172.       IF (IERR.NE.0) RETURN
  173. * KEXTO en nouvelle numérotation
  174.       DO 33 IEL=1,travk.nvcou
  175. *anc         DO 330 INO=1,KEXTO.NUM(/1)
  176.          DO 330 INO=1,IDIMP
  177.             IP=JEXTO.NUM(INO,IEL)
  178. *            JP=ICPR(IP)
  179.             JP=ICPRX.LECTX(IP)
  180.             IF (JP.NE.0) THEN
  181.                KEXTO.NUM(INO,IEL)=JP
  182.             ELSE
  183.                WRITE(IOIMP,*) 'Erreur de programmation'
  184.                GOTO 9999
  185.             ENDIF
  186.  330     CONTINUE
  187.  33   CONTINUE
  188. *tmp      if (iveri.ge.2.and.lchang) then
  189. *tmp         call vetopi(travk,'opto3 : Apres extension elem travk')
  190. *tmp         if (ierr.ne.0) return
  191. *tmp      endif
  192.       IF (IMPR.GE.4) THEN
  193.          nnod=kexto.num(/1)
  194.          nel=travk.nvcou
  195.          write(cform,125) nnod
  196.          write(ioimp,cform)
  197.      $        'opto3 : topologie ext. en coord tlocales nnod,nel=',nnod
  198.      $        ,nel,((kexto.num(ino,iel),ino=1,nnod),iel=1,nel)
  199. *         write(ioimp,*)
  200. *     $        'opto3.eso : topologie ext. en coord tlocales : '
  201. *         call ecmai1(kexto,0)
  202. *         segact kexto*mod
  203.       ENDIF
  204. *     Normalement plus besoin de transformer les noeuds virtuels :
  205. *     Ce sont juste les NKPVIR premiers
  206.       TRAVK.PVIRT=NKPVIR
  207. *      IF (IMPR.GE.4) THEN
  208. *         write(ioimp,*)
  209. *     $        'opto3.eso : noeuds virtuels en coord tlocales : ',NKPVIR
  210. *      ENDIF
  211. * Vérifier que JELEM n'a qu'un élément et en supprimer les noeuds nuls ?
  212.       NBELEM=JELEM.NUM(/2)
  213.       IF (NBELEM.NE.1) THEN
  214.          write(ioimp,*) 'on veut que jelem n''ait qu''un element'
  215.          goto 9999
  216.       ENDIF
  217.       nbnn=0
  218.       do ino=1,jelem.num(/1)
  219.          ip=jelem.NUM(ino,1)
  220.          if (ip.ne.0) then
  221.             IF (IP.LT.0.OR.IP.GT.NPTBAS) THEN
  222.                write(ioimp,*) 'ip=',ip,' noeud hors bornes dans jelem'
  223.                goto 9999
  224.             ELSE
  225.                JP=ICPRX.LECTX(IP)
  226.                IF (JP.EQ.0) then
  227. *                  WRITE(IOIMP,*)
  228. *     $                 'La topologie extraite devrait contenir',
  229. *     $                 'le noeud de ielem ip,jp=',ip,jp
  230. *     goto 9999
  231.                ELSE
  232.                   nbnn=nbnn+1
  233.                   kelemx.numx(nbnn,1)=jp
  234.                ENDIF
  235.             ENDIF
  236.          endif
  237.       enddo
  238.       kelemx.nncou=nbnn
  239.  
  240. *
  241.       IF (IMPR.GE.4) THEN
  242.          nnod=kelemx.numx(/1)
  243.          nel=1
  244.          write(cform,125) nnod
  245.          write(ioimp,cform)
  246.      $        'opto3 : element en coord tlocales nnod,nel=',nnod,nel
  247.      $        ,((kelemx.numx(ino,iel),ino=1,nnod),iel=1,nel)
  248. *         write(ioimp,*) 'opto3.eso : element en coord tlocales : '
  249. *         call ecmelx(kelemx,0)
  250.       ENDIF
  251. * Passage des coordonnées en tlocale
  252. *     NBPTS=NPTINI
  253. *      NBPTS=TRAVK.NPMAX
  254. *      SEGINI,KCOORD
  255. *     TRAVK.COORD=KCOORD
  256.  
  257.       CALL TOPADP(TRAVK,NPTINI,1,lchang,'opto3 : TRAVK')
  258.       IF (IERR.NE.0) RETURN
  259.  
  260.       DO 53 IPL=1,TRAVK.NPCOU
  261.          IREFL=IDIMP*(IPL-1)
  262. *         IP=IDCP(IPL)
  263.          IP=IDCPX.LECTX(IPL)
  264.          IREF=IDIMP*(IP-1)
  265.          DO 530 IC=1,IDIMP
  266.             KCOORD.XCOOR(IREFL+IC)=JCOORD.XCOOR(IREF+IC)
  267.  530     CONTINUE
  268.          IF (JCMETR.NE.0) THEN
  269.             DO 540 ININ=1,JCMETR.XIN(/1)
  270.                KCMETR.XIN(ININ,IPL)=JCMETR.XIN(ININ,IP)
  271.  540        CONTINUE
  272.          ENDIF
  273.  53   CONTINUE
  274. *tmp      if (iveri.ge.2.and.lchang) then
  275. *tmp         call vetopi(travk,'opto3 : Apres extension noeud travk')
  276. *tmp         if (ierr.ne.0) return
  277. *tmp      endif
  278.  
  279. *tst      WRITE(IOIMP,185) 'SEGMENT KCOORD ',KCOORD
  280. *tst         WRITE(FORMA,FMT='("(1(",I1,"(1PG12.5,2X)))")') IDIMP
  281. *tst         write(ioimp,*) 'forma=',forma
  282. *tst         write(ioimp,*) 'XCOOR'
  283. *tst         write(ioimp,forma)  (kcoord.xcoor(I),I=1,kcoord.xcoor(/1))
  284. *
  285. *     La numérotation globale devient la locale dans ce bloc  !!!
  286.       MCOORD=KCOORD
  287.       CALL TOPV2(TRAVK,KELEMX,IALGO,IAJNO,XVTOL,QTOL,IMET,IMOMET,XDENS
  288.      $     ,INCMA,ISTMA,KTBES,JCAND,JNASCM,iveri,impr,TRAVL,LCHTOP
  289.      $     ,jcritq,pcritq,qcritq)
  290. *     write(ioimp,*) 'jcand=',jcand
  291. *     Point de branchement si erreur pendant le bloc en numérotation
  292. *     locale
  293.  555  CONTINUE
  294. *     On rétablit la numérotation globale originelle et on rajoute les
  295. *     noeuds nouvellement créés
  296.       MCOORD=JCOORD
  297. *     On part en erreur après le rétablissement du MCOORD global
  298.       IF (IERR.NE.0) RETURN
  299. *!!!changé      NPTFIN=KCOORD.XCOOR(/1)/IDIMP
  300.       NPTFIN=travk.npcou
  301. *      write(ioimp,*) 'NPTINI,NPTFIN=',NPTINI,NPTFIN
  302.       IF (NPTINI.NE.NPTFIN) THEN
  303.          if (impr.ge.4) write(ioimp,'(A,I5,1X,A)') 'opto3 :',nptfin
  304.      $        -nptini,' nouveaux noeuds crees'
  305. *         if (ktbes.eq.kexto) then
  306.          if (.not.lchtop) then
  307.             write(ioimp,*) 'kexto non ameliore mais...'
  308.             write(ioimp,*) nptfin-nptini,' nouveaux noeuds crees'
  309.             write(ioimp,*) 'pas logique...'
  310.             goto 9999
  311.          endif
  312.          NBPTA=NPTBAS
  313.          NPCOUN=NBPTA+NPTFIN-NPTINI
  314. *
  315.          CALL TOPADP(TRAVJ,NPCOUN,1,lchang,'opto3 : TRAVJ')
  316.          if (ierr.ne.0) return
  317. *
  318.          DO 5000 I=NPTINI+1,NPTFIN
  319.             DO 5010 J=1,IDIMP
  320.                JCOORD.XCOOR(NBPTA*IDIMP+J)=KCOORD.XCOOR((I-1)*IDIMP
  321.      $              +J)
  322.  5010       CONTINUE
  323.             IF (JCMETR.NE.0) THEN
  324.                DO 5020 ININ=1,JCMETR.XIN(/1)
  325.                   JCMETR.XIN(ININ,NBPTA+1)=KCMETR.XIN(ININ,I)
  326.  5020          CONTINUE
  327.             ENDIF
  328.             NBPTA=NBPTA+1
  329.  5000    CONTINUE
  330.          if (iveri.ge.2.and.lchang) then
  331.             call vetopi(travj,'opto3 : Apres extension travj')
  332.             if (ierr.ne.0) return
  333.          endif
  334.       ENDIF
  335. *
  336. *      IF (KTBES.EQ.KEXTO) THEN
  337.       IF (.not.lchtop) THEN
  338.          JTBES=JEXTO
  339.       ELSE
  340. * En place
  341.          JTBES=KTBES
  342.          SEGACT JTBES*MOD
  343.          DO 63 IEL=1,JTBES.NUM(/2)
  344.             DO 630 INO=1,JTBES.NUM(/1)
  345.                IPL=JTBES.NUM(INO,IEL)
  346.                IF (IPL.LE.NPTINI) THEN
  347. *                  IP=IDCP(IPL)
  348.                   IP=IDCPX.LECTX(IPL)
  349.                ELSE
  350.                   IP=IPL-NPTINI+NPTBAS
  351.                ENDIF
  352.                JTBES.NUM(INO,IEL)=IP
  353.  630        CONTINUE
  354.  63      CONTINUE
  355.          IF (IMPR.GE.4) THEN
  356.             nnod=jtbes.num(/1)
  357.             nel=jtbes.num(/2)
  358.             write(cform,125) nnod
  359.             write(ioimp,cform)
  360.             write(ioimp,*) 'opto3 : topologie amelioree nnod,nel=',nnod
  361.      $        ,nel,((jtbes.num(ino,iel),ino=1,nnod),iel=1,nel)
  362. *            call ecmai1(jtbes,0)
  363. *            segact jtbes*mod
  364.          ENDIF
  365.       ENDIF
  366. *     SEGSUP KELEM
  367.       kelemx.nncou=0
  368. *
  369. * RAZ IDCPX et ICPRX
  370. *     SEGSUP IDCP
  371. *     SEGSUP ICPR
  372.       DO 1500 I=1,IDCPX.JGCOU
  373.          ICPRX.LECTX(IDCPX.LECTX(I))=0
  374.          IDCPX.LECTX(I)=0
  375.  1500 CONTINUE
  376.       CALL mtxadj(ICPRX,0,lchang,'opto3 : ICPRX_sup')
  377.       if (ierr.ne.0) return
  378.       CALL mtxadj(IDCPX,0,lchang,'opto3 : IDCPX_sup')
  379.       if (ierr.ne.0) return
  380. *
  381. * Normal termination
  382. *
  383.       RETURN
  384. *
  385. * Format handling
  386. *
  387.  187  FORMAT (5X,10I8)
  388.  188  FORMAT ('Apres point IELEM(',I2,',1)=',I6,' ; NBL=')
  389.  189  FORMAT ('Le noeud ',I2,'/',I2,' de IELEM de numero',I6
  390.      $     ,' a le plus petit nb de voisins :',I3)
  391. *
  392. * Error handling
  393. *
  394.  9999 CONTINUE
  395.       MOTERR(1:8)='OPTO3  '
  396. * 349 2
  397. *Problème non prévu dans le s.p. %m1:8 contactez votre assistance
  398.       CALL ERREUR(349)
  399.       RETURN
  400. *
  401. * End of subroutine OPTO3
  402. *
  403.       END
  404.  
  405.  

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