Télécharger la plaquette Cast3M


fr en

Télécharger faced2.eso

Retour à la liste

Numérotation des lignes : 

faced2

  1. C FACED2    SOURCE    PASCAL    21/02/11    21:15:11     10890          
  2. C TRACE D'ISOVALEUR EN COMMENCANT PAR CELLES DE DERRIERE
  3. C
  4. C  SG 2016/07/18 Programmation comme faced, envvo2 pour gerer les faces TRI7/QUA9
  5. C
  6.       SUBROUTINE FACED2(MELEME,XPROJ,ICPR,VCHC,VCPCHA,PTI,NISO,IVU,
  7.      # MCOUP,mcham,BLOK)
  8. *  si mchaml est non nul on trace un chamelem aux noeuds
  9. *  sinon c'est un champoint
  10.       IMPLICIT INTEGER(I-N)
  11. -INC PPARAM
  12. -INC CCOPTIO
  13. -INC SMELEME
  14. -INC CCGEOME
  15. -INC SMCHAML
  16.       DIMENSION XX(13),YY(13),ZZ(13),VV(13)
  17.       DIMENSION XR(3),YR(3),ZR(3),VR(3)
  18.       SEGMENT XPROJ(3,1)
  19.       SEGMENT ICPR(1)
  20.       SEGMENT VCPCHA(nbpts)
  21.       SEGMENT IPOI1(1,NPOI1)
  22.       SEGMENT ISEG2(2,NSEG2)
  23.       SEGMENT ISEG3(3,NSEG3)
  24.       SEGMENT XPOI1(1,NPOI1)
  25.       SEGMENT XSEG2(2,NSEG2)
  26.       SEGMENT XSEG3(3,NSEG3)
  27.       logical lmvid
  28. *
  29. *     Type de faces prises en compte: T3, Q4, T6, Q8, POLY, T7, Q9
  30. *     Numero dans KDFAC               1   2   3   4     6    7   8
  31. * Pour ne pas se prendre la tête,  on numerote pareil que KDFAC
  32. * Pour les 5 (non utilisé), 6 (polygone) et >8, ca restera à 0
  33. *  NTYFAC=20, codé en dur dans CCGEOME pour KDFAC
  34.       PARAMETER (NTYFAC=20)
  35. * Nb de faces de chaque type, sert également de compteur
  36.       SEGMENT NBFAC(NTYFAC)
  37. * Un segment pointant sur les IFACI et les XFACI
  38.       SEGMENT IPOFAC(2,NTYFAC)
  39. *     Segment IFACI contenant les noeuds et si la face d'un
  40. *     type donné est vue ou non
  41.       SEGMENT IFACI(NNODE+1,NFACI)
  42. *     Segment XFACI contenant les coordonnees noeuds, la couleur et si la face d'un
  43. *     type donné est vue ou non
  44.       SEGMENT XFACI(NNODE,NFACI)
  45. *     Nombre de noeuds max pour les polygones
  46.       PARAMETER (NNOMAX=14)
  47. *
  48.       REAL BLOK
  49. *
  50.       SEGMENT NSOMP(NFACP)
  51.       SEGMENT TFAC(NFAC)
  52.       SEGMENT KFAK(NFAC)
  53.       SEGMENT NAUX(max(2,NFAC))
  54.       SEGMENT IVU(0)
  55.       SEGMENT MCOUP(0)
  56.       REAL VCHC
  57.       DIMENSION VCHC(*)
  58. *     Poids pour le calcul de la valeur centrale 6 valeurs TRI6 8
  59. *     valeurs QUA8. La somme des poids est egale a 1
  60.       PARAMETER (XUS3=1./3.,XMUS4=-0.25,XUS2=0.5)
  61.       REAL XPOIDS(14)
  62. *                    TRI6                
  63.       DATA XPOIDS/0.,XUS3,0.,XUS3,0.,XUS3,
  64. *                    QUA8
  65.      $     XMUS4,XUS2, XMUS4,XUS2, XMUS4,XUS2, XMUS4,XUS2/
  66. *      
  67. *dbg      write(ioimp,*) 'coucou faced2 mcham=',mcham
  68.       SEGACT MELEME
  69.  
  70. *     ipt1=lisous(1)
  71. *     segact ipt1
  72. *     SEGACT XPROJ,ICPR
  73.       MELSAU=MELEME
  74.       IPT1=MELEME
  75.       melval=ipt1
  76.       lmvid=.false.
  77. *     write (6,*) ' faced2 velche '
  78. *     write (6,*)  (velche(ix,1),ix=1,20)
  79.       if (mcham.eq.0) then
  80.         CALL ECROBJ('MAILLAGE',IPT1)
  81.         CALL ENVELO
  82.         CALL LIROBJ('MAILLAGE',MELEME,1,IRETOU)
  83.       else
  84.         CALL ENVEL1(ipt1,meleme,mcoup)
  85.       endif
  86.       IF (IERR.NE.0) RETURN
  87.  
  88.       SEGACT MELEME
  89.       SEGACT XPROJ,ICPR
  90.       NBPOIN=XPROJ(/2)
  91.       TMIN=1E30
  92.       TMAX=-1E30
  93.       DO 1 I=1,NBPOIN
  94.          TMIN=MIN(TMIN,XPROJ(3,I))
  95.          TMAX=MAX(TMAX,XPROJ(3,I))
  96.    1  CONTINUE
  97.       TDIST=TMAX-TMIN
  98.       NPOI1=0
  99.       NSEG2=0
  100.       NSEG3=0
  101. c
  102. c     on compte le nombre d elements dont les faces sont de type 1 2 3 4
  103. c     6 7 8 dans NBFAC, attention à 6 : gestion des polygones
  104.       SEGINI NBFAC
  105.       IPT1=MELEME
  106.       SEGACT MELEME
  107.       DO 10 IOB=1,MAX(1,LISOUS(/1))
  108.          IF (LISOUS(/1).NE.0) THEN
  109.             IPT1=LISOUS(IOB)
  110.             SEGACT IPT1
  111.          ENDIF
  112.          NBELEM=IPT1.NUM(/2)
  113.          ILTEL=LTEL(1,IPT1.ITYPEL)
  114.          IF (ILTEL.EQ.0) GOTO 12
  115.          ILTAD=LTEL(2,IPT1.ITYPEL)
  116.          DO 13 IF=1,ILTEL
  117.             IFT=LDEL(1,ILTAD+IF-1)
  118.             NBFAC(IFT)=NBFAC(IFT)+NBELEM
  119.  13      CONTINUE
  120.  12      CONTINUE
  121.          if (ipt1.itypel.eq.1) npoi1=npoi1+nbelem
  122.          if (ipt1.itypel.eq.2) nseg2=nseg2+nbelem
  123.          if (ipt1.itypel.eq.3) nseg3=nseg3+nbelem
  124.  10   CONTINUE
  125. *
  126. *      WRITE(IOIMP,*) 'NBFAC'
  127. *      WRITE (IOIMP,9111) (NBFAC(III),III=1,NTYFAC)
  128. * 9111 FORMAT(5(2X,I6))
  129.       NFACP=NBFAC(6)
  130.       SEGINI NSOMP
  131.       SEGINI ipoi1,iseg2,iseg3
  132.       if (mcham.ne.0) SEGINI xpoi1,xseg2,xseg3
  133. c==== CREATION DES FACES ==============================================
  134. * Initialisation des IFACI,XFACI
  135.       SEGINI IPOFAC
  136.       DO ITYFAC=1,NTYFAC
  137.          NNODE=KDFAC(1,ITYFAC)
  138. * Polygone
  139.          IF (ITYFAC.EQ.6) NNODE=NNOMAX
  140.          IF (NNODE.GT.0) THEN
  141.             NFACI=NBFAC(ITYFAC)
  142.             SEGINI IFACI
  143.             IPOFAC(1,ITYFAC)=IFACI
  144.             if (mcham.ne.0) then 
  145.                segini xfaci
  146.                IPOFAC(2,ITYFAC)=xfaci
  147.             endif
  148.          ENDIF
  149.       ENDDO
  150.       NPOI1=0
  151.       NSEG2=0
  152.       NSEG3=0
  153. c     NBFAC sert maintenant de compteur
  154.       DO ITYFAC=1,NTYFAC
  155.          NBFAC(ITYFAC)=0
  156.       ENDDO
  157.       ICOUPE=0
  158.       melval = meleme
  159.       npasnul = 0
  160.       DO 50 IOB=1,MAX(1,LISOUS(/1))
  161.          IF (LISOUS(/1).NE.0) THEN
  162.             IPT1=LISOUS(IOB)
  163.             ivm1=0
  164.             ivm2=0
  165.             melval=meleme
  166.             if (mcham.ne.0) then
  167.                melval=lisref(iob)
  168.                if (melval.eq.0) then
  169.                   write (6,*) 'reference nulle dans faced2',iob
  170.                   goto 50
  171.                else
  172.                  npasnul = npasnul + 1
  173.                endif
  174.                ivm1=velche(/1)
  175.                ivm2=velche(/2)
  176.                lmvid=(ivm1*ivm2.eq.0)
  177. *sg 2016/08/23 ancienne programmation brutale qui sort de la subroutine
  178. *     sans mettre meleme=melsau => segmentation violation plus
  179. *     tard
  180. *
  181. *               if (ivm1*ivm2.eq.0) return
  182.  
  183.             endif
  184.             ICOUPE=0
  185.             IF (IOB.EQ.LISOUS(/1).AND.MCOUP.NE.0) ICOUPE=1
  186.          ENDIF
  187.          NBELEM=IPT1.NUM(/2)
  188.          NBNN  =IPT1.NUM(/1)
  189.          ILTEL=LTEL(1,IPT1.ITYPEL)
  190.          IF (ILTEL.EQ.0) GOTO 52
  191.          ILTAD=LTEL(2,IPT1.ITYPEL)
  192.          DO 60 IF=1,ILTEL
  193.             ITYFAC=LDEL(1,ILTAD+IF-1)
  194.             IAD=LDEL(2,ILTAD+IF-1)
  195.             NNODE=KDFAC(1,ITYFAC)
  196.             NNODF=NNODE
  197. * Polygone
  198.             IF (ITYFAC.EQ.6) THEN
  199.                NBNN=IPT1.NUM(/1)
  200. *  23 1
  201. *Erreur dans le module de trace
  202.                IF (NBNN.GT.NNOMAX) THEN
  203.                   CALL ERREUR(23)
  204.                   RETURN
  205.                ENDIF
  206.                NNODE=NBNN
  207.                NNODF=NNOMAX
  208.             ENDIF
  209.             IF (NNODE.GT.0) THEN 
  210. *               WRITE(IOIMP,*) 'ITYFAC=',ITYFAC
  211.                IFACI=IPOFAC(1,ITYFAC)
  212.                xfaci=IPOFAC(2,ITYFAC)
  213. *  optimiseur
  214.                if(xfaci.eq.0) xfaci=ipt1
  215.                DO 80 IEL=1,NBELEM
  216. *                  WRITE(IOIMP,*) 'IEL=',IEL
  217.                   do 70 inn=1,nbnn
  218.                      if (icpr(ipt1.num(inn,iel)).eq.0) then
  219.                         call erreur(23)
  220.                         goto  80
  221.                      endif
  222.  70               continue
  223.                   NBFAC(ITYFAC)=NBFAC(ITYFAC)+1
  224.                   j=NBFAC(ITYFAC)
  225. * Polygone
  226.                   IF (ITYFAC.EQ.6) NSOMP(j)=NNODE
  227.                   IFACI(NNODF+1,j)=1
  228.                   DO i=1,NNODE
  229.                      IFACI(i,j)=IPT1.NUM(LFAC(IAD+i-1),IEL)
  230.                      if (mcham.ne.0.AND.(.not.lmvid)) then
  231.                         xfaci(i,j)=velche(min(ivm1,LFAC(IAD+i-1)),
  232.      >                       min(ivm2,IEL))
  233. *                        WRITE(IOIMP,*) 'INODE,IFACI,XFACI=',i,ifaci(i,j)
  234. *     $                       ,xfaci(i,j)
  235.                      endif
  236.                      IF (IVU(ICPR(IFACI(i,j))).NE.1) IFACI(NNODF+1,j)=0
  237.                   ENDDO
  238. * TRI3 cas des coupes
  239.                   IF (ITYFAC.EQ.1) THEN
  240.                      IF (ICOUPE.EQ.1) THEN
  241.                         IF (MCOUP(IEL)/8.EQ.1) IFACI(NNODF+1,j)=2
  242.                         IF (MCOUP(IEL)/16.EQ.1) IFACI(NNODF+1,j)=3
  243.                      ENDIF
  244.                   ENDIF
  245.  80            CONTINUE
  246.             ENDIF
  247.  60      CONTINUE
  248.          IF (LISOUS(/1).NE.0) SEGDES IPT1
  249.          goto 50
  250.  52      CONTINUE
  251.          do 68 iel=1,nbelem
  252.             if (ipt1.itypel.eq.1) then
  253.                npoi1=npoi1+1
  254.                IPOI1(1,NPOI1)=IPT1.NUM(1,IEL)
  255.             elseif (ipt1.itypel.eq.2) then
  256.                nseg2=nseg2+1
  257.                ISEG2(1,NSEG2)=IPT1.NUM(1,IEL)
  258.                ISEG2(2,NSEG2)=IPT1.NUM(2,IEL)
  259.             elseif (ipt1.itypel.eq.3) then
  260.                nseg3=nseg3+1
  261.                ISEG3(1,NSEG3)=IPT1.NUM(1,IEL)
  262.                ISEG3(2,NSEG3)=IPT1.NUM(2,IEL)
  263.                ISEG3(3,NSEG3)=IPT1.NUM(3,IEL)
  264.             endif
  265.             if (mcham.ne.0.AND.(.not.lmvid)) then
  266.                if (ipt1.itypel.eq.1) then
  267.                   xPOI1(1,NPOI1)=velche(min(ivm1,1),min(ivm2,IEL))
  268.                elseif (ipt1.itypel.eq.2) then
  269.                   xSEG2(1,NSEG2)=velche(min(ivm1,1),min(ivm2,IEL))
  270.                   xSEG2(2,NSEG2)=velche(min(ivm1,2),min(ivm2,IEL))
  271.                elseif (ipt1.itypel.eq.3) then
  272.                   xSEG3(1,NSEG3)=velche(min(ivm1,1),min(ivm2,IEL))
  273.                   xSEG3(2,NSEG3)=velche(min(ivm1,2),min(ivm2,IEL))
  274.                   xSEG3(3,NSEG3)=velche(min(ivm1,3),min(ivm2,IEL))
  275.                endif
  276.             endif
  277.  68      continue
  278.          IF (LISOUS(/1).NE.0) SEGDES IPT1
  279.  50   CONTINUE
  280.       SEGDES MELEME
  281.  
  282. C SI que des pointeurs nuls dans LISREF : ERREUR
  283.       IF (mcham.ne.0.AND.npasnul.EQ.0) THEN
  284.         CALL ERREUR(23)
  285.         RETURN
  286.       ENDIF
  287. C  IF FAUT MAINTENANT RETASSER ET CLASSER LES TABLEAUX DES FACES
  288. C  PROBLEME ELLES NE SONT PAS FORCEMENT DECRITES DE LA MEME FACON
  289.       NFAC=0
  290.       DO ITYFAC=1,NTYFAC
  291.          NFAC=NFAC+NBFAC(ITYFAC)
  292.       ENDDO
  293.       NFAC=NFAC+NPOI1+NSEG2+NSEG3
  294.       SEGINI TFAC,KFAK
  295.       IFAC=0
  296.       ITYEL=0
  297. * D'abord les faces ensuite les segments
  298.       DO ITYFAC=1,NTYFAC
  299.          NNODE=KDFAC(1,ITYFAC)
  300.          IF (ITYFAC.EQ.6) THEN
  301.             NNODF=NNOMAX
  302.          ELSE
  303.             NNODF=NNODE
  304.          ENDIF
  305.          IF (NNODE.GT.0.OR.ITYFAC.EQ.6) THEN
  306.             IFACI=IPOFAC(1,ITYFAC)
  307.             DO I=1,NBFAC(ITYFAC)
  308.                IFAC=IFAC+1
  309. * Polygone
  310.                IF (ITYFAC.EQ.6) NNODE=NSOMP(I)
  311.                XXXX = 0.
  312.                DO J=1,NNODE
  313.                   XXXX = XXXX + XPROJ(3,ICPR(IFACI(J,I)))
  314.                ENDDO
  315.                XXXX=XXXX/NNODE
  316.                TFAC(IFAC)=XXXX
  317.                IF (IFACI(NNODF+1,I).EQ.1) TFAC(IFAC)=TFAC(IFAC)-TDIST
  318.                KFAK(IFAC)=I+(ITYEL*NFAC)
  319. * TRI3/coupe
  320.                IF (ITYFAC.EQ.1) THEN
  321.                   IF (IFACI(NNODF+1,I).EQ.2) 
  322.      $                 TFAC(IFAC)=TFAC(IFAC)-2*TDIST
  323.                   IF (IFACI(NNODF+1,I).EQ.3) KFAK(IFAC)=0
  324.                ENDIF
  325.             ENDDO
  326.          ENDIF
  327.          ITYEL=ITYEL+1
  328.       ENDDO
  329.       DO I=1,NPOI1
  330.          IFAC=IFAC+1
  331.          TFAC(IFAC)=XPROJ(3,ICPR(IPOI1(1,I)))
  332.          KFAK(IFAC)=I+(ITYEL*NFAC)
  333.       ENDDO      
  334.       ITYEL=ITYEL+1
  335. * SEG2 si NSEG2=0, la boucle est sautee
  336.       DO I=1,NSEG2
  337.          IFAC=IFAC+1
  338.          TFAC(IFAC)=(XPROJ(3,ICPR(ISEG2(1,I)))+
  339.      $               XPROJ(3,ICPR(ISEG2(2,I))))/2.
  340.          KFAK(IFAC)=I+(ITYEL*NFAC)
  341. *dbg         write(ioimp,*) 'ifac,tfac,kfak=',ifac,tfac(ifac),kfak(ifac)
  342.       ENDDO
  343.       ITYEL=ITYEL+1
  344. * SEG3 dans l'ancien faced2, il y avait /2. ????
  345.       DO I=1,NSEG3
  346.          IFAC=IFAC+1
  347.          TFAC(IFAC)=(XPROJ(3,ICPR(ISEG3(1,I)))+
  348.      #               XPROJ(3,ICPR(ISEG3(2,I)))+
  349.      $               XPROJ(3,ICPR(ISEG3(3,I))))/3.
  350.          KFAK(IFAC)=I+(ITYEL*NFAC)
  351.       ENDDO
  352. C  IL N'Y A PLUS QU'A TRIER ET RETASSER KFAK SUIVANT TFAC
  353. *dbg      WRITE (IOIMP,9111) (KFAK(I),TFAC(I),I=1,NFAC)
  354.  9111 FORMAT(5(2X,I6,1X,E12.5))
  355.  
  356.       SEGINI NAUX
  357.       IF (IDIM.EQ.2) GOTO 209
  358.       NAUX(1)=1
  359.       NAUX(2)=NFAC
  360.       IZ=2
  361.  208  IZ=IZ-1
  362.       IF (IZ.LE.0) GOTO 209
  363.       IPB=NAUX(IZ*2-1)
  364.       IPH=NAUX(IZ*2)
  365.       IF(IPB.GE.IPH) GOTO 208
  366.       JPB=IPB-1
  367.       JPH=IPH+1
  368. C  CALCUL DU PIVOT
  369.       PV=0.
  370. *     DO 207 J=IPB,IPH
  371. *     PV=PV+TFAC(J)
  372. *207  CONTINUE
  373. *     PV=PV/(IPH-IPB+1)
  374.       PV=(TFAC(IPB)+TFAC(IPH))/2.
  375.  242  JPB=JPB+1
  376.       IF (JPH.EQ.JPB) GOTO 245
  377.       IF (TFAC(JPB).LT.PV) GOTO 243
  378.       GOTO 242
  379.  243  JPH=JPH-1
  380.       IF (JPH.EQ.JPB) GOTO 245
  381.       IF (TFAC(JPH).GT.PV) GOTO 244
  382.       GOTO 243
  383.  244  IAUX=KFAK(JPB)
  384.       KFAK(JPB)=KFAK(JPH)
  385.       KFAK(JPH)=IAUX
  386.       TAUX=TFAC(JPB)
  387.       TFAC(JPB)=TFAC(JPH)
  388.       TFAC(JPH)=TAUX
  389.       GOTO 242
  390.  245  IF (JPB.GE.IPB) GOTO 247
  391.       JPB=JPB+1
  392.       JPH=JPH+2
  393.       GOTO 248
  394.  247  IF (JPH.LE.IPH) GOTO 249
  395.       JPB=JPB-2
  396.       JPH=JPH-1
  397.       GOTO 248
  398.  249  IF (TFAC(JPB).LE.PV) GOTO 250
  399.       IF (JPH.EQ.IPH) GOTO 251
  400.  252  JPH=JPH+1
  401.       GOTO 248
  402.  250  IF (JPB.EQ.IPB) GOTO 252
  403.  251  JPB=JPB-1
  404.  248  IF (JPB.EQ.IPB) GOTO 253
  405.       NAUX(2*IZ)=JPB
  406.       IZ=IZ+1
  407.  253  IF (JPH.EQ.IPH) GOTO 208
  408.       NAUX(2*IZ)=IPH
  409.       NAUX(2*IZ-1)=JPH
  410.       IZ=IZ+1
  411.       GOTO 208
  412.  209  CONTINUE
  413. C  LES FACES SONT CLASSEES DANS KFAK LES FACES ON ETE ELIMINEES PAR
  414. C  ENVELO .   IL NE RESTE PLUS QU'A TRACER
  415. *dbg      WRITE (IOIMP,9111) (KFAK(I),TFAC(I),I=1,NFAC)
  416.       DO 300 I=1,NFAC
  417.          IFF=KFAK(I)
  418.          IF (IFF.EQ.0) GOTO 300
  419.          IT=(IFF-1)/NFAC
  420.          IF=IFF-IT*NFAC
  421.          IT=IT+1
  422.          IOK=0
  423. *     
  424.          ITYFAC=IT
  425. *dbg         WRITE(IOIMP,*) 'IFAC,ITYFAC=',I,ITYFAC
  426.          IF (ITYFAC.GE.1.AND.ITYFAC.LE.NTYFAC) THEN
  427.             NNODE=KDFAC(1,ITYFAC)
  428.             NNODF=NNODE
  429. * Polygone
  430.             IF (ITYFAC.EQ.6) THEN
  431.                NNODE=NSOMP(IF)
  432.                NNODF=NNOMAX
  433.             ENDIF
  434. * A cette etape on doit avoir nnode.gt.0
  435.             IF (NNODE.LE.0) THEN
  436.                CALL ERREUR(5)
  437.                RETURN
  438.             ENDIF
  439.             IFACI=IPOFAC(1,ITYFAC)
  440.             xfaci=IPOFAC(2,ITYFAC)
  441.             DO J=1,NNODE
  442.                IF (IVU(ICPR(IFACI(J,IF))).EQ.1) IOK=1
  443.             ENDDO
  444.             IF (IOK.EQ.1) THEN
  445. * Cas du TRI3
  446.                IF (ITYFAC.EQ.1) THEN
  447.                   DO J=1,NNODE
  448.                      NUPT=IFACI(J, IF)
  449.                      IDPT=ICPR(NUPT)
  450.                      XX(J)=XPROJ(1,IDPT)
  451.                      YY(J)=XPROJ(2,IDPT)
  452.                      ZZ(J)=XPROJ(3,IDPT)
  453.                      if (mcham.eq.0) then
  454.                         VV(J)=VCPCHA(NUPT)
  455.                      else
  456.                         VV(J)=xfaci(j,IF)
  457.                      endif
  458.                   ENDDO
  459.                   CALL TRISO(VCHC,XX,YY,ZZ,VV,NNODE,NISO)
  460. *     Cas des autres faces. Elles ont pour particularité d'avoir une
  461. *     valeur centrale et un contour, on trace avec TCISO par defaut
  462. *     On peut aussi faire une boucle de TRISO mais ca donne un
  463. *     Postscript plus gros.
  464.                ELSE
  465.                   ICONT=0
  466.                   ITCISO=1
  467. *     La valeur centrale est la moyenne sur les autres noeuds QUA4, POLY
  468.                   IF (ITYFAC.EQ.2.OR.ITYFAC.EQ.6) THEN
  469.                      XXM = 0.
  470.                      YYM = 0.
  471.                      ZZM = 0.
  472.                      VVM = 0.
  473.                      DO J=1,NNODE
  474.                         NUPT=IFACI(J, IF)
  475.                         IDPT=ICPR(NUPT)
  476.                         XXM = XPROJ(1,IDPT) + XXM
  477.                         YYM = XPROJ(2,IDPT) + YYM
  478.                         ZZM = XPROJ(3,IDPT) + ZZM
  479.                         IF (MCHAM.EQ.0) THEN
  480.                            VVM = VCPCHA(NUPT) + VVM
  481.                         ELSE
  482.                            VVM = xfaci(J, IF) + VVM
  483.                         ENDIF
  484.                      ENDDO
  485.                      XXM=XXM/NNODE
  486.                      YYM=YYM/NNODE
  487.                      ZZM=ZZM/NNODE
  488.                      VVM=VVM/NNODE
  489. *     La valeur centrale est une moyenne pondérée
  490.                   ELSEIF (ITYFAC.EQ.3.OR.ITYFAC.EQ.4) THEN
  491.                      XXM = 0.
  492.                      YYM = 0.
  493.                      ZZM = 0.
  494.                      VVM = 0.
  495.                      IDEC=0
  496.                      IF (ITYFAC.EQ.4) IDEC=6
  497.                      DO J=1,NNODE
  498.                         NUPT=IFACI(J, IF)
  499.                         IDPT=ICPR(NUPT)
  500.                         XXM = XPROJ(1,IDPT)*XPOIDS(J+IDEC) + XXM
  501.                         YYM = XPROJ(2,IDPT)*XPOIDS(J+IDEC) + YYM
  502.                         ZZM = XPROJ(3,IDPT)*XPOIDS(J+IDEC) + ZZM
  503.                         IF (MCHAM.EQ.0) THEN
  504.                            VVM = VCPCHA(NUPT)*XPOIDS(J+IDEC) + VVM
  505.                         ELSE
  506.                            VVM = xfaci(J, IF)*XPOIDS(J+IDEC) + VVM
  507.                         ENDIF
  508.                      ENDDO
  509. *     La valeur centrale est celle du dernier noeud (faces TRI7/QUA9)
  510.                   ELSEIF (ITYFAC.EQ.7.OR.ITYFAC.EQ.8) THEN
  511. *                     write(ioimp,*) 'coucou ityfac=',ityfac
  512.                      NUPT=IFACI(NNODE, IF)
  513.                      IDPT=ICPR(NUPT)
  514.                      XXM = XPROJ(1,IDPT)
  515.                      YYM = XPROJ(2,IDPT)
  516.                      ZZM = XPROJ(3,IDPT)
  517.                      IF (MCHAM.EQ.0) THEN
  518.                         VVM = VCPCHA(NUPT)
  519. *                        write(ioimp,*) 'vvm=',vvm
  520.                      ELSE
  521.                         VVM = xfaci(NNODE, IF)
  522. *                        write(ioimp,*) 'vvm=',vvm
  523.                      ENDIF
  524. * On met ICONT à 1 pour ne pas mettre le dernier noeud dans le contour
  525.                      ICONT=1
  526. *!!!
  527.                      ITCISO=1
  528.                   ELSE
  529.                      write(ioimp,*) 'ITYFAC=',ityfac,' non prevu'
  530.                      call erreur(5)
  531.                      return
  532.                   ENDIF
  533.                   XX(1)=XXM
  534.                   YY(1)=YYM
  535.                   ZZ(1)=ZZM
  536.                   VV(1)=VVM
  537.                   DO J=1,NNODE-ICONT
  538.                      JP=J+1
  539.                      NUPT=IFACI(J, IF)
  540.                      IDPT=ICPR(NUPT)
  541.                      XX(JP) = XPROJ(1,IDPT)
  542.                      YY(JP) = XPROJ(2,IDPT)
  543.                      ZZ(JP) = XPROJ(3,IDPT)
  544.                      IF (MCHAM.EQ.0) THEN
  545.                         VV(JP) = VCPCHA(NUPT)
  546.                      ELSE
  547.                         VV(JP) = xfaci(J, IF)
  548.                      ENDIF
  549.                   ENDDO
  550.                   IF (ITCISO.EQ.1) THEN
  551.                      CALL TCISO(VCHC,XX,YY,ZZ,VV,NNODE-ICONT+1,NISO)
  552.                   ELSE
  553.                      XR(1)=XX(1)
  554.                      YR(1)=YY(1)
  555.                      ZR(1)=ZZ(1)
  556.                      VR(1)=VV(1)
  557.                      JMAX=NNODE-ICONT
  558.                      DO J=1,JMAX
  559.                         JP=J+1
  560.                         IF (JP.GT.JMAX) JP=1
  561.                         IA=J+1
  562.                         IB=JP+1
  563.                         XR(2)=XX(IA)
  564.                         YR(2)=YY(IA)
  565.                         ZR(2)=ZZ(IA)
  566.                         VR(2)=VV(IA)
  567.                         XR(3)=XX(IB)
  568.                         YR(3)=YY(IB)
  569.                         ZR(3)=ZZ(IB)
  570.                         VR(3)=VV(IB)
  571.                         CALL TRISO(VCHC,XR,YR,ZR,VR,3,NISO)
  572.                      ENDDO
  573.                   ENDIF
  574.                ENDIF
  575.             ENDIF
  576. * Cas des POI1
  577.          ELSEIF (ITYFAC.EQ.NTYFAC+1) THEN
  578.             NNODE=2
  579.             IF (IVU(ICPR(IPOI1(1,IF))).EQ.1) IOK=1
  580.             IF (IOK.EQ.1) THEN
  581.                NUPT=IPOI1(1,IF)
  582.                IDPT=ICPR(NUPT)
  583.  
  584.                if (mcham.eq.0) then
  585.                   VV(1)=VCPCHA(NUPT)
  586.                   VV(2)=VV(1)
  587.                else
  588.                   VV(1)=xpoi1(1,IF)
  589.                   VV(2)=VV(1)
  590.                endif
  591.  
  592.                XX(1)=XPROJ(1,IDPT)+BLOK
  593.                YY(1)=XPROJ(2,IDPT)
  594.                ZZ(1)=XPROJ(3,IDPT)
  595.                XX(2)=XPROJ(1,IDPT)-BLOK
  596.                YY(2)=XPROJ(2,IDPT)
  597.                ZZ(2)=XPROJ(3,IDPT)
  598.                CALL TRISO(VCHC,XX,YY,ZZ,VV,NNODE,NISO)
  599.  
  600.                XX(1)=XPROJ(1,IDPT)
  601.                YY(1)=XPROJ(2,IDPT)+BLOK
  602.                ZZ(1)=XPROJ(3,IDPT)
  603.                XX(2)=XPROJ(1,IDPT)
  604.                YY(2)=XPROJ(2,IDPT)-BLOK
  605.                ZZ(2)=XPROJ(3,IDPT)
  606.                CALL TRISO(VCHC,XX,YY,ZZ,VV,NNODE,NISO)
  607.             ENDIF
  608. * Cas des SEG2
  609.          ELSEIF (ITYFAC.EQ.NTYFAC+2) THEN
  610.             NNODE=2
  611.             DO J=1,NNODE
  612.                IF (IVU(ICPR(ISEG2(J,IF))).EQ.1) IOK=1
  613.             ENDDO
  614.             IF (IOK.EQ.1) THEN
  615.                DO J=1,NNODE
  616.                   NUPT=ISEG2(J, IF)
  617.                   IDPT=ICPR(NUPT)
  618.                   XX(J)=XPROJ(1,IDPT)
  619.                   YY(J)=XPROJ(2,IDPT)
  620.                   ZZ(J)=XPROJ(3,IDPT)
  621.                   if (mcham.eq.0) then
  622.                      VV(J)=VCPCHA(NUPT)
  623.                   else
  624.                      VV(J)=xseg2(j,IF)
  625.                   endif
  626.                ENDDO
  627.                CALL TRISO(VCHC,XX,YY,ZZ,VV,NNODE,NISO)
  628.             ENDIF
  629. * Cas des SEG3
  630.          ELSEIF (ITYFAC.EQ.NTYFAC+3) THEN
  631. * 2 SEG2 !
  632.             DO ISEG=1,2
  633.                isegm=iseg-1
  634.                IOK=0
  635.                DO J=1,2
  636.                   IF (IVU(ICPR(ISEG3(J+isegm,IF))).EQ.1) IOK=1
  637.                ENDDO
  638.                IF (IOK.EQ.1) THEN
  639.                   DO J=1,2
  640.                      NUPT=ISEG3(J+isegm, IF)
  641.                      IDPT=ICPR(NUPT)
  642.                      XX(J)=XPROJ(1,IDPT)
  643.                      YY(J)=XPROJ(2,IDPT)
  644.                      ZZ(J)=XPROJ(3,IDPT)
  645.                      if (mcham.eq.0) then
  646.                         VV(J)=VCPCHA(NUPT)
  647.                      else
  648.                         VV(J)=xseg3(j+isegm,IF)
  649.                      endif
  650.                   ENDDO
  651.                   CALL TRISO(VCHC,XX,YY,ZZ,VV,2,NISO)
  652.                ENDIF
  653.             ENDDO
  654.          ELSE
  655.             write(ioimp,*) 'ITYFAC=',ITYFAC,' non prevu'
  656.             call erreur(5)
  657.             return
  658.          ENDIF
  659.  300  CONTINUE
  660. C 'EST FINI
  661.       SEGACT MELEME
  662.       IF (LISOUS(/1).NE.0) THEN
  663.          NBSOUS=LISOUS(/1)
  664.          IF (MCOUP.NE.0) NBSOUS=NBSOUS-1
  665.          DO 490 IO=1,NBSOUS
  666.             IPT2=LISOUS(IO)
  667.             segact ipt2
  668.             if (ipt2.itypel.gt.3.AND.ipt2.itypel.NE.32) then
  669.                SEGSUP IPT2
  670.                if (mcham.ne.0) then
  671.                   melval=lisref(io)
  672.                   if (melval.ne.0) segsup melval
  673.                endif
  674.             else
  675.                segdes ipt2
  676.             endif
  677.  490     CONTINUE
  678.       ENDIF
  679.       if (itypel.eq.0) SEGSUP MELEME
  680.       MELEME=MELSAU
  681.       SEGSUP TFAC,KFAK
  682.       SEGSUP NAUX
  683.       SEGSUP IPOI1,ISEG2,ISEG3
  684.       if (mcham.ne.0) SEGSUP XPOI1,XSEG2,XSEG3
  685.       SEGSUP NSOMP
  686.       DO ITYFAC=1,NTYFAC
  687.          IFACI=IPOFAC(1,ITYFAC)
  688.          IF (IFACI.NE.0) THEN
  689.             SEGSUP IFACI
  690.          ENDIF
  691.          xfaci=IPOFAC(2,ITYFAC)
  692.          IF (XFACI.NE.0) THEN
  693.             SEGSUP XFACI
  694.          ENDIF
  695.       ENDDO
  696.       SEGSUP IPOFAC
  697.       SEGSUP NBFAC
  698.       RETURN
  699.       END
  700.  
  701.  
  702.  
  703.  
  704.  
  705.  
  706.  
  707.  
  708.  
  709.  

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