Télécharger la plaquette Cast3M


fr en

Télécharger envvo3.eso

Retour à la liste

Numérotation des lignes : 

envvo3

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

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