Télécharger la plaquette Cast3M


fr en

Télécharger topv3.eso

Retour à la liste

Numérotation des lignes : 

topv3

  1. C TOPV3     SOURCE    GOUNAND   21/04/06    21:15:34     10940          
  2. * On préférerait KEXTO à la place de TRAVK mais KEXTO n'est pas autoporteur.
  3. *ijob      SUBROUTINE TOPV3(TRAVK,KELEM,IJOB,TRAVL,ICBES
  4. *kelemx      SUBROUTINE TOPV3(TRAVK,KELEM,IAJNO,TRAVL,INCMA,ISTMA,
  5.       SUBROUTINE TOPV3(TRAVK,KELEMX,IAJNO,TRAVL,INCMA,ISTMA,
  6.      $     JNASCM,ICBES,IPOPL2,iveri,impr) 
  7.       IMPLICIT REAL*8 (A-H,O-Z)
  8.       IMPLICIT INTEGER (I-N)
  9. C***********************************************************************
  10. C NOM         : TOPV3
  11. C DESCRIPTION :
  12. *
  13. *     Génération des topologies candidates (stockage dans LMCANS indexé
  14. *     par LIDXCA) Issue de topv2_nettoie_final.eso 
  15. *      
  16. C     
  17. C     
  18. C LANGAGE     : ESOPE
  19. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  20. C               mél : gounand@semt2.smts.cea.fr
  21. C***********************************************************************
  22. C VERSION    : v1, 09/11/2017, version initiale
  23. C HISTORIQUE : v1, 09/11/2017, création
  24. C HISTORIQUE :
  25. C HISTORIQUE :
  26. C***********************************************************************
  27. -INC PPARAM
  28. -INC CCOPTIO
  29. *-INC TMATOP2
  30. -INC CCREEL
  31. -INC SMELEME
  32. *      POINTEUR KELEM.MELEME
  33.       POINTEUR KEXTO.MELEME
  34.       POINTEUR IBTLOC.MELEME
  35.       POINTEUR IPBTL2.MELEME
  36.       POINTEUR KTBES.MELEME
  37.       POINTEUR KTBES2.MELEME
  38.       POINTEUR IMCAND.MELEME
  39. -INC TMATOP1
  40. *-INC SMELEMX
  41.       POINTEUR LMCANS.MELEMX
  42.       POINTEUR IPBTL.MELEMX
  43.       POINTEUR KELEMX.MELEMX
  44. -INC SMLENTI
  45.       POINTEUR KNNO.MLENTI
  46.       POINTEUR LIDXCA.MLENTI
  47. *      POINTEUR LOKVOL.MLENTI
  48. *      POINTEUR LNQUAL.MLENTI
  49. *      POINTEUR LINDI.MLENTI
  50. *      POINTEUR LINDJ.MLENTI
  51. -INC SMLREEL
  52. *      POINTEUR IQUAL.MLREEL
  53. *      POINTEUR LQUALS.MLREEL
  54. -INC SMCOORD
  55. *-INC STRAVJ
  56.       POINTEUR TRAVK.TRAVJ
  57. *-INC STRAVL
  58. *     
  59.       LOGICAL LOK
  60.       LOGICAL LTOIBO
  61.       LOGICAL LTOIBA
  62.       LOGICAL LLIMCA
  63.       INTEGER JCAND
  64.       LOGICAL LCHANG
  65.       LOGICAL LCHTOP
  66. * Liste de topologies de maillages candidates
  67. *      SEGMENT ITCAND(0)
  68. * Liste de topologies de maillages candidats de plus petit volume non nul
  69. *      SEGMENT ITVOL(JG)
  70. * Liste de topologies de maillages candidats de plus petit volume
  71. * et de meilleure qualité
  72. *      SEGMENT ILQUAL(JG)
  73. *      SEGMENT ILIND(JG)
  74. *      SEGMENT JLIND(JG)
  75. *     
  76. * Executable statements
  77. *     
  78. *      WRITE(IOIMP,*) 'coucou topv3'
  79. *      IDIMP1=IDIM+1
  80.       KEXTO=TRAVK.TOPO
  81.       KPVIRT=TRAVK.PVIRT
  82. *
  83.       LMCANS=TRAVL.MCANS
  84.       LIDXCA=TRAVL.IDXCA
  85. *      LOKVOL=TRAVL.OKVOL
  86. *      LQUALS=TRAVL.QUALS
  87. *      LNQUAL=TRAVL.NQUAL
  88. *      LINDI=TRAVL.INDI
  89. *      LINDJ=TRAVL.INDJ
  90.       IPBTL=TRAVL.PBTL
  91. *     Les noeud S et S' de Gruau p.42
  92. *kelemx      IARET=KELEM.NUM(/1)
  93.       IARET=KELEMX.NNCOU
  94. *      
  95. *      IS=KELEM.NUM(1,1)
  96.       IS=KELEMX.NUMX(1,1)
  97.       ISP=0
  98.       IS3=0
  99.       IS4=0
  100. *kelemx      IF (IARET.EQ.2) ISP=KELEM.NUM(2,1)
  101. *kelemx      IF (IARET.EQ.3) IS3=KELEM.NUM(3,1)
  102. *kelemx      IF (IARET.EQ.4) IS4=KELEM.NUM(4,1)
  103.       IF (IARET.EQ.2) ISP=KELEMX.NUMX(2,1)
  104.       IF (IARET.EQ.3) IS3=KELEMX.NUMX(3,1)
  105.       IF (IARET.EQ.4) IS4=KELEMX.NUMX(4,1)
  106. *      
  107. * Le premier candidat est toujours l'original qui n'est pas forcément un étoilement
  108. *anc      SEGINI ITCAND
  109. *      ITCAND(**)=KEXTO
  110. *      
  111.       NCCOUO=TRAVL.NCCOU
  112.       NLCOUO=LMCANS.NLCOU
  113.       NNC=NCCOUO+1
  114. *      NNL=NLCOUO+KEXTO.NUM(/2)
  115.       NNL=NLCOUO+TRAVK.NVCOU
  116.       CALL TRLADJ(TRAVL,NNC,NNL,lchang,'topv3 : TRAVL')
  117.       if (ierr.ne.0) return
  118.       IDX=LIDXCA.LECT(NNC)
  119. *      DO IEL=1,KEXTO.NUM(/2)
  120.       DO IEL=1,TRAVK.NVCOU
  121.          DO INO=1,KEXTO.NUM(/1)
  122.             LMCANS.NUMX(INO,IDX)=KEXTO.NUM(INO,IEL)
  123.          ENDDO
  124.          IDX=IDX+1
  125.       ENDDO
  126.       LIDXCA.LECT(NNC+1)=IDX
  127.       ICBES=1
  128.       if (iveri.ge.2) then
  129.          call trlver(travl,'topv3 : Apres initialisation KEXTO')
  130.          if (ierr.ne.0) return
  131.       endif
  132. * Extraction du bord (contour ou enveloppe)
  133. *      write(ioimp,*) 'Avant extraction bord'
  134.       IF (IDIM.EQ.2) THEN
  135. c$$$*         call ecmai1(kexto,0)
  136. c$$$         CALL ECRCHA('NOID')
  137. c$$$         CALL ECRCHA('TOUT')
  138. c$$$         CALL ECROBJ('MAILLAGE',KEXTO)
  139. c$$$         CALL PRCONT
  140. c$$$         CALL LIROBJ('MAILLAGE',IBTLOC,1,IRET)
  141. c$$$         IF (IERR.NE.0) RETURN
  142. c$$$         SEGACT KEXTO*MOD
  143. c$$$         SEGACT IBTLOC
  144. c$$$*     prcont le passe en SEGACT non *MOD ce qui pose ensuite problème à
  145. c$$$*     baryc4.eso
  146. c$$$         SEGACT MCOORD*MOD
  147. *
  148.          IELDEB=1
  149.          IELFIN=TRAVK.NVCOU
  150.          ICPR=0
  151.          IDCP=0
  152.          NPLOC=TRAVK.NPCOU
  153. *         ITYCON=1
  154.          ITYCON=3
  155.          INOID=1
  156. *         CALL CONTOO(KEXTO,IELDEB,IELFIN,ICPR,IDCP,NPLOC,ITYCON,INOID
  157. *     $        ,IBTLOC) 
  158.          CALL CONTOU(KEXTO,IELDEB,IELFIN,ICPR,IDCP,NPLOC,ITYCON,INOID
  159.      $        ,IBTLOC) 
  160.          IF (IERR.NE.0) RETURN
  161.          SEGACT IBTLOC
  162.       ELSEIF (IDIM.EQ.3) THEN
  163. *
  164.          IELDEB=1
  165.          IELFIN=TRAVK.NVCOU
  166.          ICLE=0
  167.          INOID=1
  168.          CALL ENVVO3(KEXTO,IELDEB,IELFIN,ICLE,INOID,IBTLOC)
  169.          IF (IERR.NE.0) RETURN
  170.       ELSE
  171. * 709 2
  172. *Fonction indisponible en dimension %i1.
  173.          INTERR(1)=IDIM
  174.          CALL ERREUR(709)
  175.       ENDIF
  176.       IF (IERR.NE.0) RETURN
  177.  
  178.  
  179. *      WRITE(IOIMP,*) 'KEXTO'
  180. *      CALL ECMAI1(kexto,0)
  181.       if (impr.ge.4) then
  182.          write(ioimp,*) 'KPVIRT=',KPVIRT
  183.          write(ioimp,*) 'Apres extraction bord IBTLOC=',IBTLOC
  184.          WRITE(IOIMP,*) 'IBTLOC'
  185.          CALL ECMAI1(ibtloc,0)
  186.          SEGACT IBTLOC
  187.       endif
  188. *
  189. *     On supprimme la référence au contour car on supprimme le contour
  190. *     après...
  191. *     Pas utile ici car on est en numérotation locale pour KEXTO...
  192. * 
  193.       NLBTL=IBTLOC.NUM(/2)
  194. * Il arrive quelquefois que la topologie locale n'ait pas de bord
  195.       IF (NLBTL.GT.0) THEN
  196. * Si la topologie locale n'a qu'un seul élément, il n'est pas nécessaire
  197. * de l'étoiler
  198.          NLTLOC=TRAVK.NVCOU
  199. *         
  200.          LTOIBO=(NLTLOC.GT.1)
  201. *ijob         LTOIBA=(IJOB.EQ.1.OR.IJOB.EQ.2)
  202.          LTOIBA=(IAJNO.EQ.1)
  203. * Si on doit etoiler, on contruit le maillage des points du bord
  204. * = chan IBTLOC 'POI1'
  205. * on applique ici une méthode locale en O(n^2) ce qui suppose que IBTLOC
  206. * n'a pas trop de points...
  207.          IF (LTOIBO.OR.LTOIBA) THEN
  208.             KNNO=TRAVK.NNO
  209.             NBELEM=IBTLOC.NUM(/2)
  210. *         WRITE(IOIMP,*) 'NBELEM=',NBELEM
  211. *         IF (NBELEM.GT.100) THEN
  212. *            WRITE(IOIMP,*) 'KEXTO'
  213. *            CALL ECMAI1(KEXTO,0)
  214. *            SEGACT KEXTO
  215. *            WRITE(IOIMP,*) 'IBTLOC'
  216. *            CALL ECMAI1(IBTLOC,0)
  217. *            SEGACT IBTLOC
  218. *            STOP 16
  219. *         ENDIF
  220.             NBNN=IBTLOC.NUM(/1)
  221.             IK=0
  222.             DO IBELEM=1,NBELEM
  223.                DO IBNN=1,NBNN
  224.                   INO=IBTLOC.NUM(IBNN,IBELEM)
  225.                   if (ino.eq.0) then
  226.                      write(ioimp,*) 'Noeud nul détecté !!!!'
  227.                      WRITE(IOIMP,*) 'KEXTO'
  228.                      call ecmai1(kexto,0)
  229.                      WRITE(IOIMP,*) 'IBTLOC'
  230.                      CALL ECMAI1(ibtloc,0)
  231.                      goto 9999
  232.                   endif
  233.                   IF (KNNO.LECT(INO).EQ.0) THEN
  234.                      IK=IK+1
  235.                      KNNO.LECT(INO)=IK
  236.                   ENDIF
  237.                ENDDO
  238.             ENDDO
  239.             CALL mlxadl(IPBTL,IK,lchang,'topv3 : IPBTL_IK')
  240.             if (ierr.ne.0) return
  241.             if (iveri.ge.2) then
  242.                call vemelx(ipbtl,'topv3 : Apres requisition ipbtl')
  243.                if (ierr.ne.0) return
  244.             endif           
  245. * On regarde également si IS ou ISP font partie du bord
  246.             IS2=IS
  247.             ISP2=ISP
  248.             IS32=IS3
  249.             IS42=IS4
  250.             DO IIPO=1,TRAVK.NPCOU
  251. *     DO IIPO=1,IPO(/1)
  252.                INLOC=KNNO.LECT(IIPO)
  253.                IF (INLOC.NE.0) THEN
  254.                   IPBTL.NUMX(1,INLOC)=IIPO
  255.                   IF (IS2.EQ.IIPO) IS2=0
  256.                   IF (ISP2.EQ.IIPO) ISP2=0
  257.                   IF (IS32.EQ.IIPO) IS32=0
  258.                   IF (IS42.EQ.IIPO) IS42=0
  259. *     Nettoyage de KNNO
  260.                   KNNO.LECT(IIPO)=0
  261.                ENDIF
  262.             ENDDO
  263. * Vérification du nettoyage de KNNO
  264.             if (iveri.ge.2) then
  265.                call vetopi(travk,
  266.      $              'topv3 : Apres creation points du bord')
  267.                if (ierr.ne.0) return
  268.             endif            
  269.             IF (IVERI.GE.2.and..false.) THEN
  270. * à corriger pour le nouveau ipbtl en melemx
  271.                IPBTL2=IBTLOC
  272.                CALL CHANGE(IPBTL2,1)
  273.                SEGACT IBTLOC
  274.                CALL OUEXCL(IPBTL,IPBTL2,IPT3)
  275.                IF (IERR.NE.0) RETURN
  276.                SEGACT IPBTL
  277.                SEGACT MCOORD*MOD
  278.                IF (IPT3.NE.0) THEN
  279.                   WRITE(IOIMP,*) 'IPT3 pour IPBTL'
  280.                   CALL ECMAI1(IPT3,0)
  281.                   IF (IERR.NE.0) RETURN
  282.                   WRITE(IOIMP,*) 'NEL1=',IPBTL.NLCOU
  283.                   CALL ECMELX(IPBTL,0)
  284.                   SEGACT IPBTL2
  285.                   WRITE(IOIMP,*) 'NEL2=',IPBTL2.NUM(/2)
  286.                   CALL ECMAI1(IPBTL2,0)
  287.                   CALL ERREUR(5)
  288.                   RETURN
  289.                ENDIF
  290.             ENDIF
  291. * On étoile à partir des éléments du bord
  292. * SEGACT en trop ? On est méfiant...
  293. *            SEGACT IBTLOC
  294.             IF (LTOIBO) THEN
  295. * On étoile à partir de S ou S' s'ils ne font pas partie du bord
  296.                DO IBIS=1,4
  297.                   IF (IBIS.EQ.1) THEN
  298.                      NODE=IS2
  299.                      MOTERR(1:4)='IS2 '
  300.                   ELSEIF (IBIS.EQ.2) THEN
  301.                      NODE=ISP2
  302.                      MOTERR(1:4)='ISP2'
  303.                   ELSEIF (IBIS.EQ.3) THEN
  304.                      NODE=IS32
  305.                      MOTERR(1:4)='IS32'
  306.                   ELSEIF (IBIS.EQ.4) THEN
  307.                      NODE=IS42
  308.                      MOTERR(1:4)='IS42'
  309.                   ELSE
  310.                      write(ioimp,*) 'pb boucle ibis'
  311.                      goto 9999
  312.                   ENDIF
  313.                   IF (NODE.NE.0) THEN
  314. *
  315.                      CALL ETOIL2(NODE,IBTLOC,TRAVL)
  316.                      IF (IERR.NE.0) RETURN
  317.                      if (iveri.ge.2) then
  318.                         call trlver(travl
  319.      $                       ,'topv3 : Apres etoil2, IBIS')
  320.                         if (ierr.ne.0) return
  321.                      endif
  322.                      ncc=travl.nccou
  323.                      if (lidxca.lect(ncc+1).eq.lidxca.lect(ncc)) goto
  324.      $                    666
  325.                   ENDIF
  326.                ENDDO
  327. *anc!!!  NPBTL=IPBTL.NUM(/2)
  328.                NPBTL=IPBTL.NLCOU
  329. *     WRITE(IOIMP,*) 'NPBTL=',NPBTL
  330. *               NPBTLM=10000000
  331.                IF (NPBTL.GT.INCMA) THEN
  332.                   LLIMCA=.TRUE.
  333.                   JNASCM=JNASCM+1
  334.                   IF (ISTMA.EQ.0) THEN
  335.                      NPBTLR=0
  336.                      LTOIBA=.FALSE.
  337.                   ELSEIF (ISTMA.EQ.1) THEN
  338.                      NPBTLR=1
  339.                      JNPBTL=(NPBTL+1)/2
  340.                   ELSEIF (ISTMA.EQ.2) THEN
  341. * Attention overflow potentiel...                     
  342.                      NPBTLR=MAX(1,NINT(INCMA*(DBLE(INCMA)/DBLE(NPBTL))))
  343.                      JNPBTL=(NPBTL+1)/2
  344.                   ELSE
  345.                      WRITE(IOIMP,*) 'ISTMA=',ISTMA,' non prevu'
  346.                      GOTO 9999
  347.                   ENDIF
  348.                   if (impr.ge.2) then 
  349.                      write(ioimp,*) 'topv3 : reduction nb cand de '
  350.      $                    ,NPBTL,' a ',NPBTLR
  351.                   endif
  352.                ELSE
  353.                   LLIMCA=.FALSE.
  354.                   NPBTLR=NPBTL
  355.                ENDIF
  356.  
  357.                DO INPBTL=1,NPBTLR
  358.                   IF (.NOT.LLIMCA) THEN
  359.                      NODE=IPBTL.NUMX(1,INPBTL)
  360.                   ELSE
  361.                      IF (ISTMA.EQ.1) THEN
  362.                         NODE=IPBTL.NUMX(1,JNPBTL)
  363.                      ELSEIF (ISTMA.EQ.2) THEN
  364.                         IF (NPBTLR.NE.1) JNPBTL=1+NINT((NPBTLR-1)
  365.      $                       *(DBLE(INPBTL-1)/DBLE(NPBTLR-1))) 
  366.                         NODE=IPBTL.NUMX(1,JNPBTL)
  367.                      ELSE
  368.                         WRITE(IOIMP,*) 'ISTMA=',ISTMA,' non prevu 2'
  369.                         GOTO 9999
  370.                      ENDIF
  371.                   ENDIF
  372.  
  373. *            WRITE(IOIMP,*) 'INPBTL=',INPBTL,' NODE=',NODE
  374.                   MOTERR(1:4)='NBOR'
  375. *
  376. *                  write(ioimp,*) 'lmcans avant'
  377. *                  call ecmelx(lmcans,0)
  378.                   CALL ETOIL2(NODE,IBTLOC,TRAVL)
  379.                   IF (IERR.NE.0) RETURN
  380. *                  write(ioimp,*) 'lmcans apres'
  381. *                  call ecmelx(lmcans,0)
  382.                   if (iveri.ge.2) then
  383.                      call trlver(travl
  384.      $                    ,'topv3 : Apres etoil2, INPBTL')
  385.                      if (ierr.ne.0) return
  386.                   endif
  387.                   ncc=travl.nccou
  388.                   if (lidxca.lect(ncc+1).eq.lidxca.lect(ncc)) goto
  389.      $                 666
  390.                ENDDO
  391.             ENDIF
  392.             IF (LTOIBA) THEN 
  393. * Cas 1 : on étoile avec l'isobarycentre du contour
  394.                IF (IARET.EQ.1) THEN
  395.                   CALL BARYC5(IPBTL,KPVIRT,TRAVK,NODE)
  396.                   MOTERR(1:4)='BARC'
  397. * Cas 2 : on étoile avec l'isobarycentre de S et S'
  398.                ELSEIF (IARET.EQ.2) THEN
  399. *kelemx                  CALL BARYC4(KELEM,0,TRAVK,NODE)
  400.                   CALL BARYC5(KELEMX,0,TRAVK,NODE)
  401.                   MOTERR(1:4)='BARS'
  402. * Cas 3 ajout 2017/08/22                  
  403.                ELSEIF (IARET.EQ.3) THEN
  404. *kelemx                  CALL BARYC4(KELEM,0,TRAVK,NODE)
  405.                   CALL BARYC5(KELEMX,0,TRAVK,NODE)
  406.                   MOTERR(1:4)='BAR3'
  407. * Cas 4 ajout 2017/08/22                  
  408.                ELSEIF (IARET.EQ.4) THEN
  409. *kelemx                  CALL BARYC4(KELEM,0,TRAVK,NODE)
  410.                   CALL BARYC5(KELEMX,0,TRAVK,NODE)
  411.                   MOTERR(1:4)='BAR4'
  412.                ELSE
  413.                   Write(ioimp,*) 'iaret=',iaret
  414.                   call erreur(5)
  415.                   return
  416.                ENDIF
  417.                IF (IERR.NE.0) RETURN
  418. *
  419.                CALL ETOIL2(NODE,IBTLOC,TRAVL)
  420.                IF (IERR.NE.0) RETURN
  421.                if (iveri.ge.2) then
  422.                   call trlver(travl
  423.      $                 ,'topv3 : Apres etoil2, BARYC')
  424.                   if (ierr.ne.0) return
  425.                endif
  426.                ipopl2=travl.nccou
  427.                ncc=travl.nccou
  428.                if (lidxca.lect(ncc+1).eq.lidxca.lect(ncc)) goto
  429.      $              666
  430.             ENDIF
  431. *     SEGSUP IPBTL
  432. * Ne le faire que si iveri=1 ?
  433.             if (iveri.ge.1) then
  434.                DO IZER=1,IPBTL.NLCOU
  435.                   IPBTL.NUMX(1,IZER)=0
  436.                ENDDO
  437.             endif
  438.             CALL mlxadl(IPBTL,0,lchang,'topv3 : IPBTL_0')
  439.             if (ierr.ne.0) return
  440.             if (iveri.ge.2) then
  441.                call vemelx(ipbtl,'topv3 : Apres nettoyage ipbtl')
  442.                if (ierr.ne.0) return
  443.             endif           
  444.          ENDIF
  445.       ENDIF
  446.       SEGSUP IBTLOC
  447.       RETURN
  448. *     
  449. *     
  450. *     
  451.  9999 CONTINUE
  452.       MOTERR(1:8)='TOPV3   '
  453. * 349 2
  454. *Problème non prévu dans le s.p. %m1:8 contactez votre assistance
  455.       CALL ERREUR(349)
  456.       RETURN
  457.  666  CONTINUE
  458.       WRITE(IOIMP,*) 'topv3 : Pb candidat ',MOTERR(1:4)
  459. *a upgrader      CALL ECMAI1(IMCAND,0)
  460.       WRITE(IOIMP,*) 'KEXTO'
  461.       CALL ECMAI1(KEXTO,0)
  462.       WRITE(IOIMP,*) 'IBTLOC'
  463.       CALL ECMAI1(IBTLOC,0)
  464.       WRITE(IOIMP,*) 'IPBTL'
  465.       CALL ECMELX(IPBTL,0)
  466.       WRITE(IOIMP,*) 'NODE=',NODE
  467.       CALL ERREUR(5)
  468.       RETURN
  469. *     
  470. * End of subroutine TOPV3
  471. *     
  472.       END
  473.  
  474.  
  475.  

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