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   26/06/09    21:15:19     12566          
  2.       SUBROUTINE TOPV2(TRAVK,KELEMX,IALGO,IAJNO,XVTOL,QTOL,IMET,IMOMET
  3.      $     ,XDENS,INCMA,ISTMA,KTBES,JCAND,JNASCM,IVERI,impr,TRAVL
  4.      $     ,lchtop,jcritq,pcritq,qcritq)
  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.       LOGICAL LCQUAL
  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.       CHARACTER*40 CFORM,CFORM2
  91. *
  92. * Executable statements
  93. *
  94.       LDBG=.FALSE.
  95. *      WRITE(IOIMP,*) 'entree topv2'
  96.       IDIMP1=IDIM+1
  97.       KCOORD=TRAVK.COORD
  98.       NKPVIR=TRAVK.PVIRT
  99.       KEXTO=TRAVK.TOPO
  100.       KCMETR=TRAVK.CMETR
  101. *
  102.       LMCANS=TRAVL.MCANS
  103.       LIDXCA=TRAVL.IDXCA
  104.       LOKVOL=TRAVL.OKVOL
  105.       LQUALS=TRAVL.QUALS
  106.       LNQUAL=TRAVL.NQUAL
  107. *      LINDI=TRAVL.INDI
  108. *      LINDJ=TRAVL.INDJ
  109. *      IPBTL=TRAVL.PBTL
  110.       LCQUAL=.FALSE.
  111. *
  112. * Initialisations SEGMENT TRAVL : nb candidats=0
  113. *                                 nbel maillage candidats=0
  114. *
  115.       CALL TRLADJ(TRAVL,0,0,lchang,'topv2 : TRAVL')
  116.       IF (IERR.NE.0) RETURN
  117. *
  118.       IPOPL2=0
  119.       KTBES=KEXTO
  120.       JCAND=0
  121. * Tests sur le maillage initial : uniquement des éléments à 3 noeuds en 2D
  122. * et des éléments à 4 noeuds en 3D
  123. *      SEGACT KEXTO
  124. *!      NLTLOC=KEXTO.NUM(/2)
  125.       NLTLOC=TRAVK.NVCOU
  126. *      WRITE(IOIMP,*) 'NLTLOC=',NLTLOC
  127.       IF (NLTLOC.EQ.0) GOTO 7
  128.       NBSOUS=KEXTO.LISOUS(/1)
  129.       IF (NBSOUS.NE.0) THEN
  130.          WRITE(IOIMP,*) 'We want only elementary meshes'
  131.          GOTO 9999
  132.       ENDIF
  133.       NBNN=KEXTO.NUM(/1)
  134.       IF (NBNN.NE.IDIM+1) THEN
  135.          WRITE(IOIMP,*) 'Only simplices are allowed'
  136.          GOTO 9999
  137.       ENDIF
  138. *
  139. *      IARET=KELEM.NUM(/1)
  140. *      NBELEM=KELEM.NUM(/2)
  141.       IARET=KELEMX.NNCOU
  142.       NBELEM=KELEMX.NLCOU
  143. *      write(ioimp,*) 'IARET,NBELEM=',IARET,NBELEM
  144. *      segprt,KELEMX
  145. * Une vérif qu'on aurait dû faire depuis longtemps
  146. *      IF (IARET.LT.0.OR.NBELEM.LT.0) THEN
  147. *      IF (IARET.LE.0.OR.NBELEM.LE.0) THEN
  148.       IF (IARET.LE.0.OR.NBELEM.NE.1) THEN
  149.          WRITE(IOIMP,*) 'KELEMX bizarre IARET, NBELEM=',IARET,NBELEM
  150.          goto 9999
  151.       ENDIF
  152.       LOK=.TRUE.
  153.       DO IIAR=1,IARET
  154.          INOD=KELEMX.NUMX(IIAR,1)
  155.          IF (INOD.EQ.0) THEN
  156.             WRITE(IOIMP,*) 'KELEMX noeud nul iiar=',IIAR
  157.             LOK=.FALSE.
  158.          ENDIF
  159.       ENDDO
  160.       IF (.NOT.LOK) THEN
  161.          WRITE(IOIMP,*) 'Ce cas nest pas bon du tout.'
  162.          goto 9999
  163.       ENDIF
  164. * On remet les noeuds KPVIR au barycentre de KELEM (modifie le MCOORD !!!)
  165. * D'après Gruau p.32, il faudrait plutôt le remettre sur S.
  166. *     Mais on pense qu'il est mieux au barycentre de KELEM
  167. *     Aussi, on lui met la métrique
  168. *! 2025119 On ne le fait plus      IF (NKPVIR.NE.0) THEN
  169.       IF (.FALSE.) THEN
  170.          DO L=1,IDIM
  171.             XGRAV(L)=0.D0
  172.          ENDDO
  173.          IF (KCMETR.NE.0) THEN
  174.             DO ININ=1,KCMETR.XIN(/1)
  175.                XMET(ININ)=0.D0
  176.             ENDDO
  177.          ENDIF
  178.          NPOIN=0
  179.          DO 31 IIAR=1,IARET
  180.             INO=KELEMX.NUMX(IIAR,1)
  181.             IF (INO.LE.NKPVIR) GOTO 31
  182.             NPOIN=NPOIN+1
  183.             IREF=IDIMP1*(INO-1)
  184.             DO 5 L=1,IDIM
  185.                XGRAV(L)=XGRAV(L)+KCOORD.XCOOR(IREF+L)
  186.  5          CONTINUE
  187.             IF (KCMETR.NE.0) THEN
  188.                if (ldbg) then
  189.                   write(ioimp,*) 'iiar,ino=',iiar,ino
  190.                   write(ioimp,*) '  coo=',(KCOORD.XCOOR(IREF+L),L=1
  191.      $                 ,IDIM)
  192.                   write(ioimp,*) '  met=',(KCMETR.XIN(L,INO),L=1
  193.      $                 ,KCMETR.XIN(/1))
  194.                endif
  195.                DO 6 ININ=1,KCMETR.XIN(/1)
  196.                   XMET(ININ)=XMET(ININ)+KCMETR.XIN(ININ,INO)
  197.  6             CONTINUE
  198.             ENDIF
  199.  31      CONTINUE
  200.          if (ldbg) write(ioimp,*) 'NPOIN=',NPOIN
  201.          if (ldbg) then
  202.             write(ioimp,*) 'Centre de gravite'
  203.             write(ioimp,*) '  coo=',(xgrav(L)/NPOIN,L=1,IDIM)
  204.             if (kcmetr.ne.0) then
  205.                write(ioimp,*) '  met=',(XMET(L)/NPOIN,L=1,KCMETR.XIN(
  206.      $              /1))
  207.             endif
  208.          endif
  209.          DO IKPVIR=1,NKPVIR
  210.             iref=idimp1*(ikpvir-1)
  211.             do L=1,idim
  212.                kcoord.xcoor(iref+L)=xgrav(L)/NPOIN
  213.             enddo
  214.             IF (KCMETR.NE.0) THEN
  215.                DO 12 ININ=1,KCMETR.XIN(/1)
  216.                   KCMETR.XIN(ININ,IKPVIR)=XMET(ININ)/NPOIN
  217.  12            CONTINUE
  218.             ENDIF
  219.          ENDDO
  220.       ENDIF
  221. * Calcul du volume de la topologie initiale
  222. *  !! pour calculer le volume, il ne faut pas utiliser la métrique !
  223.       IELDEB=1
  224.       IELFIN=TRAVK.NVCOU
  225.       CALL VOMSI2(KEXTO,IELDEB,IELFIN,NKPVIR,XVTLOC,XVTLOS,XVTLOV)
  226.       IF (IERR.NE.0) RETURN
  227.       IF (IMPR.GE.4) THEN
  228.          WRITE(IOIMP,'(A,3(2X,A8,E10.3))') 'topv2 : vomsi2','XVTLOC='
  229.      $        ,XVTLOC,'XVTLOS=',XVTLOS,'XVTLOV=',XVTLOV
  230.       ENDIF
  231.       IF (NKPVIR.GT.0) THEN
  232.          if (ABS(XVTLOV).GE.XVTOL) THEN
  233.             WRITE(IOIMP,*) '!! Vomsi2 XVTLOV=',XVTLOV,' XVTOL=',XVTOL
  234.             WRITE(IOIMP,*) '          XVTLOC=',XVTLOC,' XVTLOS=',XVTLOS
  235.             write(ioimp,*) 'NVCOU=',TRAVK.NVCOU,'IELFIN=',IELFIN
  236.             write(ioimp,*) 'NKPVIR=',NKPVIR
  237. *            segprt,kexto
  238.             GOTO 7
  239.          endif
  240.       ENDIF
  241.       IF (IALGO.EQ.1) THEN
  242.          if (ABS(XVTLOC-XVTLOS).GE.XVTOL) THEN
  243.             WRITE(IOIMP,*) '!! Vomsi2 XVTLOC=',XVTLOC,' XVTLOS=',XVTLOS
  244.             WRITE(IOIMP,*) '          XVTLOV=',XVTLOV,' XVTOL=',XVTOL
  245.             write(ioimp,*) 'NVCOU=',TRAVK.NVCOU,'IELFIN=',IELFIN
  246.             write(ioimp,*) 'NKPVIR=',NKPVIR
  247. *            segprt,kexto
  248.             GOTO 7
  249.          endif
  250.       ENDIF
  251.  
  252. *      JCAND=0
  253. * Note : on saute les topologies locales ayant un volume nul
  254. * car on n'arrivera pas à les améliorer. En fait, rarement mais si.
  255. *!      IF (XVTLOC.LE.XVTOL) GOTO 7
  256. *!      impr=4
  257. *     On se demandait si on devait sauter les topologies avec des
  258. *     retournements...      XVTLOC.NE.XVTLOS
  259. *
  260. *     Génération des topologies candidates (stockage dans TRAVL: MCANS
  261. *     indexé par IDXCA)
  262. *     ICBES est le meilleur candidat par défaut (souvent=1, la topologie initiale)
  263. *     IPOPL2 est le candidat où on a créé un nouveau noeud (souvent le
  264. *     dernier, égal à travl.nccou)
  265. *
  266. *      write(ioimp,*) 'Entree topv3'
  267.       CALL TOPV3(TRAVK,KELEMX,IAJNO,TRAVL,INCMA,ISTMA,JNASCM,ICBES
  268.      $     ,IPOPL2,iveri,impr)
  269.       IF (IERR.NE.0) RETURN
  270. *      write(ioimp,*) 'Sortie topv3'
  271. *
  272. * On a les candidats dans ITCAND
  273. * On sélectionne ceux de volume minimum non nul
  274. *
  275.       JCAND=TRAVL.NCCOU
  276.       if (impr.ge.4) then
  277.          write(ioimp,'(A,I5)') 'topv2 : apres gen candidat n=',jcand
  278.       endif
  279. * Sortie anticipée s'il n'y a qu'un candidat
  280. *      WRITE(IOIMP,*) 'JCAND=',JCAND
  281.       IF (JCAND.EQ.1) GOTO 8
  282. *
  283.       XVMIN=XVTLOC
  284.       if (impr.ge.4) then
  285.          write(ioimp,'(A,E10.3)') 'topv2 : volume vise XVMIN=',XVMIN
  286.       endif
  287. *
  288. *      write(ioimp,*) 'Entree topv4'
  289.       CALL TOPV4(TRAVL,IALGO,XVMIN,XVTOL,NKPVIR,impr)
  290.       IF (IERR.NE.0) RETURN
  291. *      write(ioimp,*) 'Sortie topv4'
  292. * S'il n'en reste pas (ne doit pas se produire)
  293.       IF (TRAVL.NVOCOU.LE.0) THEN
  294.          write(ioimp,*) 'Plus de candidats apres topv4 ???'
  295.          goto 9999
  296. *         ICBES=0
  297. *         GOTO 666
  298.       ENDIF
  299. * S'il n'en reste qu'un, on peut sauter le calcul des qualités ??????
  300.       IF (TRAVL.NVOCOU.EQ.1) THEN
  301.          ICBES=LOKVOL.LECT(NVOCOU)
  302.          GOTO 8
  303.       ENDIF
  304. *
  305. *     Les candidats ayant le bon volume sont dans ITVOL
  306. *     On calcule les qualités de chaque élément des candidats et on ordonne.
  307. *
  308. *      write(ioimp,*) 'Entree topv5'
  309.       CALL TOPV5(TRAVL,XVTOL,IMET,IMOMET,XDENS,KCMETR,NKPVIR,jcritq
  310.      $     ,pcritq,qcritq)
  311.       IF (IERR.NE.0) RETURN
  312. *      write(ioimp,*) 'Sortie topv5'
  313.  
  314.  123  format (A10,I3,3X,A10,E10.3)
  315.  124  format ('(A10,10(2X,',I1,'(E10.3)))')
  316.  125  format ('(A,I3,3X,A10,20(2X,',I1,'(F7.3)))')
  317.  
  318.       if (impr.ge.4) then
  319. *     if (.true.) then
  320.          write(cform,125) 1
  321. *         write(ioimp,*) 'Apres calcul qualite candidats'
  322.          DO i=1,NVOCOU
  323. *            Write(ioimp,*) 'Candidat i=',i
  324.             ICAND=LOKVOL.LECT(i)
  325. *            write(ioimp,123) 'Cand i=',icand,'Vol=',XVMIN
  326.             NDQC=LNQUAL.LECT(ICAND)
  327.             IELDEB=LIDXCA.LECT(ICAND)
  328.             jeldeb=ieldeb-1
  329.             write(ioimp,cform) 'topv2 : candidat i=',icand,'qualites='
  330.      $           ,(lquals.prog(jeldeb+k),k=1,ndqc)
  331.          enddo
  332.       endif
  333. *
  334. * Calcul des meilleurs candidats par maximum lexical
  335. *
  336. *      write(ioimp,*) 'Entree topv6'
  337.       CALL TOPV6(TRAVL,QTOL,ICBES)
  338. *      write(ioimp,*) 'Sortie topv6'
  339. *
  340.       if (impr.ge.4) then
  341. *      if (.true.) then
  342.          write(ioimp,'(A,I5)')
  343.      $        'topv2 : apres maximum lexical, le meilleur est :',icbes
  344.       endif
  345.       LCQUAL=.TRUE.
  346.  8    CONTINUE
  347. *
  348.       if (.false.) then
  349.          IF (LCQUAL) THEN
  350.             write(cform,124) 1
  351. *     write(ioimp,*) 'cform=',cform
  352.          ENDIF
  353.          if (icbes.ne.1) THEN
  354.             ICAND=1
  355.             write(ioimp,123) 'Av i=',icand,'Vol=',XVTLOC
  356.             IF (LCQUAL) THEN
  357.                IELDEB=LIDXCA.LECT(ICAND)
  358.                NDQC=LNQUAL.LECT(ICAND)
  359.                jeldeb=ieldeb-1
  360.                write(ioimp,cform) 'Qualite=',(lquals.prog(jeldeb+k),k
  361.      $              =1,ndqc)
  362.             ENDIF
  363.             ICAND=ICBES
  364.             write(ioimp,123) 'Ap i=',icand,'Vol=',XVMIN
  365.             IF (LCQUAL) THEN
  366.                IELDEB=LIDXCA.LECT(ICAND)
  367.                NDQC=LNQUAL.LECT(ICAND)
  368.                jeldeb=ieldeb-1
  369.                write(ioimp,cform) 'Qualite=',(lquals.prog(jeldeb+k),k
  370.      $              =1,ndqc)
  371.             ENDIF
  372.          endif
  373.       endif
  374. *anc      KTBES2=KTBES
  375.       ICAND=ICBES
  376. *dbg      write(ioimp,*) 'icand,nccou=',icand,travl.nccou
  377.       IF (ICAND.NE.1) THEN
  378.          IELDEB=LIDXCA.LECT(ICAND)
  379.          IELFIN=LIDXCA.LECT(ICAND+1)-1
  380.          NBNN=LMCANS.NNCOU
  381.          NBELEM=IELFIN-IELDEB+1
  382.          NBSOUS=0
  383.          NBREF=0
  384. *            write(ioimp,*) 'icand,ieldeb,ielfin,nbnn,nbelem=',icand
  385. *     $           ,ieldeb,ielfin,nbnn,nbelem
  386.          SEGINI KTBES
  387.          KTBES.ITYPEL=LMCANS.ITYPEX
  388.          IDX=1
  389.          DO IEL=IELDEB,IELFIN
  390.             DO INO=1,NBNN
  391.                KTBES.NUM(INO,IDX)=LMCANS.NUMX(INO,IEL)
  392.             ENDDO
  393.             IDX=IDX+1
  394.          ENDDO
  395.       ENDIF
  396. *
  397.  666  CONTINUE
  398. * Si on n'a pas sélectionné le candidat avec le noeud supplémentaire
  399. * on peut enlever le noeud
  400. *ijob      IF (IJOB.EQ.1.OR.IJOB.EQ.2) THEN
  401.       IF (IAJNO.NE.0) THEN
  402.          IF (IPOPL2.NE.0.AND.ICBES.NE.IPOPL2) THEN
  403.             NPCOUN=TRAVK.NPCOU-1
  404. *     Remise à zéro : nécessaire uniquement pour la vérification ??
  405.             if (iveri.ge.1) then
  406.                IREF=(TRAVK.NPCOU-1)*IDIMP1
  407.                DO 111 I=1,IDIMP1
  408.                   KCOORD.XCOOR(IREF+I)=0.D0
  409.  111           CONTINUE
  410.                IF (KCMETR.NE.0) THEN
  411.                   DO 112 ININ=1,KCMETR.XIN(/1)
  412.                      KCMETR.XIN(ININ,TRAVK.NPCOU)=0.D0
  413.  112              CONTINUE
  414.                ENDIF
  415.             endif
  416.             TRAVK.NPCOU=NPCOUN
  417. * Désactivation temporaire car pas de ménage...
  418. *tmp            if (iveri.ge.2) then
  419. *tmp            call vetopi(travk,'topv2 : Apres reduction travk')
  420. *tmp               if (ierr.ne.0) return
  421. *tmp            endif
  422.          ELSE
  423. *            write(ioimp,*) 'Nouveau noeud cree ','IPOPL2,ICBES,NPCOU='
  424. *     $           ,IPOPL2,ICBES,TRAVK.NPCOU
  425. *            IREF=(TRAVK.NPCOU-1)*IDIMP1
  426. *            write(ioimp,*)  (kcoord.xcoor(iref+iii),iii=1,idimp1)
  427.          ENDIF
  428.       ENDIF
  429.  7    CONTINUE
  430.       IF (KEXTO.EQ.KTBES) JCAND=-JCAND
  431.       IF (KEXTO.EQ.KTBES) THEN
  432.          LCHTOP=.FALSE.
  433.       ELSE
  434.          LCHTOP=.TRUE.
  435.       ENDIF
  436. *      WRITE(IOIMP,*) 'sortie topv2'
  437.       RETURN
  438. *
  439. *
  440. *
  441.  9999 CONTINUE
  442.       MOTERR(1:8)='TOPV2   '
  443. * 349 2
  444. *Problème non prévu dans le s.p. %m1:8 contactez votre assistance
  445.       CALL ERREUR(349)
  446.       RETURN
  447. *
  448. * End of subroutine TOPV2
  449. *
  450.       END
  451.  
  452.  

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