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   21/04/06    21:15:19     10940          
  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 KELEM.MELEME
  61.       POINTEUR JEXTO.MELEME,KEXTO.MELEME
  62.       POINTEUR JTBES.MELEME,KTBES.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. *
  77.  
  78. *
  79. * Passage de numerotation globale -> locale
  80. *   et locale -> globale    
  81. *      SEGMENT ICPR(XCOOR(/1)/(IDIM+1))
  82. *      SEGMENT IDCP(NPTINI)
  83. *-INC SMLENTX
  84.       POINTEUR ICPRX.MLENTX
  85.       POINTEUR IDCPX.MLENTX
  86. *-INC SMELEMX
  87.       POINTEUR KELEMX.MELEMX
  88.  
  89.       logical lchang,lchtop
  90. *     
  91. * Executable statements
  92. *
  93.       if (impr.ge.3) WRITE(IOIMP,*) 'Entrée dans opto3.eso'
  94. *
  95. * Initialisation et extension des segments JTOPO et JCOORD
  96. *
  97.       IDIMP=IDIM+1
  98. *
  99.       JCOORD=TRAVJ.COORD
  100.       JCMETR=TRAVJ.CMETR
  101.       JEXTO=TRAVX.TOPO
  102.       KEXTO=TRAVK.TOPO
  103.       KCOORD=TRAVK.COORD
  104.       KCMETR=TRAVK.CMETR
  105. * Correspondances de numérotation
  106.       JGDONN=XCOOR(/1)/(IDIM+1)
  107.       CALL mtxadj(ICPRX,JGDONN,lchang,'opto3 : ICPRX_dim')
  108.       if (ierr.ne.0) return
  109. *     Ici, on met le IDCP à la même dimension pour simplifier le
  110. *     remplissage (fait en une passe)
  111.       CALL mtxadj(IDCPX,JGDONN,lchang,'opto3 : IDCPX_dim1')
  112.       if (ierr.ne.0) return
  113. *      SEGINI ICPR
  114.       IK=0
  115.       DO 23 IEL=1,TRAVX.NVCOU
  116.          DO 230 INO=1,JEXTO.NUM(/1)
  117.             IP=JEXTO.NUM(INO,IEL)
  118.             IF (ICPRX.LECTX(IP).EQ.0) THEN
  119.                IK=IK+1
  120. *               ICPR(IP)=IK
  121.                ICPRX.LECTX(IP)=IK
  122.                IDCPX.LECTX(IK)=IP
  123.             ENDIF
  124.  230     CONTINUE
  125.  23   CONTINUE
  126. *
  127.       NPTINI=IK
  128.       CALL mtxadj(IDCPX,NPTINI,lchang,'opto3 : IDCPX_dim2')
  129.       if (ierr.ne.0) return
  130. *      SEGINI IDCP
  131.       NPTBAS=TRAVJ.NPCOU
  132. **      write(ioimp,*) 'mcoord,jcoord,npt,npcou,ik=',mcoord,jcoord,xcoor(
  133. **     $     /1)/idimp,npcou,ik
  134. *      DO 500 I=1,NPTBAS
  135. **         if (icpr(i).ne.0) IDCP(ICPR(I))=I
  136. *         if (icprx.lectx(i).ne.0) IDCPx.lectx(ICPRx.lectx(I))=I
  137. * 500  CONTINUE
  138.       if (impr.ge.4) then 
  139.          write(ioimp,*) 'Nb noeud locaux, tlocaux=',NPTBAS,IK
  140. *         write(ioimp,*) 'ICPR'
  141. *         write(ioimp,187) (ICPR(I),I=1,ICPR(/1))
  142. *         write(ioimp,*) 'IDCP'
  143. *         write(ioimp,187) (IDCP(I),I=1,IDCP(/1))
  144.          write(ioimp,*) 'IDCPX'
  145.          write(ioimp,187) (IDCPX.LECTX(I),I=1,IDCPX.JGCOU)
  146.       endif
  147.       IF (IMPR.GE.4) THEN
  148.          write(ioimp,*)
  149.      $        'opto3.eso : topologie ext. en coord locales : '
  150.          call ecmai1(jexto,0)
  151.          segact jexto*mod
  152.       ENDIF
  153. *
  154. * Melemes en coordonnées tlocales
  155. *
  156.       NBLEXT=TRAVX.NVCOU
  157.       CALL TOPADV(TRAVK,NBLEXT,1,lchang,'opto3 : TRAVK')
  158.       IF (IERR.NE.0) RETURN
  159. * KEXTO en nouvelle numérotation      
  160.       DO 33 IEL=1,travk.nvcou
  161. *anc         DO 330 INO=1,KEXTO.NUM(/1)
  162.          DO 330 INO=1,IDIMP
  163.             IP=JEXTO.NUM(INO,IEL)
  164. *            JP=ICPR(IP)
  165.             JP=ICPRX.LECTX(IP)
  166.             IF (JP.NE.0) THEN
  167.                KEXTO.NUM(INO,IEL)=JP
  168.             ELSE
  169.                WRITE(IOIMP,*) 'Erreur de programmation'
  170.                GOTO 9999
  171.             ENDIF
  172.  330     CONTINUE
  173.  33   CONTINUE
  174. *tmp      if (iveri.ge.2.and.lchang) then
  175. *tmp         call vetopi(travk,'opto3 : Apres extension elem travk')
  176. *tmp         if (ierr.ne.0) return
  177. *tmp      endif            
  178.       IF (IMPR.GE.4) THEN
  179.          write(ioimp,*)
  180.      $        'opto3.eso : topologie ext. en coord tlocales : '
  181.          call ecmai1(kexto,0)
  182.          segact kexto*mod
  183.       ENDIF
  184. *     Noeud virtuel en coordonnées locales
  185.       JPVIRT=TRAVJ.PVIRT
  186.       IF (JPVIRT.NE.0) THEN
  187. *         KPVIRT=ICPR(JPVIRT)
  188.          KPVIRT=ICPRX.LECTX(JPVIRT)
  189. * Ici, c'est normal de ne pas être inclus         
  190. *         IF (KPVIRT.EQ.0) THEN
  191. *            write(ioimp,*)
  192. *     $           'Noeud virtuel non inclus dans la topologie ?'
  193. *            goto 9999
  194. *         ENDIF
  195.       ELSE
  196.          KPVIRT=0
  197.       ENDIF
  198.       TRAVK.PVIRT=KPVIRT
  199.       IF (IMPR.GE.4) THEN
  200.          write(ioimp,*)
  201.      $        'opto3.eso : noeud virtuel en coord tlocales : ',KPVIRT
  202.       ENDIF
  203. * Vérifier que JELEM n'a qu'un élément et en supprimer les noeuds nuls ?
  204.       NBELEM=JELEM.NUM(/2)
  205.       IF (NBELEM.NE.1) THEN
  206.          write(ioimp,*) 'on veut que jelem n''ait qu''un element'
  207.          goto 9999
  208.       ENDIF
  209.       nbnn=0
  210.       do ino=1,jelem.num(/1)
  211.          ip=jelem.NUM(ino,1)
  212.          if (ip.ne.0) then
  213.             IF (IP.LT.0.OR.IP.GT.NPTBAS) THEN
  214.                write(ioimp,*) 'ip=',ip,' noeud hors bornes dans jelem'
  215.                goto 9999
  216.             ELSE
  217. *               JP=ICPR(IP)
  218.                JP=ICPRX.LECTX(IP)
  219.                IF (JP.EQ.0) then
  220. *                  WRITE(IOIMP,*)
  221. *     $                 'La topologie extraite devrait contenir',
  222. *     $                 'le noeud de ielem ip,jp=',ip,jp
  223. *     goto 9999
  224.                ELSE   
  225.                   nbnn=nbnn+1
  226.                   kelemx.numx(nbnn,1)=jp
  227.                ENDIF
  228.             ENDIF
  229.          endif
  230.       enddo
  231.       kelemx.nncou=nbnn
  232. *      nbsous=0
  233. *      nbref=0
  234. *     Element autour duquel on extrait
  235. *      SEGINI KELEM
  236. *      ibnn=0
  237. *      DO INO=1,JELEM.NUM(/1)
  238. *         IP=JELEM.NUM(INO,1)
  239. *         IF (IP.NE.0) THEN
  240. **            JP=ICPR(IP)
  241. *            JP=ICPRX.LECTX(IP)
  242. *            IF (JP.NE.0) then
  243. *               ibnn=ibnn+1
  244. *               KELEM.NUM(ibnn,1)=jp
  245. *            ENDIF
  246. *         ENDIF
  247. *      ENDDO
  248. *                 
  249.       IF (IMPR.GE.4) THEN
  250.          write(ioimp,*) 'opto3.eso : element en coord tlocales : '
  251. *         call ecmai1(kelem,0)
  252. *         segact kelem*mod
  253.          call ecmelx(kelemx,0)
  254.       ENDIF         
  255. * Passage des coordonnées en tlocale
  256. *     NBPTS=NPTINI
  257. *      NBPTS=TRAVK.NPMAX
  258. *      SEGINI,KCOORD
  259. *     TRAVK.COORD=KCOORD
  260.  
  261.       CALL TOPADP(TRAVK,NPTINI,1,lchang,'opto3 : TRAVK')
  262.       IF (IERR.NE.0) RETURN
  263.  
  264.       DO 53 IPL=1,TRAVK.NPCOU
  265.          IREFL=IDIMP*(IPL-1)
  266. *         IP=IDCP(IPL)
  267.          IP=IDCPX.LECTX(IPL)
  268.          IREF=IDIMP*(IP-1)
  269.          DO 530 IC=1,IDIMP
  270.             KCOORD.XCOOR(IREFL+IC)=JCOORD.XCOOR(IREF+IC)
  271.  530     CONTINUE
  272.          IF (JCMETR.NE.0) THEN
  273.             DO 540 ININ=1,JCMETR.XIN(/1)
  274.                KCMETR.XIN(ININ,IPL)=JCMETR.XIN(ININ,IP)
  275.  540        CONTINUE
  276.          ENDIF
  277.  53   CONTINUE
  278. *tmp      if (iveri.ge.2.and.lchang) then
  279. *tmp         call vetopi(travk,'opto3 : Apres extension noeud travk')
  280. *tmp         if (ierr.ne.0) return
  281. *tmp      endif            
  282.  
  283. *tst      WRITE(IOIMP,185) 'SEGMENT KCOORD ',KCOORD
  284. *tst         WRITE(FORMA,FMT='("(1(",I1,"(1PG12.5,2X)))")') IDIMP
  285. *tst         write(ioimp,*) 'forma=',forma
  286. *tst         write(ioimp,*) 'XCOOR'
  287. *tst         write(ioimp,forma)  (kcoord.xcoor(I),I=1,kcoord.xcoor(/1))
  288. *
  289. *     La numérotation globale devient la locale dans ce bloc  !!!
  290.       MCOORD=KCOORD
  291. *ijob      CALL TOPV2(TRAVK,KELEM,IJOB,XVTOL,QTOL,IMET,XDENS,
  292. *      CALL TOPV2(TRAVK,KELEM,IALGO,IAJNO,XVTOL,QTOL,IMET,XDENS,INCMA
  293.       CALL TOPV2(TRAVK,KELEMX,IALGO,IAJNO,XVTOL,QTOL,IMET,IMOMET,XDENS
  294.      $     ,INCMA,ISTMA,KTBES,JCAND,JNASCM,iveri,impr,TRAVL,LCHTOP)
  295. *     write(ioimp,*) 'jcand=',jcand
  296. *     Point de branchement si erreur pendant le bloc en numérotation
  297. *     locale
  298.  555  CONTINUE
  299. *     On rétablit la numérotation globale originelle et on rajoute les
  300. *     noeuds nouvellement créés
  301.       MCOORD=JCOORD
  302. *     On part en erreur après le rétablissement du MCOORD global
  303.       IF (IERR.NE.0) RETURN
  304. *!!!changé      NPTFIN=KCOORD.XCOOR(/1)/IDIMP
  305.       NPTFIN=travk.npcou
  306.       IF (NPTINI.NE.NPTFIN) THEN
  307.          if (impr.ge.4)
  308.      $        write(ioimp,*) nptfin-nptini,' nouveaux noeuds crees'
  309. *         if (ktbes.eq.kexto) then
  310.          if (.not.lchtop) then
  311.             write(ioimp,*) 'kexto non ameliore mais...'
  312.             write(ioimp,*) nptfin-nptini,' nouveaux noeuds crees'
  313.             write(ioimp,*) 'pas logique...'
  314.             goto 9999
  315.          endif
  316.          NBPTA=NPTBAS
  317.          NPCOUN=NBPTA+NPTFIN-NPTINI
  318. *
  319.          CALL TOPADP(TRAVJ,NPCOUN,1,lchang,'opto3 : TRAVJ')
  320.          if (ierr.ne.0) return
  321. *
  322.          DO 5000 I=NPTINI+1,NPTFIN
  323.             DO 5010 J=1,IDIMP
  324.                JCOORD.XCOOR(NBPTA*IDIMP+J)=KCOORD.XCOOR((I-1)*IDIMP
  325.      $              +J)
  326.  5010       CONTINUE
  327.             IF (JCMETR.NE.0) THEN
  328.                DO 5020 ININ=1,JCMETR.XIN(/1)
  329.                   JCMETR.XIN(ININ,NBPTA+1)=KCMETR.XIN(ININ,I)
  330.  5020          CONTINUE
  331.             ENDIF
  332.             NBPTA=NBPTA+1
  333.  5000    CONTINUE
  334.          if (iveri.ge.2.and.lchang) then
  335.             call vetopi(travj,'opto3 : Apres extension travj')
  336.             if (ierr.ne.0) return
  337.          endif            
  338.       ENDIF
  339. *      
  340. *      IF (KTBES.EQ.KEXTO) THEN
  341.       IF (.not.lchtop) THEN
  342.          JTBES=JEXTO
  343.       ELSE
  344. * En place
  345.          JTBES=KTBES
  346.          SEGACT JTBES*MOD
  347.          DO 63 IEL=1,JTBES.NUM(/2)
  348.             DO 630 INO=1,JTBES.NUM(/1)
  349.                IPL=JTBES.NUM(INO,IEL)
  350.                IF (IPL.LE.NPTINI) THEN
  351. *                  IP=IDCP(IPL)
  352.                   IP=IDCPX.LECTX(IPL)
  353.                ELSE
  354.                   IP=IPL-NPTINI+NPTBAS
  355.                ENDIF
  356.                JTBES.NUM(INO,IEL)=IP
  357.  630        CONTINUE
  358.  63      CONTINUE
  359.          IF (IMPR.GE.4) THEN
  360.             write(ioimp,*) 'opto3.eso : topologie amelioree : '
  361.             call ecmai1(jtbes,0)
  362.             segact jtbes*mod
  363.          ENDIF
  364.       ENDIF
  365. *     SEGSUP KELEM
  366.       kelemx.nncou=0
  367. *
  368. * RAZ IDCPX et ICPRX
  369. *      SEGSUP IDCP
  370. *     SEGSUP ICPR
  371.       DO 1500 I=1,IDCPX.JGCOU
  372.          ICPRX.LECTX(IDCPX.LECTX(I))=0
  373.          IDCPX.LECTX(I)=0
  374.  1500 CONTINUE
  375.       CALL mtxadj(ICPRX,0,lchang,'opto3 : ICPRX_sup')
  376.       if (ierr.ne.0) return
  377.       CALL mtxadj(IDCPX,0,lchang,'opto3 : IDCPX_sup')
  378.       if (ierr.ne.0) return
  379. *
  380. * Normal termination
  381. *
  382.       RETURN
  383. *
  384. * Format handling
  385. *
  386.  186  FORMAT ('Segment ',A6,' ',A6,' ajusté de ',I6,' à ',I6)
  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.  
  406.  

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