Télécharger la plaquette Cast3M


fr en

Télécharger topv2.eso

Retour à la liste

Numérotation des lignes : 

topv2

  1. C TOPV2     SOURCE    GOUNAND   25/11/24    21:15:20     12406          
  2.       SUBROUTINE TOPV2(TRAVK,KELEMX,IALGO,IAJNO,XVTOL,QTOL,IMET,IMOMET
  3.      $     ,XDENS,INCMA,ISTMA,KTBES,JCAND,JNASCM,IVERI,impr,TRAVL
  4.      $     ,lchtop)
  5.       IMPLICIT REAL*8 (A-H,O-Z)
  6.       IMPLICIT INTEGER (I-N)
  7. C***********************************************************************
  8. C NOM         : TOPV2
  9. C DESCRIPTION : IJOB=0
  10. C               Minimise le volume d'une topologie de maillage
  11. C               en le maintenant supérieur à 0
  12. C               IJOB=1
  13. C               Minimise le volume, mais on a le droit d'ajouter des
  14. C               noeuds internes
  15. C               IJOB=2
  16. C               La topologie de maillage est supposée être un maillage
  17. C               On essaie de l'améliorer en conservant son volume
  18. C               mais en augmentant sa qualité grace a l'adjonction
  19. C               de noeuds internes
  20. C
  21. * 2017/11/30 : On remplace par ialgo (0 ou 1 : génération ou
  22. *              optimisation de maillage) et iajno (autorise-t-on
  23. *              l'algorithme à ajouter des noeuds.)
  24. C
  25. C
  26. C LANGAGE     : ESOPE
  27. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  28. C               mél : gounand@semt2.smts.cea.fr
  29. C***********************************************************************
  30. C VERSION    : v1, 05/05/2013, version initiale
  31. C HISTORIQUE : v1, 05/05/2013, création
  32. C HISTORIQUE :
  33. C HISTORIQUE :
  34. C***********************************************************************
  35. -INC PPARAM
  36. -INC CCOPTIO
  37. -INC CCREEL
  38. -INC SMELEME
  39. *      POINTEUR KELEM.MELEME
  40.       POINTEUR KEXTO.MELEME
  41.       POINTEUR IBTLOC.MELEME
  42.       POINTEUR IPBTL2.MELEME
  43.       POINTEUR KTBES.MELEME
  44.       POINTEUR KTBES2.MELEME
  45. * anc  POINTEUR IMCAND.MELEME
  46. -INC TMATOP1
  47. *-INC SMELEMX
  48.       POINTEUR LMCANS.MELEMX
  49.       POINTEUR IPBTL.MELEMX
  50.       POINTEUR KELEMX.MELEMX
  51. -INC SMLENTI
  52. *anc      POINTEUR KNNO.MLENTI
  53.       POINTEUR LIDXCA.MLENTI
  54.       POINTEUR LOKVOL.MLENTI
  55.       POINTEUR LNQUAL.MLENTI
  56.       POINTEUR LINDI.MLENTI
  57.       POINTEUR LINDJ.MLENTI
  58. -INC SMLREEL
  59.       POINTEUR IQUAL.MLREEL
  60.       POINTEUR LQUALS.MLREEL
  61. -INC SMCOORD
  62.       POINTEUR KCOORD.MCOORD
  63. *-INC SMETRIQ
  64.       POINTEUR KCMETR.METRIQ
  65. *-INC STRAVJ
  66.       POINTEUR TRAVK.TRAVJ
  67. *-INC STRAVL
  68.  
  69. *
  70.       LOGICAL LOK,LDBG
  71. *anc      LOGICAL LTOIBO
  72. *anc      LOGICAL LTOIBA
  73.       INTEGER JCAND
  74.       LOGICAL LCHANG
  75.       LOGICAL LCHTOP
  76.  
  77. * Liste de topologies de maillages candidates
  78. *      SEGMENT ITCAND(0)
  79. * Liste de topologies de maillages candidats de plus petit volume non nul
  80. *      SEGMENT ITVOL(JG)
  81. * Liste de topologies de maillages candidats de plus petit volume
  82. * et de meilleure qualité
  83. *      SEGMENT ILQUAL(JG)
  84. *      SEGMENT ILIND(JG)
  85. *      SEGMENT JLIND(JG)
  86.       PARAMETER(NGRAV=3)
  87.       DIMENSION XGRAV(NGRAV)
  88.       PARAMETER(NMET=6)
  89.       DIMENSION XMET(NMET)
  90. *
  91. * Executable statements
  92. *
  93.       LDBG=.FALSE.
  94. *      WRITE(IOIMP,*) 'entree topv2'
  95.       IDIMP1=IDIM+1
  96.       KCOORD=TRAVK.COORD
  97.       NKPVIR=TRAVK.PVIRT
  98.       KEXTO=TRAVK.TOPO
  99.       KCMETR=TRAVK.CMETR
  100. *
  101.       LMCANS=TRAVL.MCANS
  102.       LIDXCA=TRAVL.IDXCA
  103.       LOKVOL=TRAVL.OKVOL
  104.       LQUALS=TRAVL.QUALS
  105. *      LNQUAL=TRAVL.NQUAL
  106. *      LINDI=TRAVL.INDI
  107. *      LINDJ=TRAVL.INDJ
  108. *      IPBTL=TRAVL.PBTL
  109. *
  110. * Initialisations SEGMENT TRAVL : nb candidats=0
  111. *                                 nbel maillage candidats=0
  112. *
  113.       CALL TRLADJ(TRAVL,0,0,lchang,'topv2 : TRAVL')
  114.       IF (IERR.NE.0) RETURN
  115. *
  116.       IPOPL2=0
  117.       KTBES=KEXTO
  118.       JCAND=0
  119. * Tests sur le maillage initial : uniquement des éléments à 3 noeuds en 2D
  120. * et des éléments à 4 noeuds en 3D
  121. *      SEGACT KEXTO
  122. *!      NLTLOC=KEXTO.NUM(/2)
  123.       NLTLOC=TRAVK.NVCOU
  124. *      WRITE(IOIMP,*) 'NLTLOC=',NLTLOC
  125.       IF (NLTLOC.EQ.0) GOTO 7
  126.       NBSOUS=KEXTO.LISOUS(/1)
  127.       IF (NBSOUS.NE.0) THEN
  128.          WRITE(IOIMP,*) 'We want only elementary meshes'
  129.          GOTO 9999
  130.       ENDIF
  131.       NBNN=KEXTO.NUM(/1)
  132.       IF (NBNN.NE.IDIM+1) THEN
  133.          WRITE(IOIMP,*) 'Only simplices are allowed'
  134.          GOTO 9999
  135.       ENDIF
  136. *
  137. *      IARET=KELEM.NUM(/1)
  138. *      NBELEM=KELEM.NUM(/2)
  139.       IARET=KELEMX.NNCOU
  140.       NBELEM=KELEMX.NLCOU
  141. *      write(ioimp,*) 'IARET,NBELEM=',IARET,NBELEM
  142. *      segprt,KELEMX
  143. * Une vérif qu'on aurait dû faire depuis longtemps
  144. *      IF (IARET.LT.0.OR.NBELEM.LT.0) THEN
  145. *      IF (IARET.LE.0.OR.NBELEM.LE.0) THEN
  146.       IF (IARET.LE.0.OR.NBELEM.NE.1) THEN
  147.          WRITE(IOIMP,*) 'KELEMX bizarre IARET, NBELEM=',IARET,NBELEM
  148.          goto 9999
  149.       ENDIF
  150.       LOK=.TRUE.
  151.       DO IIAR=1,IARET
  152.          INOD=KELEMX.NUMX(IIAR,1)
  153.          IF (INOD.EQ.0) THEN
  154.             WRITE(IOIMP,*) 'KELEMX noeud nul iiar=',IIAR
  155.             LOK=.FALSE.
  156.          ENDIF
  157.       ENDDO
  158.       IF (.NOT.LOK) THEN
  159.          WRITE(IOIMP,*) 'Ce cas nest pas bon du tout.'
  160.          goto 9999
  161.       ENDIF
  162. * On remet les noeuds KPVIR au barycentre de KELEM (modifie le MCOORD !!!)
  163. * D'après Gruau p.32, il faudrait plutôt le remettre sur S.
  164. *     Mais on pense qu'il est mieux au barycentre de KELEM
  165. *     Aussi, on lui met la métrique
  166. *! 2025119 On ne le fait plus      IF (NKPVIR.NE.0) THEN
  167.       IF (.FALSE.) THEN
  168.          DO L=1,IDIM
  169.             XGRAV(L)=0.D0
  170.          ENDDO
  171.          IF (KCMETR.NE.0) THEN
  172.             DO ININ=1,KCMETR.XIN(/1)
  173.                XMET(ININ)=0.D0
  174.             ENDDO
  175.          ENDIF
  176.          NPOIN=0
  177.          DO 31 IIAR=1,IARET
  178.             INO=KELEMX.NUMX(IIAR,1)
  179.             IF (INO.LE.NKPVIR) GOTO 31
  180.             NPOIN=NPOIN+1
  181.             IREF=IDIMP1*(INO-1)
  182.             DO 5 L=1,IDIM
  183.                XGRAV(L)=XGRAV(L)+KCOORD.XCOOR(IREF+L)
  184.  5          CONTINUE
  185.             IF (KCMETR.NE.0) THEN
  186.                if (ldbg) then
  187.                   write(ioimp,*) 'iiar,ino=',iiar,ino
  188.                   write(ioimp,*) '  coo=',(KCOORD.XCOOR(IREF+L),L=1
  189.      $                 ,IDIM)
  190.                   write(ioimp,*) '  met=',(KCMETR.XIN(L,INO),L=1
  191.      $                 ,KCMETR.XIN(/1))
  192.                endif
  193.                DO 6 ININ=1,KCMETR.XIN(/1)
  194.                   XMET(ININ)=XMET(ININ)+KCMETR.XIN(ININ,INO)
  195.  6             CONTINUE
  196.             ENDIF
  197.  31      CONTINUE
  198.          if (ldbg) write(ioimp,*) 'NPOIN=',NPOIN
  199.          if (ldbg) then
  200.             write(ioimp,*) 'Centre de gravite'
  201.             write(ioimp,*) '  coo=',(xgrav(L)/NPOIN,L=1,IDIM)
  202.             if (kcmetr.ne.0) then
  203.                write(ioimp,*) '  met=',(XMET(L)/NPOIN,L=1,KCMETR.XIN(
  204.      $              /1))
  205.             endif
  206.          endif
  207.          DO IKPVIR=1,NKPVIR
  208.             iref=idimp1*(ikpvir-1)
  209.             do L=1,idim
  210.                kcoord.xcoor(iref+L)=xgrav(L)/NPOIN
  211.             enddo
  212.             IF (KCMETR.NE.0) THEN
  213.                DO 12 ININ=1,KCMETR.XIN(/1)
  214.                   KCMETR.XIN(ININ,IKPVIR)=XMET(ININ)/NPOIN
  215.  12            CONTINUE
  216.             ENDIF
  217.          ENDDO
  218.       ENDIF
  219. * Calcul du volume de la topologie initiale
  220. *  !! pour calculer le volume, il ne faut pas utiliser la métrique !
  221.       IELDEB=1
  222.       IELFIN=TRAVK.NVCOU
  223.       CALL VOMSI2(KEXTO,IELDEB,IELFIN,NKPVIR,XVTLOC,XVTLOS,XVTLOV)
  224.       IF (IERR.NE.0) RETURN
  225. *      WRITE(IOIMP,*) 'Vomsi2 XVTLOC=',XVTLOC,' XVTLOS=',XVTLOS,' XVTLOV
  226. *     $     =',XVTLOV
  227.       IF (NKPVIR.GT.0) THEN
  228.          if (ABS(XVTLOV).GE.XVTOL) THEN
  229.             WRITE(IOIMP,*) '!! Vomsi2 XVTLOV=',XVTLOV,' XVTOL=',XVTOL
  230.             write(ioimp,*) 'NVCOU=',TRAVK.NVCOU,'IELFIN=',IELFIN
  231.             write(ioimp,*) 'NKPVIR=',NKPVIR
  232.             segprt,kexto
  233.             GOTO 7
  234.          endif
  235.       ENDIF
  236.       IF (IALGO.EQ.1) THEN
  237.          if (ABS(XVTLOC-XVTLOS).GE.XVTOL) THEN
  238.             WRITE(IOIMP,*) '!! Vomsi2 XVTLOC=',XVTLOC,' XVTLOS=',XVTLOS
  239.             write(ioimp,*) 'NVCOU=',TRAVK.NVCOU,'IELFIN=',IELFIN
  240.             write(ioimp,*) 'NKPVIR=',NKPVIR
  241.             segprt,kexto
  242.             GOTO 7
  243.          endif
  244.       ENDIF
  245.  
  246. *      JCAND=0
  247. * Note : on saute les topologies locales ayant un volume nul
  248. * car on n'arrivera pas à les améliorer
  249.       IF (XVTLOC.LE.XVTOL) GOTO 7
  250. *     On se demandait si on devait sauter les topologies avec des
  251. *     retournements...      XVTLOC.NE.XVTLOS
  252. *
  253. *     Génération des topologies candidates (stockage dans TRAVL: MCANS
  254. *     indexé par IDXCA)
  255. *     ICBES est le meilleur candidat par défaut (souvent=1, la topologie initiale)
  256. *     IPOPL2 est le candidat où on a créé un nouveau noeud (souvent le
  257. *     dernier, égal à travl.nccou)
  258. *
  259. *      write(ioimp,*) 'Entree topv3'
  260.       CALL TOPV3(TRAVK,KELEMX,IAJNO,TRAVL,INCMA,ISTMA,JNASCM,ICBES
  261.      $     ,IPOPL2,iveri,impr)
  262.       IF (IERR.NE.0) RETURN
  263. *      write(ioimp,*) 'Sortie topv3'
  264. *
  265. * On a les candidats dans ITCAND
  266. * On sélectionne ceux de volume minimum non nul
  267. *
  268.       JCAND=TRAVL.NCCOU
  269.       if (impr.ge.4) then
  270.          write(ioimp,*) 'Apres gen candidat n=',jcand
  271.       endif
  272. * Sortie anticipée s'il n'y a qu'un candidat
  273. *      WRITE(IOIMP,*) 'JCAND=',JCAND
  274.       IF (JCAND.EQ.1) GOTO 8
  275. *
  276.       XVMIN=XVTLOC
  277.       if (impr.ge.4) then
  278.          write(ioimp,*) 'volume vise XVMIN=',XVMIN
  279.       endif
  280. *
  281. *      write(ioimp,*) 'Entree topv4'
  282.       CALL TOPV4(TRAVL,IALGO,XVMIN,XVTOL,NKPVIR,impr)
  283.       IF (IERR.NE.0) RETURN
  284. *      write(ioimp,*) 'Sortie topv4'
  285. * S'il n'en reste pas (ne doit pas se produire)
  286.       IF (TRAVL.NVOCOU.LE.0) THEN
  287.          write(ioimp,*) 'Plus de candidats apres topv4 ???'
  288.          goto 9999
  289. *         ICBES=0
  290. *         GOTO 666
  291.       ENDIF
  292. * S'il n'en reste qu'un, on peut sauter le calcul des qualités ??????
  293.       IF (TRAVL.NVOCOU.EQ.1) THEN
  294.          ICBES=LOKVOL.LECT(NVOCOU)
  295.          GOTO 8
  296.       ENDIF
  297. *
  298. *     Les candidats ayant le bon volume sont dans ITVOL
  299. *     On calcule les qualités de chaque élément des candidats et on ordonne.
  300. *
  301. *      write(ioimp,*) 'Entree topv5'
  302.       CALL TOPV5(TRAVL,XVTOL,IMET,IMOMET,XDENS,KCMETR,NKPVIR)
  303.       IF (IERR.NE.0) RETURN
  304. *      write(ioimp,*) 'Sortie topv5'
  305.  
  306.       if (impr.ge.4) then
  307.          write(ioimp,*) 'Apres calcul qualité candidats'
  308.          DO i=1,NVOCOU
  309.             Write(ioimp,*) 'Candidat i=',i
  310.             ICAND=LOKVOL.LECT(i)
  311.             IELDEB=LIDXCA.LECT(ICAND)
  312.             IELFIN=LIDXCA.LECT(ICAND+1)-1
  313.             Write(ioimp,*) (lquals.prog(iii),iii=ieldeb,ielfin)
  314. *     $           ,lnqual.lect(icand))
  315.          enddo
  316.       endif
  317. *
  318. * Calcul des meilleurs candidats par maximum lexical
  319. *
  320. *      write(ioimp,*) 'Entree topv6'
  321.       CALL TOPV6(TRAVL,QTOL,ICBES)
  322. *      write(ioimp,*) 'Sortie topv6'
  323. *
  324.       if (impr.ge.4) then
  325.          write(ioimp,*) 'Apres tri qualité candidats, le meilleur est :'
  326.      $        ,icbes
  327.       endif
  328.  
  329.  8    CONTINUE
  330. *
  331. *anc      KTBES2=KTBES
  332.       ICAND=ICBES
  333. *dbg      write(ioimp,*) 'icand,nccou=',icand,travl.nccou
  334.       IF (ICAND.NE.1) THEN
  335.          IELDEB=LIDXCA.LECT(ICAND)
  336.          IELFIN=LIDXCA.LECT(ICAND+1)-1
  337.          NBNN=LMCANS.NNCOU
  338.          NBELEM=IELFIN-IELDEB+1
  339.          NBSOUS=0
  340.          NBREF=0
  341. *            write(ioimp,*) 'icand,ieldeb,ielfin,nbnn,nbelem=',icand
  342. *     $           ,ieldeb,ielfin,nbnn,nbelem
  343.          SEGINI KTBES
  344.          KTBES.ITYPEL=LMCANS.ITYPEX
  345.          IDX=1
  346.          DO IEL=IELDEB,IELFIN
  347.             DO INO=1,NBNN
  348.                KTBES.NUM(INO,IDX)=LMCANS.NUMX(INO,IEL)
  349.             ENDDO
  350.             IDX=IDX+1
  351.          ENDDO
  352.       ENDIF
  353. *
  354.  666  CONTINUE
  355. * Si on n'a pas sélectionné le candidat avec le noeud supplémentaire
  356. * on peut enlever le noeud
  357. *ijob      IF (IJOB.EQ.1.OR.IJOB.EQ.2) THEN
  358.       IF (IAJNO.NE.0) THEN
  359.          IF (IPOPL2.NE.0.AND.ICBES.NE.IPOPL2) THEN
  360.             NPCOUN=TRAVK.NPCOU-1
  361. *     Remise à zéro : nécessaire uniquement pour la vérification ??
  362.             if (iveri.ge.1) then
  363.                IREF=(TRAVK.NPCOU-1)*IDIMP1
  364.                DO 111 I=1,IDIMP1
  365.                   KCOORD.XCOOR(IREF+I)=0.D0
  366.  111           CONTINUE
  367.                IF (KCMETR.NE.0) THEN
  368.                   DO 112 ININ=1,KCMETR.XIN(/1)
  369.                      KCMETR.XIN(ININ,TRAVK.NPCOU)=0.D0
  370.  112              CONTINUE
  371.                ENDIF
  372.             endif
  373.             TRAVK.NPCOU=NPCOUN
  374. * Désactivation temporaire car pas de ménage...
  375. *tmp            if (iveri.ge.2) then
  376. *tmp            call vetopi(travk,'topv2 : Apres reduction travk')
  377. *tmp               if (ierr.ne.0) return
  378. *tmp            endif
  379.          ELSE
  380. *            write(ioimp,*) 'Nouveau noeud cree ','IPOPL2,ICBES,NPCOU='
  381. *     $           ,IPOPL2,ICBES,TRAVK.NPCOU
  382. *            IREF=(TRAVK.NPCOU-1)*IDIMP1
  383. *            write(ioimp,*)  (kcoord.xcoor(iref+iii),iii=1,idimp1)
  384.          ENDIF
  385.       ENDIF
  386.  7    CONTINUE
  387.       IF (KEXTO.EQ.KTBES) JCAND=-JCAND
  388.       IF (KEXTO.EQ.KTBES) THEN
  389.          LCHTOP=.FALSE.
  390.       ELSE
  391.          LCHTOP=.TRUE.
  392.       ENDIF
  393. *      WRITE(IOIMP,*) 'sortie topv2'
  394.       RETURN
  395. *
  396. *
  397. *
  398.  9999 CONTINUE
  399.       MOTERR(1:8)='TOPV2   '
  400. * 349 2
  401. *Problème non prévu dans le s.p. %m1:8 contactez votre assistance
  402.       CALL ERREUR(349)
  403.       RETURN
  404. *
  405. * End of subroutine TOPV2
  406. *
  407.       END
  408.  
  409.  

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