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/01/09    21:16:01     12442          
  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.       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. *      WRITE(IOIMP,*) 'Vomsi2 XVTLOC=',XVTLOC,' XVTLOS=',XVTLOS,' XVTLOV
  228. *     $     =',XVTLOV
  229.       IF (NKPVIR.GT.0) THEN
  230.          if (ABS(XVTLOV).GE.XVTOL) THEN
  231. *            WRITE(IOIMP,*) '!! Vomsi2 XVTLOV=',XVTLOV,' XVTOL=',XVTOL
  232. *            WRITE(IOIMP,*) '          XVTLOC=',XVTLOC,' XVTLOS=',XVTLOS
  233. *            write(ioimp,*) 'NVCOU=',TRAVK.NVCOU,'IELFIN=',IELFIN
  234. *            write(ioimp,*) 'NKPVIR=',NKPVIR
  235. *            segprt,kexto
  236.             GOTO 7
  237.          endif
  238.       ENDIF
  239.       IF (IALGO.EQ.1) THEN
  240.          if (ABS(XVTLOC-XVTLOS).GE.XVTOL) THEN
  241. *            WRITE(IOIMP,*) '!! Vomsi2 XVTLOC=',XVTLOC,' XVTLOS=',XVTLOS
  242. *            WRITE(IOIMP,*) '          XVTLOV=',XVTLOV,' XVTOL=',XVTOL
  243. *            write(ioimp,*) 'NVCOU=',TRAVK.NVCOU,'IELFIN=',IELFIN
  244. *            write(ioimp,*) 'NKPVIR=',NKPVIR
  245. *            segprt,kexto
  246.             GOTO 7
  247.          endif
  248.       ENDIF
  249.  
  250. *      JCAND=0
  251. * Note : on saute les topologies locales ayant un volume nul
  252. * car on n'arrivera pas à les améliorer. En fait, rarement mais si.
  253. *!      IF (XVTLOC.LE.XVTOL) GOTO 7
  254. *!      impr=4
  255. *     On se demandait si on devait sauter les topologies avec des
  256. *     retournements...      XVTLOC.NE.XVTLOS
  257. *
  258. *     Génération des topologies candidates (stockage dans TRAVL: MCANS
  259. *     indexé par IDXCA)
  260. *     ICBES est le meilleur candidat par défaut (souvent=1, la topologie initiale)
  261. *     IPOPL2 est le candidat où on a créé un nouveau noeud (souvent le
  262. *     dernier, égal à travl.nccou)
  263. *
  264. *      write(ioimp,*) 'Entree topv3'
  265.       CALL TOPV3(TRAVK,KELEMX,IAJNO,TRAVL,INCMA,ISTMA,JNASCM,ICBES
  266.      $     ,IPOPL2,iveri,impr)
  267.       IF (IERR.NE.0) RETURN
  268. *      write(ioimp,*) 'Sortie topv3'
  269. *
  270. * On a les candidats dans ITCAND
  271. * On sélectionne ceux de volume minimum non nul
  272. *
  273.       JCAND=TRAVL.NCCOU
  274.       if (impr.ge.4) then
  275.          write(ioimp,*) 'Apres gen candidat n=',jcand
  276.       endif
  277. * Sortie anticipée s'il n'y a qu'un candidat
  278. *      WRITE(IOIMP,*) 'JCAND=',JCAND
  279.       IF (JCAND.EQ.1) GOTO 8
  280. *
  281.       XVMIN=XVTLOC
  282.       if (impr.ge.4) then
  283.          write(ioimp,*) 'volume vise XVMIN=',XVMIN
  284.       endif
  285. *
  286. *      write(ioimp,*) 'Entree topv4'
  287.       CALL TOPV4(TRAVL,IALGO,XVMIN,XVTOL,NKPVIR,impr)
  288.       IF (IERR.NE.0) RETURN
  289. *      write(ioimp,*) 'Sortie topv4'
  290. * S'il n'en reste pas (ne doit pas se produire)
  291.       IF (TRAVL.NVOCOU.LE.0) THEN
  292.          write(ioimp,*) 'Plus de candidats apres topv4 ???'
  293.          goto 9999
  294. *         ICBES=0
  295. *         GOTO 666
  296.       ENDIF
  297. * S'il n'en reste qu'un, on peut sauter le calcul des qualités ??????
  298.       IF (TRAVL.NVOCOU.EQ.1) THEN
  299.          ICBES=LOKVOL.LECT(NVOCOU)
  300.          GOTO 8
  301.       ENDIF
  302. *
  303. *     Les candidats ayant le bon volume sont dans ITVOL
  304. *     On calcule les qualités de chaque élément des candidats et on ordonne.
  305. *
  306. *      write(ioimp,*) 'Entree topv5'
  307.       CALL TOPV5(TRAVL,XVTOL,IMET,IMOMET,XDENS,KCMETR,NKPVIR,ISTRID)
  308.       IF (IERR.NE.0) RETURN
  309. *      write(ioimp,*) 'Sortie topv5'
  310.  
  311.  
  312.       if (impr.ge.4) then
  313. *      if (.true.) then
  314.          write(cform,124) ISTRID
  315.          write(ioimp,*) 'Apres calcul qualite candidats'
  316.          DO i=1,NVOCOU
  317. *            Write(ioimp,*) 'Candidat i=',i
  318.             ICAND=LOKVOL.LECT(i)
  319.             write(ioimp,123) 'Cand i=',icand,'Vol=',XVMIN
  320.             NDQC=LNQUAL.LECT(ICAND)
  321.             IELDEB=LIDXCA.LECT(ICAND)
  322.             jeldeb=(ieldeb-1)*istrid
  323. *            write(ioimp,*) 'i,icand,ndqc,ieldeb,jeldeb=',i,icand,ndqc
  324. *     $           ,ieldeb,jeldeb
  325.             write(ioimp,cform) 'Qualite=',(lquals.prog(jeldeb+k),k
  326.      $           =1,ndqc*istrid)
  327. *            write(ioimp,cform) 'Cand i=',i,(lquals.prog(jeldeb+k),k=1
  328. *     $           ,ndqc*istrid)
  329. *            do iii=ieldeb,ielfin
  330. *               jjj=(iii-1)*istrid
  331. *               Write(ioimp,*) (lquals.prog(jjj+k),k=1,istrid)
  332. *     $           ,lnqual.lect(icand))
  333. *            enddo
  334.          enddo
  335.       endif
  336. *
  337. * Calcul des meilleurs candidats par maximum lexical
  338. *
  339. *      write(ioimp,*) 'Entree topv6'
  340.       CALL TOPV6(TRAVL,QTOL,ISTRID,ICBES)
  341. *      write(ioimp,*) 'Sortie topv6'
  342. *
  343.       if (impr.ge.4) then
  344. *      if (.true.) then
  345.          write(ioimp,*) 'Apres maximum lexical, le meilleur est :'
  346.      $        ,icbes
  347.       endif
  348.       LCQUAL=.TRUE.
  349.  8    CONTINUE
  350. *
  351.  123     format (A10,I3,3X,A10,E10.3)
  352.  124     format ('(A10,10(2X,',I1,'(E10.3)))')
  353.       if (.false.) then
  354.          IF (LCQUAL) THEN
  355.             write(cform,124) ISTRID
  356. *     write(ioimp,*) 'cform=',cform
  357.          ENDIF
  358.          if (icbes.ne.1) THEN
  359.             ICAND=1
  360.             write(ioimp,123) 'Av i=',icand,'Vol=',XVTLOC
  361.             IF (LCQUAL) THEN
  362.                IELDEB=LIDXCA.LECT(ICAND)
  363.                NDQC=LNQUAL.LECT(ICAND)
  364.                jeldeb=(ieldeb-1)*istrid
  365.                write(ioimp,cform) 'Qualite=',(lquals.prog(jeldeb+k),k
  366.      $              =1,ndqc*istrid)
  367.             ENDIF
  368.             ICAND=ICBES
  369.             write(ioimp,123) 'Ap i=',icand,'Vol=',XVMIN
  370.             IF (LCQUAL) THEN
  371.                IELDEB=LIDXCA.LECT(ICAND)
  372.                NDQC=LNQUAL.LECT(ICAND)
  373.                jeldeb=(ieldeb-1)*istrid
  374.                write(ioimp,cform) 'Qualite=',(lquals.prog(jeldeb+k),k
  375.      $              =1,ndqc*istrid)
  376.             ENDIF
  377.          endif
  378.       endif
  379. *anc      KTBES2=KTBES
  380.       ICAND=ICBES
  381. *dbg      write(ioimp,*) 'icand,nccou=',icand,travl.nccou
  382.       IF (ICAND.NE.1) THEN
  383.          IELDEB=LIDXCA.LECT(ICAND)
  384.          IELFIN=LIDXCA.LECT(ICAND+1)-1
  385.          NBNN=LMCANS.NNCOU
  386.          NBELEM=IELFIN-IELDEB+1
  387.          NBSOUS=0
  388.          NBREF=0
  389. *            write(ioimp,*) 'icand,ieldeb,ielfin,nbnn,nbelem=',icand
  390. *     $           ,ieldeb,ielfin,nbnn,nbelem
  391.          SEGINI KTBES
  392.          KTBES.ITYPEL=LMCANS.ITYPEX
  393.          IDX=1
  394.          DO IEL=IELDEB,IELFIN
  395.             DO INO=1,NBNN
  396.                KTBES.NUM(INO,IDX)=LMCANS.NUMX(INO,IEL)
  397.             ENDDO
  398.             IDX=IDX+1
  399.          ENDDO
  400.       ENDIF
  401. *
  402.  666  CONTINUE
  403. * Si on n'a pas sélectionné le candidat avec le noeud supplémentaire
  404. * on peut enlever le noeud
  405. *ijob      IF (IJOB.EQ.1.OR.IJOB.EQ.2) THEN
  406.       IF (IAJNO.NE.0) THEN
  407.          IF (IPOPL2.NE.0.AND.ICBES.NE.IPOPL2) THEN
  408.             NPCOUN=TRAVK.NPCOU-1
  409. *     Remise à zéro : nécessaire uniquement pour la vérification ??
  410.             if (iveri.ge.1) then
  411.                IREF=(TRAVK.NPCOU-1)*IDIMP1
  412.                DO 111 I=1,IDIMP1
  413.                   KCOORD.XCOOR(IREF+I)=0.D0
  414.  111           CONTINUE
  415.                IF (KCMETR.NE.0) THEN
  416.                   DO 112 ININ=1,KCMETR.XIN(/1)
  417.                      KCMETR.XIN(ININ,TRAVK.NPCOU)=0.D0
  418.  112              CONTINUE
  419.                ENDIF
  420.             endif
  421.             TRAVK.NPCOU=NPCOUN
  422. * Désactivation temporaire car pas de ménage...
  423. *tmp            if (iveri.ge.2) then
  424. *tmp            call vetopi(travk,'topv2 : Apres reduction travk')
  425. *tmp               if (ierr.ne.0) return
  426. *tmp            endif
  427.          ELSE
  428. *            write(ioimp,*) 'Nouveau noeud cree ','IPOPL2,ICBES,NPCOU='
  429. *     $           ,IPOPL2,ICBES,TRAVK.NPCOU
  430. *            IREF=(TRAVK.NPCOU-1)*IDIMP1
  431. *            write(ioimp,*)  (kcoord.xcoor(iref+iii),iii=1,idimp1)
  432.          ENDIF
  433.       ENDIF
  434.  7    CONTINUE
  435.       IF (KEXTO.EQ.KTBES) JCAND=-JCAND
  436.       IF (KEXTO.EQ.KTBES) THEN
  437.          LCHTOP=.FALSE.
  438.       ELSE
  439.          LCHTOP=.TRUE.
  440.       ENDIF
  441. *      WRITE(IOIMP,*) 'sortie topv2'
  442.       RETURN
  443. *
  444. *
  445. *
  446.  9999 CONTINUE
  447.       MOTERR(1:8)='TOPV2   '
  448. * 349 2
  449. *Problème non prévu dans le s.p. %m1:8 contactez votre assistance
  450.       CALL ERREUR(349)
  451.       RETURN
  452. *
  453. * End of subroutine TOPV2
  454. *
  455.       END
  456.  
  457.  

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