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

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