Télécharger la plaquette Cast3M


fr en

Télécharger envel1.eso

Retour à la liste

Numérotation des lignes : 

envel1

  1. C ENVEL1    SOURCE    CB215821  19/08/20    21:17:02     10287          
  2. *   copier sur envvol avec gestion du chamelem de valeurs associes
  3. *   utilise dans trac cham
  4. C
  5. C  SG 2016/07/20 Programmation comme faced2, envvo2 pour gerer les faces TRI7/QUA9
  6. *
  7.       SUBROUTINE ENVEL1(MELEME,MELRES,mcoup)
  8.       IMPLICIT INTEGER(I-N)
  9.       IMPLICIT REAL*8 (A-H,O-Z)
  10.  
  11.  
  12.  
  13. -INC PPARAM
  14. -INC CCOPTIO
  15. -INC SMELEME
  16. -INC CCGEOME
  17. -INC SMCHAML
  18. *
  19. *     Type de faces prises en compte: T3, Q4, T6, Q8, T7, Q9
  20. *     Numero dans KDFAC               1   2   3   4    7   8
  21. * Pour ne pas se prendre la tête,  on numerote pareil que KDFAC
  22. * Pour les 5 (non utilisé), 6 (polygone) et >8, ca restera à 0
  23. *  NTYFAC=20, codé en dur dans CCGEOME pour KDFAC
  24.       PARAMETER (NTYFAC=20)
  25. * Nb de faces de chaque type, sert également de compteur
  26.       SEGMENT NBFAC(NTYFAC)
  27. * Tableau d'index de début fin dans les tableaux ???(NFAC)
  28.       SEGMENT IDXFAC(NTYFAC+1)
  29. * Pointeurs sur des segments MELEME et MELVAL par type de face
  30.       SEGMENT IPTFAC(NTYFAC)
  31.       SEGMENT MLVFAC(NTYFAC)
  32. * Un segment pointant sur les IFACI et les XFACI
  33.       SEGMENT IPOFAC(2,NTYFAC)
  34. *     Segment IFACI contenant les noeuds, la couleur et si la face d'un
  35. *     type donné est vue ou non
  36.       SEGMENT IFACI(NNODE+2,NFACI)
  37. *     Segment XFACI contenant les coordonnees noeuds, la couleur et si la face d'un
  38. *     type donné est vue ou non
  39.       SEGMENT XFACI(NNODE,NFACI)
  40. *
  41.       SEGMENT IPPOL(NTPOL)
  42.       SEGMENT IPREF(NTPOL)
  43.       SEGMENT NTFAC(NFAC)
  44.       SEGMENT KFAK(NFAC)
  45.       SEGMENT NAUX(max(2,NFAC))
  46. *SG
  47. *     Logique loquaf : pour les faces TRI7 et QUA9, normalement, le
  48. *     dernier noeud de la face est unique à la face : il peut donc
  49. *     servir de clé de hachage et on peut éviter de vérifier l'égaliteé
  50. *     de tous les autres noeuds lorsque l'on teste l'égalité des faces.
  51. *     C'est ce qu'on fait si loquaf=vrai.
  52. *
  53.       LOGICAL LOQUAF,LOPT
  54.       PARAMETER (LOQUAF=.TRUE.)
  55. * Pour chaque face dans KDFAC, le numéro d'élément associé
  56. * Ne se trouve pas dans CCGEOME, etonnant
  57.       INTEGER ITYEL(NTYFAC)
  58. *     T3, Q4, T6, Q8, ? , POLY, T7, Q9
  59.       DATA ITYEL/4,8,6,10,0,0,7,11,12*0/
  60.  
  61. *dbg      write(ioimp,*) 'coucou envel1'
  62.       n2ptel=0
  63.       n2el=0
  64.       SEGACT MELEME
  65.  
  66. c     on compte le nombre d elements dont les faces sont de type 1 2 3 4
  67. c     7 8 dans NBFAC
  68.       SEGINI NBFAC
  69.       NTPOL=0
  70.       IPT1=MELEME
  71.       SEGACT MELEME
  72.       nbsour=lisous(/1)
  73.       if (mcoup.ne.0) nbsour=nbsour-1
  74.       DO 10 IOB=1,nbsour
  75.          IPT1=LISOUS(IOB)
  76.          SEGACT IPT1
  77.          NBELEM=IPT1.NUM(/2)
  78.          ILTEL=LTEL(1,IPT1.ITYPEL)
  79.          IF (ILTEL.EQ.0) GOTO 12
  80.          ILTAD=LTEL(2,IPT1.ITYPEL)
  81.          DO 13 IF=1,ILTEL
  82.             IFT=LDEL(1,ILTAD+IF-1)
  83.             IF (IFT.EQ.6) THEN 
  84.                NTPOL=NTPOL+1
  85.             ELSE
  86.                NBFAC(IFT)=NBFAC(IFT)+NBELEM
  87.             ENDIF
  88.  13      CONTINUE
  89.  12      CONTINUE
  90.  10   CONTINUE
  91.  
  92. c==== CREATION DES FACES ==============================================
  93. * Initialisation des IFACI,XFACI
  94.       SEGINI IPOFAC
  95.       DO ITYFAC=1,NTYFAC
  96.          NNODE=KDFAC(1,ITYFAC)
  97.          IF (NNODE.GT.0) THEN
  98.             NFACI=NBFAC(ITYFAC)
  99.             SEGINI IFACI
  100.             IPOFAC(1,ITYFAC)=IFACI
  101.             SEGINI XFACI
  102.             IPOFAC(2,ITYFAC)=XFACI
  103.          ENDIF
  104.       ENDDO
  105.       SEGINI IPPOL,IPREF
  106. c     NBFAC sert maintenant de compteur
  107.       DO ITYFAC=1,NTYFAC
  108.          NBFAC(ITYFAC)=0
  109.       ENDDO
  110.       NTPOL=0
  111.       DO 50 IOB=1,nbsour
  112.          IPT1=LISOUS(IOB)
  113. *     si objet en double on saute
  114.          do 51 io2=1,iob-1
  115.             if (ipt1.eq.lisous(io2)) goto 50
  116.  51      continue
  117.          SEGACT IPT1
  118.          IELIM=1
  119.          IF (KSURF(IPT1.ITYPEL).EQ.IPT1.ITYPEL) THEN
  120. *  face non eliminable (pas un volume)
  121.             IELIM=0
  122.          ENDIF
  123.          melval=lisref(iob)
  124.          if (melval.eq.0) goto 50
  125.          segact melval
  126.          lval=velche(/1)
  127.          ival=velche(/2)
  128.          NBELEM=IPT1.NUM(/2)
  129.          ILTEL=LTEL(1,IPT1.ITYPEL)
  130.          IF (ILTEL.EQ.0) GOTO 52
  131.          ILTAD=LTEL(2,IPT1.ITYPEL)
  132.          DO 60 IF=1,ILTEL
  133.             ITYFAC=LDEL(1,ILTAD+IF-1)
  134.             IAD=LDEL(2,ILTAD+IF-1)
  135.             NNODE=KDFAC(1,ITYFAC)
  136.             IF (NNODE.GT.0) THEN
  137.                IFACI=IPOFAC(1,ITYFAC)
  138.                XFACI=IPOFAC(2,ITYFAC)
  139.                DO 66 IEL=1,NBELEM
  140.                   ielr=min(ival,iel)
  141.                   NBFAC(ITYFAC)=NBFAC(ITYFAC)+1
  142.                   j=NBFAC(ITYFAC)
  143.                   DO i=1,NNODE
  144.                      IFACI(i,j)=IPT1.NUM(LFAC(IAD+i-1),IEL)
  145.                      XFACI(i,j)=velche(min(lval,LFAC(IAD+i-1)),ielr)
  146.                   ENDDO
  147.                   IFACI(NNODE+1,j)=IPT1.ICOLOR(IEL)
  148.                   IFACI(NNODE+2,j)=IELIM
  149.  66            CONTINUE
  150.             ENDIF
  151. *     Avant ce if était après le 52 CONTINUE mais alors ITYFAC etait
  152. *     potentiellement non initialise
  153.             IF (ITYFAC.EQ.6) THEN
  154. C     Polygone
  155.                NTPOL = NTPOL+1
  156.                IPPOL(NTPOL)= IPT1
  157.                SEGINI,MELVA1 = MELVAL
  158.                IPREF(NTPOL) = MELVA1
  159.             ENDIF
  160.  60      CONTINUE
  161.  52      CONTINUE
  162.  50   CONTINUE
  163. C  IF FAUT MAINTENANT RETASSER ET CLASSER LES TABLEAUX DES FACES
  164. C  PROBLEME ELLES NE SONT PAS FORCEMENT DECRITES DE LA MEME FACON
  165. C     SG 20160712 NTFAC sert de cle de hachage, elle est égale à la
  166. C     somme des numeros de noeuds de la face
  167.       NFAC=0
  168.       SEGINI IDXFAC
  169.       IDXFAC(1)=NFAC+1
  170.       DO ITYFAC=1,NTYFAC
  171.          NFAC=NFAC+NBFAC(ITYFAC)
  172.          IDXFAC(ITYFAC+1)=NFAC+1
  173. *         write(ioimp,*) 'ityfac=',ityfac,' nbfac=',NBFAC(ITYFAC)
  174.       ENDDO
  175.       SEGINI NTFAC,KFAK
  176.       IFAC=0
  177.       DO ITYFAC=1,NTYFAC
  178.          NNODE=KDFAC(1,ITYFAC)
  179.          IF (NNODE.GT.0) THEN
  180.             LOPT=LOQUAF.AND.(ITYFAC.EQ.7.OR.ITYFAC.EQ.8)
  181.             IFACI=IPOFAC(1,ITYFAC)
  182.             DO I=1,NBFAC(ITYFAC)
  183.                IFAC=IFAC+1
  184.                IF (LOPT) THEN
  185.                   NTFAC(IFAC)=IFACI(NNODE,I)
  186.                ELSE
  187.                   DO J=1,NNODE
  188.                      NTFAC(IFAC)=NTFAC(IFAC)+IFACI(J,I)
  189.                   ENDDO
  190.                ENDIF
  191.                KFAK(IFAC)=I
  192.             ENDDO
  193.          ENDIF
  194.       ENDDO
  195. C  IL N'Y A PLUS QU'A TRIER ET RETASSER KFAK SUIVANT NTFAC
  196.       SEGINI NAUX
  197.       DO 300 ITYFAC=1,NTYFAC
  198.          IDEB=IDXFAC(ITYFAC)
  199.          IFIN=IDXFAC(ITYFAC+1)-1
  200.          IF (IFIN.LE.IDEB) GOTO 300
  201.          NAUX(1)=IDEB
  202.          NAUX(2)=IFIN
  203.          IZ=2
  204.  208     IZ=IZ-1
  205.          IF (IZ.LE.0) GOTO 209
  206.          IPB=NAUX(IZ*2-1)
  207.          IPH=NAUX(IZ*2)
  208.          IF(IPB.GE.IPH) GOTO 208
  209.          JPB=IPB-1
  210.          JPH=IPH+1
  211. C  CALCUL DU PIVOT
  212.          NPV=0
  213. *     DO 207 J=IPB,IPH
  214. *     NPV=NPV+NTFAC(J)
  215. *207  CONTINUE
  216. *     NPV=NPV/(IPH-IPB+1)
  217.          NPV=(NTFAC(IPB)+NTFAC(IPH))/2
  218.  242     JPB=JPB+1
  219.          IF (JPH.EQ.JPB) GOTO 245
  220.          IF (NTFAC(JPB).LE.NPV) GOTO 243
  221.          GOTO 242
  222.  243     JPH=JPH-1
  223.          IF (JPH.EQ.JPB) GOTO 245
  224.          IF (NTFAC(JPH).GE.NPV) GOTO 244
  225.          GOTO 243
  226.  244     IAUX=KFAK(JPB)
  227.          KFAK(JPB)=KFAK(JPH)
  228.          KFAK(JPH)=IAUX
  229.          NTAUX=NTFAC(JPB)
  230.          NTFAC(JPB)=NTFAC(JPH)
  231.          NTFAC(JPH)=NTAUX
  232.          GOTO 242
  233.  245     IF (JPB.GE.IPB) GOTO 247
  234.          JPB=JPB+1
  235.          JPH=JPH+2
  236.          GOTO 248
  237.  247     IF (JPH.LE.IPH) GOTO 249
  238.          JPB=JPB-2
  239.          JPH=JPH-1
  240.          GOTO 248
  241.  249     IF (NTFAC(JPB).LE.NPV) GOTO 250
  242.          IF (JPH.EQ.IPH) GOTO 251
  243.  252     JPH=JPH+1
  244.          GOTO 248
  245.  250     IF (JPB.EQ.IPB) GOTO 252
  246.  251     JPB=JPB-1
  247.  248     IF (JPB.EQ.IPB) GOTO 253
  248.          NAUX(2*IZ)=JPB
  249.          IZ=IZ+1
  250.  253     IF (JPH.EQ.IPH) GOTO 208
  251.          NAUX(2*IZ)=IPH
  252.          NAUX(2*IZ-1)=JPH
  253.          IZ=IZ+1
  254.          GOTO 208
  255.  209     CONTINUE
  256.  300  CONTINUE
  257. C  LES FACES SONT CLASSEES DANS KFAK IL FAUT ELIMINER LES FACES EN DOUBL
  258. C  ELLES SONT PAR TYPE LES UNES DERRIERES LES AUTRES LES PLUS HAUTES
  259. C  DEVANT
  260.       IF (IIMPI.NE.0) WRITE (IOIMP,9111) (KFAK(I),NTFAC(I),I=1,NFAC)
  261.  9111 FORMAT(5(2X,2I6))
  262.       DO 400 ITYFAC=1,NTYFAC
  263.          IDEB=IDXFAC(ITYFAC)
  264.          IFIN=IDXFAC(ITYFAC+1)-1
  265.          IF (IFIN.LE.IDEB) GOTO 400
  266.          NNODE=KDFAC(1,ITYFAC)
  267. * A cette etape on doit avoir nnode.gt.0
  268.          IF (NNODE.LE.0) THEN
  269.             CALL ERREUR(5)
  270.             RETURN
  271.          ENDIF
  272.          LOPT=LOQUAF.AND.(ITYFAC.EQ.7.OR.ITYFAC.EQ.8)
  273.          IFACI=IPOFAC(1,ITYFAC)
  274. *
  275.          IFINM=IFIN-1
  276.          DO 450 I1=IDEB,IFINM
  277.             NTI1=NTFAC(I1)
  278.             IF (NTI1.EQ.0) GOTO 450
  279.             IDEB1=I1+1
  280.             DO 460 I2=IDEB1,IFIN
  281.                NTI2=NTFAC(I2)
  282.                IF (NTI2.EQ.0) GOTO 460
  283.                IF (NTI2.NE.NTI1) GOTO 450
  284.                IR1=KFAK(I1)
  285.                IR2=KFAK(I2)
  286.                IF (IFACI(NNODE+2,IR1).EQ.0) GOTO 460
  287.                IF (IFACI(NNODE+2,IR2).EQ.0) GOTO 460
  288.                IF (.NOT.LOPT) THEN
  289.                   DO 471 J1=1,NNODE
  290.                      INU=IFACI(J1,IR1)
  291.                      DO 472 J2=1,NNODE
  292.                         IF (INU.EQ.IFACI(J2,IR2)) GOTO 471
  293.  472                 CONTINUE
  294.                      GOTO 460
  295.  471              CONTINUE
  296.                ENDIF
  297. C  DEUX FACES EGALES ON LES SUPPRIMENT
  298.                NTFAC(I1)=0
  299.                NTFAC(I2)=0
  300.                GOTO 450
  301.  460        CONTINUE
  302.  450     CONTINUE
  303.  400  CONTINUE
  304. * 
  305.       IF (IIMPI.NE.0) WRITE (IOIMP,9111) (KFAK(I),NTFAC(I),I=1,NFAC)
  306.  
  307.       SEGINI IPTFAC,MLVFAC
  308.       NBSOUS=0
  309.       NBREF=0
  310.       NBSOU2=0
  311.       DO 600 ITYFAC=1,NTYFAC
  312.          IDEB=IDXFAC(ITYFAC)
  313.          IFIN=IDXFAC(ITYFAC+1)-1
  314. *         write(ioimp,*) 'ityfac=',ityfac,' ideb=',ideb,' ifin=',ifin
  315.          IF (IFIN.LT.IDEB) GOTO 600
  316.          NNODE=KDFAC(1,ITYFAC)
  317. * A cette etape on doit avoir nnode.gt.0
  318.          IF (NNODE.LE.0) THEN
  319.             CALL ERREUR(5)
  320.             RETURN
  321.          ENDIF
  322.          IFACI=IPOFAC(1,ITYFAC)
  323.          XFACI=IPOFAC(2,ITYFAC)
  324.          NBELEM=0
  325.          DO 611 I=IDEB,IFIN
  326.             IF (NTFAC(I).NE.0) NBELEM=NBELEM+1
  327.  611     CONTINUE
  328. *         write(ioimp,*) 'nbelem=',nbelem,' nnode=',nnode
  329.          IF (NBELEM.EQ.0) GOTO 600
  330.          NBSOU2=NBSOU2+1
  331.          NBNN=NNODE
  332.          SEGINI IPT1
  333.          IPT1.ITYPEL=ITYEL(ITYFAC)
  334.          n1ptel=nnode
  335.          n1el=nbelem
  336.          segini melva1
  337.          JAUX=0
  338.          DO 612 J=IDEB,IFIN
  339.             IF (NTFAC(J).EQ.0) GOTO 612
  340.             JAUX=JAUX+1
  341.             IPT1.ICOLOR(JAUX)=IFACI(NNODE+1,KFAK(J))
  342.             DO 613 I=1,NBNN
  343.                IPT1.NUM(I,JAUX)=IFACI(I,KFAK(J))
  344.                melva1.velche(I,JAUX)=XFACI(I,KFAK(J))
  345.  613        CONTINUE
  346.  612     CONTINUE
  347.          IPTFAC(ITYFAC)=IPT1
  348. *         write(ioimp,*) 'ipt1=',ipt1 
  349.          MLVFAC(ITYFAC)=melva1
  350.  600  CONTINUE
  351. *  on rajoute les points et les segments qui pouvaient  etre dans le
  352. *  maillage initial
  353.       ipt5=0
  354.       segact meleme
  355.       ipt6=meleme
  356.       do 710 io=1,max(1,nbsour)
  357.          if (nbsour.ne.0) ipt6=lisous(io)
  358.          segact ipt6
  359.          if (ipt6.itypel.le.3) then
  360.             nbsou2=nbsou2+1
  361.             ipt5=ipt6
  362.          endif
  363.  710  continue
  364. *      write(ioimp,*) 'nbsou2=',nbsou2
  365.       NBSOUS=NBSOU2+NTPOL
  366.       if (mcoup.ne.0) nbsous=nbsous+1
  367.       IF (NBSOUS.EQ.0) CALL ERREUR(26)
  368.       IF (IERR.NE.0) RETURN
  369.       NBREF=nbsous
  370.       NBNN=0
  371.       NBELEM=0
  372.       SEGINI IPT5
  373.       I=0
  374.       DO ITYFAC=1,NTYFAC
  375.          IPT1=IPTFAC(ITYFAC)
  376.          melva1=MLVFAC(ITYFAC)
  377.          IF (IPT1.NE.0) THEN
  378.             if (melva1.eq.0) then 
  379.                call erreur(5)
  380.                return
  381.             endif
  382.             I=I+1
  383.             IPT5.LISOUS(I)=IPT1
  384.             IPT5.LISref(I)=melva1
  385.          ENDIF
  386.       ENDDO
  387.       segact meleme
  388.       ipt1=meleme
  389.       do 711 io=1,max(1,nbsour)
  390.          if (nbsour.ne.0) ipt1=lisous(io)
  391.          segact ipt1
  392.          if (ipt1.itypel.le.3) then
  393.             I=I+1
  394.             IPT5.LISOUS(I)=IPT1
  395.             IPT5.LISref(I)=lisref(io)
  396.          endif
  397.  711  continue
  398.       DO 720, IO = 1, NTPOL
  399.          I= I+1
  400.          IPT5.LISOUS(I) = IPPOL(IO)
  401.          IPT5.LISREF(I) = IPREF(IO)
  402.  720  CONTINUE
  403.       if (mcoup.ne.0) then
  404.          I= I+1
  405.          IPT5.LISOUS(I) = lisous(nbsour+1)
  406.          IPT5.LISREF(I) = lisref(nbsour+1)
  407.       endif
  408.       melres=ipt5
  409.  
  410.       SEGSUP IPTFAC,MLVFAC,NAUX,NTFAC,KFAK,IDXFAC,IPPOL,IPREF
  411.       DO ITYFAC=1,NTYFAC
  412.          IFACI=IPOFAC(1,ITYFAC)
  413.          IF (IFACI.NE.0) THEN
  414.             SEGSUP IFACI
  415.          ENDIF
  416.          XFACI=IPOFAC(2,ITYFAC)
  417.          IF (XFACI.NE.0) THEN
  418.             SEGSUP XFACI
  419.          ENDIF
  420.       ENDDO
  421.       SEGSUP IPOFAC,NBFAC
  422.       END
  423.  
  424.  
  425.  

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