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

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