Télécharger la plaquette Cast3M


fr en

Télécharger demait.eso

Retour à la liste

Numérotation des lignes : 

demait

  1. C DEMAIT    SOURCE    JC220346  16/11/29    21:15:09     9221           
  2. C|-------------------------------------------------------------------|
  3. C|                                                                   |
  4. C|                       PROGRAMME PRINCIPAL                         |
  5. C|                        MAIN  DE  DEMETER                          |
  6. C|                                                                   |
  7. C|-------------------------------------------------------------------|
  8. C
  9.        SUBROUTINE DEMAIT(IDCP,NPTBAS)
  10. C
  11. C
  12.       IMPLICIT INTEGER(I-N)
  13.       IMPLICIT REAL*8(a-h,o-z)
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17. -INC TDEMAIT
  18.       SEGMENT IDCP(NPTCOM)
  19.        REAL*8 XDEM
  20.        LOGICAL REPONS,FACET,INTER,PINT,FERME,VAL,IN,DROIT,
  22.        real*8 xval,xo(3),xa(3),xb(3),xc(3),epsi,epsj,xmesu,ymesu
  23.        character*4 mcle(8)
  24.        data mcle /'TCRI','CFAC','CDIS','TETR','EXPC','FINA','DIAC',
  25.      > 'EXPF'/
  26.       ymesu=0
  27.       NPTORI=NPTMAX
  28.       NFCORI=NFCMAX
  29. *  flags type d'operation
  30.       ipass=1
  31.       tcrit=3.
  32.         cfacei=6.0
  33.         cfacet=cfacei
  34.         cdist=0.125
  35.         tetrl=2.75
  36.         expcom=1.00
  37.         diacri=0.93
  38.         diacrd=diacri
  39.         diacre=diacri
  40.         expfac=sqrt(3.)/2.
  41.         faccri=16
  42.         volcri=0.01
  43.   10    continue
  44.         call lirmot(mcle,8,imot,0)
  45.          if (imot.eq.1) then
  46.           call lirree(xval,1,iret)
  47.           tcrit=xval
  48.          endif
  49.          if (imot.eq.2) then
  50.           call lirree(xval,1,iret)
  51.           cfacei=xval
  52.          endif
  53.          if (imot.eq.3) then
  54.           call lirree(xval,1,iret)
  55.           cdist=xval
  56.          endif
  57.          if (imot.eq.4) then
  58.           call lirree(xval,1,iret)
  59.           tetrl=xval
  60.          endif
  61.          if (imot.eq.5) then
  62.           call lirree(xval,1,iret)
  63.           expcom=xval
  64.          endif
  65.          if (imot.eq.7) then
  66.           call lirree(xval,1,iret)
  67.           diacre=xval
  68.          endif
  69.          if (imot.eq.8) then
  70.           call lirree(xval,1,iret)
  71.           expfac=xval
  72.          endif
  73.          if (imot.ne.0) goto 10
  74. C
  75. C      INITIALISATION DU TABLEAU DES VOLUMES
  76. C      -------------------------------------
  77.       NFTOT=IFUT(/1)
  78.        NVTOT=IVOL(/2)
  79.        DO 130 J=1,NVTOT
  80.           DO 120 I=1,9
  81.              IVOL(I,J)=0
  82. 120       CONTINUE
  83. 130    CONTINUE
  84. C
  85. C
  86. C      CONSTRUCTION DU TABLEAU NPF ( POINTS-FACETTES )
  87. C      -----------------------------------------------
  88.        DO 150 J=1,40
  89.           DO 150 I=1,NPTMAX
  90.               NPF(J,I)=0
  91. 150    CONTINUE
  92.           DO 141 J=1,NFCMAX
  93.              DO 140 K=1,4
  94.                IP=NFC(K,J)
  95.                 IF (IP.EQ.0) GOTO 140
  96.                 L=-NPF(40,IP)+1
  97.                   IF (L.LE.0) CALL ERREUR(126)
  98.                   IF (IERR.NE.0) RETURN
  99.                  NPF(40,IP)=-L
  100.                 NPF(L,IP)=J
  101. 140          CONTINUE
  102. 141          CONTINUE
  103.         DO 145 I=1,NPTMAX
  104.           NPF(40,I)=MAX(0,NPF(40,I))
  105. 145     CONTINUE
  106. C
  107. C      RECHERCHE DE LA TAILLE MOYENNE DE MAILLE ASSOCIEE A
  108. C      CHAQUE POINT ( 4EME COMPOSANTE DU POINT )
  109.        DO 190 I=1,NPTMAX
  110.           DD=0.
  111.           KK=0
  112.           DO 170 J=1,40
  113.              IF (NPF(J,I).EQ.0) GOTO 180
  114.              if=npf(j,i)
  115.              nc=4
  116.              if (nfc(4,if).eq.0) nc=3
  117.              jp=nfc(nc,if)
  118.              do 175 ic=1,nc
  119.              ip=nfc(ic,if)
  120.              XX=(XYZ(1,IP)-XYZ(1,JP))**2
  121.              YY=(XYZ(2,IP)-XYZ(2,JP))**2
  122.              ZZ=(XYZ(3,IP)-XYZ(3,JP))**2
  123.              DD=DD+SQRT(XX+YY+ZZ)
  124.              KK=KK+1
  125.              jp=ip
  126.  175         continue
  127. *            IP=ISUCC(NPF(J,I),I)
  128. *            XX=(XYZ(1,I)-XYZ(1,IP))**2
  129. *            YY=(XYZ(2,I)-XYZ(2,IP))**2
  130. *            ZZ=(XYZ(3,I)-XYZ(3,IP))**2
  131. *            DD=DD+SQRT(XX+YY+ZZ)
  132. *            KK=KK+1
  133. 170       CONTINUE
  134. 180       XYZ(4,I)=DD/KK
  135. 190    CONTINUE
  136. *  regularisation locale
  137. *     DO 182 I=1,NPTMAX
  138. *         DD= XYZ(4,I)
  139. *         KK=1
  140. *         DO 184 J=1,40
  141. *            IF (NPF(J,I).EQ.0) GOTO 186
  142. *            IP=ISUCC(NPF(J,I),I)
  143. *            DD=DD+XYZ(4,IP)
  144. *            KK=KK+1
  145. *84       CONTINUE
  146. *86       XYZ(4,I)=DD/KK
  147. *82   CONTINUE
  148. *  taille moyenne generale
  149.       xmoy=0
  150.       DO 181 I=1,NPTMAX
  151.       XMOY=XMOY+LOG(XYZ(4,I))
  152.  181  CONTINUE
  153.       XMOYG=EXP(XMOY/NPTMAX)
  154.       IF (IVERB.EQ.1) WRITE (6,*) ' TAILLE MOYENNE VISEE ',XMOYG
  155. C
  156. C      LE MAILLAGE DE LA SURFACE EST-IL FERME ?
  157. C      ----------------------------------------
  158.        REPONS=FERME(KKK)
  159.        IF (REPONS) GOTO 210
  160.        CALL ERREUR(127)
  161.        RETURN
  162.  210  continue
  163.       xmesu=0
  164.       do 100  if=1,nfcmax
  165.        xmesu=xmesu+vol(1,nfc(1,if),nfc(2,if),nfc(3,if))
  166.        if (nfc(4,if).ne.0) then
  167.          xmesu=xmesu+vol(1,nfc(1,if),nfc(3,if),nfc(4,if))
  168.          ipass=2
  169.        endif
  170.  100  continue
  171.        xmesu=xmesu/(-6.)
  172.       IF (IVERB.EQ.1) WRITE (6,*) ' volume de la piece   ',xmesu
  173. C
  174. C      NFACET : NOMBRE DE FACETTES DU MAILLAGE DE SURFACE
  175. C      --------------------------------------------------
  176.        NFACET=NFCMAX
  177.        NPTCOM=NPTMAX
  178.        NVOL=0
  179.        NFCPRE=0
  180.        NFCTRT=0
  181.        NARET=0
  182.        NPTOT=XYZ(/2)
  183.        NTTRAV=NFCMAX
  184.        SEGINI TRAV
  185.        YMESU=0
  186.        NVOLY=0
  187.        NPTDEB=NPTMAX
  188.        NPTDIS=1
  189.        ICRTS=0
  190. 220    CONTINUE
  191.        do 222 jvol=nvoly+1,nvol
  192.         if (ivol(9,jvol).eq.25) then
  193.           ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol),
  194.      >                    ivol(3,jvol),ivol(4,jvol))/6.
  195.         endif
  196.         if (ivol(9,jvol).eq.35) then
  197.           ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol),
  198.      >                    ivol(3,jvol),ivol(5,jvol))/6.
  199.      >               +vol(ivol(1,jvol),ivol(3,jvol),
  200.      >                    ivol(4,jvol),ivol(5,jvol))/6.
  201.         endif
  202.         if (ivol(9,jvol).eq.30) then
  203.           ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol),
  204.      >                    ivol(3,jvol),ivol(4,jvol))/6.
  205.      >               +vol(ivol(2,jvol),ivol(3,jvol),
  206.      >                    ivol(4,jvol),ivol(5,jvol))/6.
  207.      >               +vol(ivol(3,jvol),ivol(5,jvol),
  208.      >                    ivol(6,jvol),ivol(4,jvol))/6.
  209.         endif
  210.         if (ivol(9,jvol).eq.20) then
  211.           ymesu=ymesu-vol(ivol(1,jvol),ivol(3,jvol),
  212.      >                    ivol(6,jvol),ivol(8,jvol))/6.
  213.      >               -vol(ivol(5,jvol),ivol(6,jvol),
  214.      >                    ivol(8,jvol),ivol(1,jvol))/6.
  215.      >               -vol(ivol(2,jvol),ivol(6,jvol),
  216.      >                    ivol(1,jvol),ivol(3,jvol))/6.
  217.           ymesu=ymesu-vol(ivol(7,jvol),ivol(8,jvol),
  218.      >                    ivol(6,jvol),ivol(3,jvol))/6.
  219.      >               -vol(ivol(4,jvol),ivol(1,jvol),
  220.      >                    ivol(8,jvol),ivol(3,jvol))/6.
  221.         endif
  222. 222    continue
  223.        nvoly=nvol
  224.        if (ymesu.gt.xmesu*1.01) goto 340
  225.        if (ierr.ne.0) goto 340
  226. *  AJUSTEMENT EVENTUEL DES DIMENSIONS DES TABLEAUX
  227.        IF (NFCMAX+250.GE.NFTOT) THEN
  228.         NFTOT=NFTOT+300
  229.         SEGADJ NFC,NFV,IFUT,IFAT
  230.        ENDIF
  231.        IF (NVOL+10.GE.NVTOT) THEN
  232.         NVTOT=NVTOT+50
  233.         SEGADJ IVOL
  234.        ENDIF
  235.        IF (NPTMAX+50.GE.NPTOT) THEN
  236.         NPTOT=NPTOT+100
  237.         SEGADJ NPF,XYZ
  238.        ENDIF
  239. *  nouvelle methode de calcul de la taille locale
  240.        DO 221 I=NPTDEB+1,NPTMAX
  241.         call vcrit(i)
  242. 221    CONTINUE
  243.        NPTDEB=NPTMAX
  244.        IGAGNE=0
  245.        IF (NFACET.EQ.0) GOTO 370
  246.        NFCMA=NFCMAX
  247. C
  248. C
  249. C      RECHERCHE DES DIEDRES
  250. C      ---------------------
  251. *  FAIRE ICI LE MENAGE DANS NARET(de temps en temps)
  252.       IF (ICRTS.GE.100) THEN
  253.        NPTDIS=NPTMAX
  254.        ICRTS=0
  255.        IVA=0
  256.        DO 285 I=1,NARET
  257.        if (IIGARD(I).LE.0) goto 285
  258.        IF1=IF1GAR(I)
  259.        IF (IFAT(IF1).EQ.0) GOTO 285
  260.        IF2=IF2GAR(I)
  261.        IF (IFAT(IF2).EQ.0) GOTO 285
  262.        IVA=IVA+1
  263.          II=IIGARD(I)
  264.          NPTDIS=MIN(NPTDIS,II)
  265.          IIGARD(IVA)=II
  266.          JJ=JJGARD(I)
  267.          NPTDIS=MIN(NPTDIS,JJ)
  268.          JJGARD(IVA)=JJ
  269.          IF1=IF1GAR(I)
  270.          IF1GAR(IVA)=IF1
  271.          IF2=IF2GAR(I)
  272.          IF2GAR(IVA)=IF2
  273.          ANGAR(IVA)=ANGAR(I)
  274.          IF (IVA.NE.1) THEN
  275.              ANGMA(IVA)=MAX(ANGAR(IVA),ANGMA(IVA-1))
  276.          ELSE
  277.              ANGMA(1)=ANGAR(1)
  278.          ENDIF
  279.  285   CONTINUE
  280. *      write (6,*) ' demait retassement effectue ',naret,iva
  281.        NARET=IVA
  282.       ENDIF
  283.        DO 290 IF1=NFCTRT+1,NFCMAX
  284.        IF (IFAT(IF1).EQ.0) GOTO 290
  285.        NBD=4
  286.        IF (NFC(4,IF1).EQ.0) NBD=3
  287.        DO 292 IC=1,NBD
  288.         IC1=IC-1
  289.         IF (IC1.EQ.0) IC1=NBD
  290.         II=NFC(IC1,IF1)
  291.         JJ=NFC(IC,IF1)
  292.         DO 294 I=1,40
  293.          IF2=NPF(I,II)
  294.          IF (IF2.EQ.0) GOTO 292
  295.          IF (IF2.GE.IF1) GOTO 294
  296.          IF (IPRED(IF2,II).NE.JJ) GOTO 294
  297. C
  298. C      COMMENT SONT LES ELEMENTS DU DIEDRE ?
  299. C      -------------------------------------
  300. C
  301.        ANGLL=TETA(IF1,IF2,II,JJ)
  302. C
  303. *  pour gagner du temps
  304.        if (angll.le.-1d4) goto 294
  305. *      if (if1.le.nfcori.or.if2.le.nfcori) angll=angll+1d6
  306.        NARET=NARET+1
  307.        IF (NARET.GT.NTTRAV) THEN
  308.          NTTRAV=NARET+20
  309.          SEGADJ TRAV
  310.        ENDIF
  311.         IIGARD(NARET)=II
  312.         JJGARD(NARET)=JJ
  313.         IF1GAR(NARET)=IF1
  314.         IF2GAR(NARET)=IF2
  315.         ANGAR(NARET)=ANGLL
  316.          IF (NARET.NE.1) THEN
  317.              ANGMA(NARET)=MAX(ANGAR(NARET),ANGMA(NARET-1))
  318.          ELSE
  319.              ANGMA(1)=ANGAR(1)
  320.          ENDIF
  321. 294    CONTINUE
  322. 292    CONTINUE
  323. 290    CONTINUE
  324.        NFCTRT=NFCMAX
  325. *
  326. *  ON COMMENCE PAR L'ANGLE LE PLUS FERME
  327. *
  328.  315   CONTINUE
  329.        ANGMAX=-1.E30
  330.        IOK=0
  331.        DO 310 I=NARET,1,-1
  332.         IF(ANGMAX.GE.ANGMA(I)) GOTO 311
  333.         IF(IIGARD(I).LE.0) GOTO 310
  334.         IF(ANGMAX.GE.ANGAR(I)) GOTO 310
  335.         IOK=I
  336. *       if (angar(i).gt.0d0) ANGMAX=ANGAR(I)*0.99999D0
  337. *       if (angar(i).lt.0d0) ANGMAX=ANGAR(I)*1.00001D0
  338.         angmax=angar(i)-1d-6
  339.  
  340.  310   CONTINUE
  341.  311   CONTINUE
  342.        IF (IOK.EQ.0) GOTO 320
  343.         II=IIGARD(IOK)
  344.         JJ=JJGARD(IOK)
  345.         IF1=IF1GAR(IOK)
  346.         IF2=IF2GAR(IOK)
  347.         IIGARD(IOK)=-II
  348.         ICRTS=ICRTS+1
  349.         IF (IFAT(IF1).EQ.0) GOTO 315
  350.         IF (IFAT(IF2).EQ.0) GOTO 315
  351.         IGAGNE=0
  352. *      write (6,*) 'traitement diedre ',ii,jj,if1,if2,angmax
  353.        idiac=5
  354.  313   continue
  355. *  on essaie d'abord de faire les hexaedres
  356.        if (ipass.eq.2) then
  357.         IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).NE.0)
  358.      #             CALL hexa(II,JJ,IF1,IF2,IGAGNE)
  359.          GOTO 312
  360.        else
  361.         IF (NFC(4,IF1).EQ.0.AND.NFC(4,IF2).EQ.0)
  362.      #             CALL CONS33(II,JJ,IF1,IF2,IGAGNE,0)
  363.         IF (IGAGNE.EQ.1) GOTO 312
  364.         IF (NFC(4,IF1).EQ.0.AND.NFC(4,IF2).NE.0)
  365.      #             CALL CONS34(II,JJ,IF1,IF2,IGAGNE)
  366.         IF (IGAGNE.EQ.1) GOTO 312
  367.         IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).EQ.0)
  368.      #             CALL CONS34(JJ,II,IF2,IF1,IGAGNE)
  369.         IF (IGAGNE.EQ.1) GOTO 312
  370.         IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).NE.0)
  371.      #             CALL CONS44(II,JJ,IF1,IF2,IGAGNE)
  372.         IF (IGAGNE.EQ.1) GOTO 312
  373.        endif
  374. *       write (6,*) ' demait relachement de diacrd'
  375.         diacrd=1-0.3*(1-diacrd)
  376.         diacri=diacrd
  377.         cfacet=cfacet*1.5
  378.         faccri=faccri*2
  379.         cdist=0.085
  380.         tetrl=tetrl*1.3
  381.         idiac=idiac-1
  382.         if (idiac.ne.0) goto 313
  383.         diacrd=diacre
  384.         diacri=diacrd
  385.         cfacet=cfacei
  386.         faccri=16
  387.         cdist=0.125
  388.         tetrl=2.75
  389.         IF (IVERB.EQ.1) write (6,*) ' demait echec traitement diedre',
  390.      >                                           ii,jj,if1,if2,angmax
  391.         goto 315
  392.  312   CONTINUE
  393.         diacrd=diacre
  394.         diacri=diacrd
  395.         cfacet=cfacei
  396.         faccri=16
  397.         cdist=0.125
  398.         tetrl=2.75
  399.         if (ipass.eq.-2) then
  400.          ipass=-1
  401. *        cfacet=6
  402. *        cdist=0.125
  403. *        tetrl=2.75
  404. *        faccri=16
  405. *        volcri=0.2
  406. *       diacri=0.90
  407. *       diacrd=0.92
  408.         endif
  409.        GOTO 220
  410.  320  CONTINUE
  411.        IF (ipass.eq.2) then
  412.         ipass=1
  413.         IF (IVERB.EQ.1) write (6,*) ' fin generation de cube '
  414. *  on remet les compteurs a zero
  415.         nfctrt=0
  416.         naret =0
  417.         goto 220
  418.        endif
  419.        IF (ipass.eq.1) then
  420.         ipass=0
  421.         IF (IVERB.EQ.1) write (6,*) ' strategie finale 1'
  422.         cfacet=16
  423.         cdist=0.085
  424.         tetrl=9
  425.         diacri=0.95
  426.         diacrd=0.95
  427.          faccri=64
  428.          volcri=0.01
  429. *  on remet les compteurs a zero
  430.         nfctrt=0
  431.         naret =0
  432.         goto 220
  433.        endif
  434.        IF (ipass.eq.0) then
  435.         ipass=-1
  436.         IF (IVERB.EQ.1) write (6,*) ' strategie finale 2'
  437.         cfacet=100
  438.         cdist=0.050
  439.         tetrl=9
  440.         diacri=0.99
  441.         diacrd=0.99
  442.          faccri=64
  443.          volcri=0.01
  444. *  on remet les compteurs a zero
  445.         nfctrt=0
  446.         naret =0
  447.         goto 220
  448.        endif
  449.        IF (ipass.eq.-1) then
  450.         ipass=-2
  451.         IF (IVERB.EQ.1) write (6,*) ' strategie finale 3'
  452.         cfacet=1000
  453.         cdist=0.005
  454.         tetrl=9
  455.         diacri=0.999
  456.         diacrd=0.999
  457.          faccri=64
  458.          volcri=0.01
  459. *  on remet les compteurs a zero
  460.         nfctrt=0
  461.         naret =0
  462.         goto 220
  463.        endif
  464.  340   continue
  465.        DO 444 I=1,NFCMAX
  466.         IF (IFAT(I).EQ.0) GOTO 444
  467.         IF (IVERB.EQ.1) WRITE (6,*) ' FACETTE RESTANTE IF ',I,
  468.      *        NFC(1,I),NFC(2,I),NFC(3,I),NFC(4,I)
  469.  444  CONTINUE
  470.  4001  CONTINUE
  471.        IF (NFACET.NE.0) CALL ERREUR(27)
  472.  370   CONTINUE
  473.        IF (IVERB.EQ.1) WRITE (6,*) ' volume du maillage   ',ymesu
  474.        IF (IVERB.EQ.1) write (6,*) ' nb facette ',nfacet
  475.        SEGSUP TRAV
  476.         CALL OPTVOL
  477.        RETURN
  478. C      FIN DU PROGRAMME PRINCIPAL
  479.        END
  480.  
  481.  
  482.  
  483.  
  484.  
  485.  
  486.  
  487.  
  488.  

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