Télécharger la plaquette Cast3M


fr en

Télécharger quali6.eso

Retour à la liste

Numérotation des lignes : 

quali6

  1. C QUALI6    SOURCE    GOUNAND   21/04/01    21:15:07     10933          
  2.       SUBROUTINE QUALI6(MELEMX,IELDEB,IELFIN,IMET,IMOMET,XDENS,KCMETR
  3.      $     ,KPVIRT,XVTOL,MLREEL,NDQC)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5.       IMPLICIT INTEGER (I-N)
  6. C***********************************************************************
  7. C NOM         : QUALI6
  8. C DESCRIPTION : Etant donné un maillage volumique simple MELEMX, 
  9. C               on construit la qualité de chacun de ses éléments 
  10. C               dans un listreel MLREEL.
  11. C               MELEMX est supposé actif.
  12. C               MLREEL est rendu actif.
  13. C      
  14. C     Par rapport à quali2, on utilise xvtol pour mettre
  15. C     le volume d'un élément à 0 s'il est petit.
  16. C     Ceci est important car on utilise le signe pour dégrader la
  17. C     qualité d'un élément (-1)
  18. C
  19. C     Par rapport à quali3, MELEME devient un MELEMX, MLREEL est un
  20. C     segment déjà existant et on introduit les éléments de début et de
  21. C     fin IELDEB et IELFIN qui servent pour MELEMX ET MLREEL (qui sont
  22. C     supposés de même dimension cf. trlver.eso.)
  23. *     Pour MLREEL, comme on ne calcule pas la qualité des éléments
  24. *     contenant le noeud virtuel, on a NDQC qui dit le nombre de qualité
  25. *     calculés (IELDEB sert donc pour MELEMX et MLREEL mais IELFIN
  26. *     uniquement pour MELEMX et NDQC uniquement pour MLREEL).
  27. C      
  28. *     Par rapport à quali5, on essaie de simplifier et de regrouper les
  29. *     cas avec, sans métrique + un peu de ménage
  30. C     
  31. C     
  32. C
  33. C LANGAGE     : ESOPE
  34. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  35. C               mél : gounand@semt2.smts.cea.fr
  36. C***********************************************************************
  37. C VERSION    : v1, 27/11/2017, version initiale
  38. C HISTORIQUE : v1, 27/11/2017, création
  39. C HISTORIQUE :
  40. C HISTORIQUE :
  41. C***********************************************************************
  42. -INC CCGEOME
  43. -INC PPARAM
  44. -INC CCOPTIO
  45. -INC CCREEL
  46. -INC SMCOORD
  47. -INC TMATOP1      
  48. *-INC SMETRIQ
  49.       POINTEUR KCMETR.METRIQ
  50. *-INC SMELEMX
  51. -INC SMLREEL
  52.       PARAMETER(NMET=6)
  53.       DIMENSION XMET(NMET)
  54.       DIMENSION XMET2(2,2)
  55.       DIMENSION XMET3(3,3)
  56.       REAL*8 DXI(3)
  57. * Statement functions
  58. *      DISTA(A,B,C,D)=SQRT((A-C)*(A-C)+(B-D)*(B-D))
  59. *      DISTB(A,B,C,D,E,F)=SQRT((A-D)*(A-D)+(B-E)*(B-E)+(C-F)*(C-F))
  60. *      DISTA(A,B)=SQRT(A*A+B*B)
  61. *      DISTB(A,B,C)=SQRT(A*A+B*B+C*C)
  62. *      DETTRI(A11,A12,A21,A22)=A11*A22-A21*A12
  63.       DETTET(A11,A12,A13,A21,A22,A23,A31,A32,A33)=
  64.      &  A11*(A22*A33-A23*A32)
  65.      &       +A12*(A23*A31-A21*A33)
  66.      &       +A13*(A21*A32-A22*A31)
  67. *     
  68. * Executable statements
  69. *
  70.       istneg=0
  71.       xstmax=-1.*xgrand
  72.       xstmin=xgrand
  73.       ISIGN=0
  74.       NDQC=0
  75.       IDIMP1=IDIM+1
  76. *      NBNN=NUMX(/1)
  77.       NBNN=NNCOU
  78. *      NBELEM=NUMX(/2)
  79. *?      JG=PROG(/1)
  80. * Trop de verif nuit...
  81. *      IF (NBNN.NE.IDIMP1.OR.JG.NE.NBELEM) THEN
  82.       IF (NBNN.NE.IDIMP1) THEN
  83.          CALL ERREUR(5)
  84.          RETURN
  85.       ENDIF
  86.       IF
  87.      $     (.NOT.(IELDEB.GE.1.AND.IELFIN.GE.IELDEB.AND.NLCOU.GE.IELFIN
  88.      $     .AND.NUMX(/2).GE.NLCOU))  THEN
  89.          WRITE(IOIMP,*) 'coucou quali6'
  90.          write(ioimp,*) 'IELDEB=',IELDEB
  91.          write(ioimp,*) 'IELFIN=',IELFIN
  92.          write(ioimp,*) 'NLCOU=',NLCOU
  93.          write(ioimp,*) 'NUM2=',NUMX(/2)
  94.          write(ioimp,*) 'MELEMX=',MELEMX
  95.          CALL ECMELX(MELEMX,0)
  96.          CALL ERREUR(5)
  97.          RETURN
  98.       ENDIF
  99. *     NARET=(NBNN-1)*NBNN/2
  100. c$$$      nbnn=2
  101. c$$$      NARET=(NBNN-1)*(NBNN/2)
  102. c$$$      WRITE(IOIMP,*) 'NBNN,NARET=',NBNN,NARET
  103. c$$$      nbnn=3
  104. c$$$      NARET=(NBNN-1)*(NBNN/2)
  105. c$$$      WRITE(IOIMP,*) 'NBNN,NARET=',NBNN,NARET
  106. c$$$      nbnn=4
  107. c$$$      NARET=(NBNN-1)*(NBNN/2)
  108. c$$$      WRITE(IOIMP,*) 'NBNN,NARET=',NBNN,NARET
  109. c$$$      stop 16
  110. *     2017/08/28 La ligne ci-dessus peut-elle induire des bugs ?
  111. *       Après test, non mais c'est chaud quand même...      
  112.       NARET=((NBNN-1)*NBNN)/2
  113. * Coefficient de normalisation
  114.       I=IDIM
  115.       XCOF=DFACT(I)*(SQRT((DBLE(2**I))/(DBLE(I+1))))
  116. *
  117.       IF (IMET.EQ.1) XMETD=1.D0/DENSIT
  118.       IF (IMET.EQ.2) XMETD=1.D0/XDENS
  119.       IF (IMET.EQ.3) NFMET=1
  120.       IF (IMET.EQ.4) NFMET=IDIM*(IDIM+1)/2
  121.  
  122. *     DO 10 IBELEM=1,NUMX(/2)
  123.       DO 10 IBELEM=IELDEB,IELFIN
  124. *     WRITE(IOIMP,*) 'IBELEM=',IBELEM
  125. *     Calcul de la métrique moyenne et de son déterminant suivant les
  126. *     cas 
  127.          if (imet.gt.0) then
  128.             IF (IMET.GE.1.AND.IMET.LE.3) THEN
  129.                IF (IMET.EQ.3) THEN
  130.                   YDENS=0.D0
  131.                   YDMIN=xgrand
  132.                   YDMAX=-xgrand
  133.                   DO I=1,IDIMP1
  134.                      INO=NUMX(I,IBELEM)
  135.                      IF (KPVIRT.NE.0) THEN
  136.                         IF (KPVIRT.EQ.INO) GOTO 10
  137.                      ENDIF
  138. *!
  139. *! MODIF !
  140. *!
  141.                      if (imomet.eq.1) then 
  142.                         YDENS=YDENS+LOG(KCMETR.XIN(1,INO))
  143.                         YDMIN=MIN(YDMIN,LOG(KCMETR.XIN(1,INO)))
  144.                         YDMAX=MAX(YDMAX,LOG(KCMETR.XIN(1,INO)))
  145.                      else
  146.                         YDENS=YDENS+KCMETR.XIN(1,INO)
  147.                         YDMIN=MIN(YDMIN,KCMETR.XIN(1,INO))
  148.                         YDMAX=MAX(YDMAX,KCMETR.XIN(1,INO))
  149.                      endif
  150.                   ENDDO
  151. *!
  152. *! MODIF !
  153. *!
  154.                   if (imomet.eq.1) then
  155.                      YDENS=EXP(YDENS/IDIMP1)
  156.                      YDMIN=EXP(YDMIN)
  157.                      YDMAX=EXP(YDMAX)
  158.                   else
  159.                      YDENS=YDENS/IDIMP1
  160.                   endif
  161. *                  write(ioimp,*) 'ydens,ydmin,ydmax=',ydens,ydmin,ydmax
  162. *                  YDENS=YDMIN
  163.                   if (.false.) then
  164. *     ydens2=1.D0/SQRT(YDENS)
  165.                      ydens2=YDENS
  166.                      xxxx=abs(ydens2-densit)/densit
  167.                      if (xxxx.gt.xzprec) then
  168.                         write(ioimp,*) 'ydens2,densit,xxxx=',ydens2
  169.      $                       ,densit,xxxx
  170.                         call erreur(5)
  171.                         return
  172.                      endif
  173.                   endif
  174. *                  XMETD=1.D0/YDENS
  175.                   XMETD=SQRT(YDENS)
  176.                ENDIF
  177.                XDETMD=XMETD**IDIM
  178.             ELSEIF (IMET.EQ.4) THEN
  179.                DO J=1,NFMET
  180.                   XMET(J)=0.D0
  181.                ENDDO
  182.                DO I=1,IDIMP1
  183.                   INO=NUMX(I,IBELEM)
  184.                   IF (KPVIRT.NE.0) THEN
  185.                      IF (KPVIRT.EQ.INO) GOTO 10
  186.                   ENDIF
  187.                   DO J=1,NFMET
  188.                      XMET(J)=XMET(J)+KCMETR.XIN(J,INO)
  189.                   ENDDO
  190.                ENDDO
  191.                DO J=1,NFMET
  192.                   XMET(J)=XMET(J)/IDIMP1
  193.                ENDDO
  194.                IF (IDIM.EQ.2) THEN
  195.                   XMET2(1,1)=XMET(1)
  196.                   XMET2(2,1)=XMET(2)
  197.                   XMET2(1,2)=XMET(2)
  198.                   XMET2(2,2)=XMET(3)
  199.                   XDETMD=XMET2(1,1)*XMET2(2,2)-XMET2(2,1)*XMET2(1,2)
  200.                   XDETMD=SQRT(XDETMD)
  201.                ELSEIF (IDIM.EQ.3) THEN
  202.                   XMET3(1,1)=XMET(1)
  203.                   XMET3(2,1)=XMET(2)
  204.                   XMET3(1,2)=XMET(2)
  205.                   XMET3(2,2)=XMET(3)
  206.                   XMET3(3,1)=XMET(4)
  207.                   XMET3(1,3)=XMET(4)
  208.                   XMET3(3,2)=XMET(5)
  209.                   XMET3(2,3)=XMET(5)
  210.                   XMET3(3,3)=XMET(6)
  211.                   XDETMD=DETTET(XMET3(1,1),XMET3(1,2),XMET3(1,3),XMET3(2
  212.      $                 ,1),XMET3(2,2),XMET3(2,3),XMET3(3,1),XMET3(3,2)
  213.      $                 ,XMET3(3,3))   
  214.                   XDETMD=SQRT(XDETMD)
  215.                ELSE
  216.                   WRITE(IOIMP,*) 'quali6 idim=',IDIM
  217.                   INTERR(1)=IDIM
  218.                   CALL ERREUR(709)
  219.                   RETURN
  220.                ENDIF
  221.             ELSE
  222.                WRITE(IOIMP,*) 'quali6 imet=',IMET
  223.                CALL ERREUR(5)
  224.                RETURN
  225.             ENDIF
  226.          endif
  227. * Calcul des volumes
  228. * Volume d'un triangle
  229.          IF (IDIM.EQ.2) THEN 
  230.             I0=NUMX(1,IBELEM)
  231.             I1=NUMX(2,IBELEM)
  232.             I2=NUMX(3,IBELEM)
  233.             IF (KPVIRT.NE.0) THEN
  234.                IF (KPVIRT.EQ.I0.OR.KPVIRT.EQ.I1.OR.KPVIRT.EQ.I2) 
  235.      $              goto 10
  236.             ENDIF
  237.             IP0=(I0-1)*IDIMP1
  238.             IP1=(I1-1)*IDIMP1
  239.             IP2=(I2-1)*IDIMP1
  240.             X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  241.             Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  242.             X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  243.             Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  244. *     XVOL=ABS(DETTRI(X10,Y10,X20,Y20))/2.D0
  245. *     XVOL=ABS(X10*Y20-X20*Y10)/2.D0
  246.             XVOLO=(X10*Y20-X20*Y10)/2.D0
  247. *     sigvol=sign(1.d0,xvolo)
  248. *     Volume d'un tétraèdre
  249.          ELSEIF (IDIM.EQ.3) THEN
  250.             I0=NUMX(1,IBELEM)
  251.             I1=NUMX(2,IBELEM)
  252.             I2=NUMX(3,IBELEM)
  253.             I3=NUMX(4,IBELEM)
  254.             IF (KPVIRT.NE.0) THEN
  255.                IF (KPVIRT.EQ.I0.OR.KPVIRT.EQ.I1.OR.KPVIRT.EQ.I2
  256.      $              .OR.KPVIRT.EQ.I3) goto 10
  257.             ENDIF
  258.             IP0=(I0-1)*IDIMP1
  259.             IP1=(I1-1)*IDIMP1
  260.             IP2=(I2-1)*IDIMP1
  261.             IP3=(I3-1)*IDIMP1
  262.             X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  263.             Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  264.             Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
  265.             X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  266.             Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  267.             Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
  268.             X30=XCOOR(IP3+1)-XCOOR(IP0+1)
  269.             Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
  270.             Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
  271. *     XVOL=ABS(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  272. *     $              /6.D0
  273.             XVOLO=(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  274.      $           /6.D0
  275.          ENDIF
  276. * Déterminant de la métrique
  277.          if (imet.gt.0) then
  278.             XVOLO=XVOLO*XDETMD
  279.          endif            
  280.          sigvol=sign(1.d0,xvolo+xvtol)
  281.          if (sigvol.lt.0.d0) then
  282.             istneg=istneg+1
  283.             xstmax=max(xstmax,xvolo)
  284.             xstmin=min(xstmin,xvolo)
  285.          endif
  286.          xvol=abs(xvolo)
  287. *     Calcul des longueurs
  288.          XLAR=0.D0
  289.          XLAR2=0.D0
  290.          IMEXP=IDIM
  291.             XLAMAX=XPETIT**(1.D0/IMEXP)
  292.             XLAMIN=XGRAND**(1.D0/IMEXP)
  293.          DO IBNN=1,NBNN
  294.             DO JBNN=IBNN+1,NBNN
  295.                I0=NUMX(IBNN,IBELEM)
  296.                I1=NUMX(JBNN,IBELEM)
  297.                IP0=(I0-1)*IDIMP1
  298.                IP1=(I1-1)*IDIMP1
  299.                DO IIDIM=1,IDIM
  300.                   DXI(IIDIM)=XCOOR(IP1+IIDIM)-XCOOR(IP0+IIDIM)
  301.                   IF (IMET.GE.1.AND.IMET.LE.3) THEN
  302. *!
  303. *! MODIF !
  304. *!
  305. *!                     DXI(IIDIM)=DXI(IIDIM)*(1.D0/XMETD)
  306.                      DXI(IIDIM)=DXI(IIDIM)*XMETD
  307.                   ENDIF
  308.                ENDDO
  309.                DXLAR2=0.D0
  310.                IF (IMET.EQ.4) THEN
  311.                   IF (IDIM.EQ.2) THEN
  312.                      DO J=1,IDIM
  313.                         DO I=1,IDIM
  314.                            DXLAR2=DXLAR2+XMET2(I,J)*DXI(I)*DXI(J)
  315.                         ENDDO
  316.                      ENDDO
  317.                   ELSEIF (IDIM.EQ.3) THEN
  318.                      DO J=1,IDIM
  319.                         DO I=1,IDIM
  320.                            DXLAR2=DXLAR2+XMET3(I,J)*DXI(I)*DXI(J)
  321.                         ENDDO
  322.                      ENDDO
  323.                   ELSE
  324.                      WRITE(IOIMP,*) 'quali6 idim=',IDIM
  325.  
  326.                      CALL ERREUR(5)
  327.                      RETURN
  328.                   ENDIF
  329.                ELSE
  330.                   DO I=1,IDIM
  331.                      DXLAR2=DXLAR2+DXI(I)*DXI(I)
  332.                   ENDDO
  333.                ENDIF
  334.                DXLAR=SQRT(DXLAR2)
  335.                XLAR=XLAR+DXLAR
  336.                XLAR2=XLAR2+DXLAR2
  337. *               XLAR2=XLAR2+(DXLAR**2)
  338.                XLAMAX=MAX(XLAMAX,DXLAR)
  339.                XLAMIN=MIN(XLAMIN,DXLAR)
  340.             ENDDO
  341.          ENDDO
  342.          XLAMOY=XLAR/NARET
  343.          if (imet.eq.0) then
  344.             XQUAL=XVOL/(XLAMOY**IDIM)
  345.             XQUALN=XQUAL*XCOF
  346.          else            
  347.             XLAMO2=SQRT(XLAR2*(2.D0/(DBLE(IDIM)*DBLE(IDIM+1))))
  348. **tst 2017/08/29 XLAMO3 ne marche pas bien            
  349. **     XLAMO3=SQRT(XLAMAX*XLAMIN)
  350.             XLAD=XLAMO2**IMEXP
  351. *            XLAD=XLAMOY**IDIM
  352. *
  353.             XLAMAD=XLAMAX**IMEXP
  354.             XLAMID=XLAMIN**IMEXP
  355.             XQUAL=XVOL/XLAD
  356.             XQUALN=XQUAL*XCOF
  357. *            PROG(IBELEM)=XQUALN
  358. *     PROG(IBELEM)=MIN(XQUALN,XLAD,1.D0/XLAD)
  359. *      Je n'aime pas trop ce critère            
  360. * 2017/08/28            PROG(**)=MIN(XQUALN,XLAD,1.D0/XLAD)
  361. *     write(ioimp,*) 'xqualn,xqual2=',xqualn,xqual2
  362. *            write(ioimp,*) 'xlad=',xlad
  363. *      xqual2=min(XLAD,1.D0/XLAD)
  364.             xqualn=MIN(XQUALN,XLAD,1.D0/XLAD)
  365. *            xqualn=MIN(XQUALN,XQUAL2)
  366. *            write(ioimp,*) 'xqualn,xlamax,xlamin,xlamad,xlamid=',xqualn
  367. *     $           ,xlamax,xlamin,min(XLAMAD,1.D0/XLAMAD),min(XLAMID,1.D0
  368. *     $           /XLAMID)
  369. *            xqual2=min(XLAMAD,1.D0/XLAMAD,XLAMID,1.D0/XLAMID)
  370. *            xqual2=min(XLAMID,1.D0/XLAMID)
  371. *            write(ioimp,*) 'xqualn,xqual2=',xqualn,xqual2
  372. *!
  373. *! MODIF !
  374. *!
  375. *            xqualn=MIN(XQUALN,XLAMAD,1.D0/XLAMAD,XLAMID,1.D0/XLAMID)
  376. *            xqualn=MIN(XQUALN,XQUAL2)
  377.          endif
  378.          if (isign.eq.1) then
  379.             if (sigvol.lt.0.d0) xqualn=xqualn-1.d0
  380.          endif
  381. *     PROG(**)=XQUALN
  382. * faux            PROG(IBELEM)=XQUALN
  383.          PROG(IELDEB+NDQC)=XQUALN
  384.          NDQC=NDQC+1
  385.  10   CONTINUE
  386.  
  387. *      write(ioimp,*) 'quali6 : istneg=',istneg
  388. *      write(ioimp,*) 'xstmin',xstmin,'xstmax',xstmax
  389.  
  390.       RETURN
  391. *
  392. * formats      
  393. *
  394.   188  FORMAT (2X,12(A6,'=',1PG12.5,2X))
  395. *     
  396. * End of subroutine QUALI6
  397. *
  398.       END
  399.  
  400.  
  401.  

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