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   24/09/27    21:15:18     12019          
  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.       DIMENSION A(3,3),D(3),R(3,3)
  58.       DIMENSION IDXSYM(3,3,3)
  59.       DIMENSION INVSYM(2,6,3)
  60. *
  61. * Statement functions
  62. *      DISTA(A,B,C,D)=SQRT((A-C)*(A-C)+(B-D)*(B-D))
  63. *      DISTB(A,B,C,D,E,F)=SQRT((A-D)*(A-D)+(B-E)*(B-E)+(C-F)*(C-F))
  64. *      DISTA(A,B)=SQRT(A*A+B*B)
  65. *      DISTB(A,B,C)=SQRT(A*A+B*B+C*C)
  66. *      DETTRI(A11,A12,A21,A22)=A11*A22-A21*A12
  67.       DETTET(A11,A12,A13,A21,A22,A23,A31,A32,A33)=
  68.      &  A11*(A22*A33-A23*A32)
  69.      &       +A12*(A23*A31-A21*A33)
  70.      &       +A13*(A21*A32-A22*A31)
  71. *
  72.       DATA ((A(I,J),I=1,3),J=1,3) /9*0.D0/
  73.       DATA (D(I),I=1,3) /3*0.D0/
  74.       DATA ((R(I,J),I=1,3),J=1,3) /9*0.D0/
  75.       DATA (((IDXSYM(I,J,K),I=1,1),J=1,1),K=1,1) /1/
  76.       DATA (((IDXSYM(I,J,K),I=1,2),J=1,2),K=2,2) /1,2,2,3/
  77.       DATA (((IDXSYM(I,J,K),I=1,3),J=1,3),K=3,3) /1,2,4,2,3,5,4,5,6/
  78. *
  79.       DATA (((INVSYM(I,J,K),I=1,2),J=1,1),K=1,1) /1,1/
  80.       DATA (((INVSYM(I,J,K),I=1,2),J=1,3),K=2,2) /1,1,2,1,2,2/
  81.       DATA (((INVSYM(I,J,K),I=1,2),J=1,6),K=3,3)
  82.      $     /1,1,2,1,2,2,3,1,3,2,3,3/
  83. *
  84. * Executable statements
  85. *
  86.       istneg=0
  87.       xstmax=-1.*xgrand
  88.       xstmin=xgrand
  89.       ISIGN=0
  90.       NDQC=0
  91.       IDIMP1=IDIM+1
  92. *      NBNN=NUMX(/1)
  93.       NBNN=NNCOU
  94. *      NBELEM=NUMX(/2)
  95. *?      JG=PROG(/1)
  96. * Trop de verif nuit...
  97. *      IF (NBNN.NE.IDIMP1.OR.JG.NE.NBELEM) THEN
  98.       IF (NBNN.NE.IDIMP1) THEN
  99.          CALL ERREUR(5)
  100.          RETURN
  101.       ENDIF
  102.       IF
  103.      $     (.NOT.(IELDEB.GE.1.AND.IELFIN.GE.IELDEB.AND.NLCOU.GE.IELFIN
  104.      $     .AND.NUMX(/2).GE.NLCOU))  THEN
  105.          WRITE(IOIMP,*) 'coucou quali6'
  106.          write(ioimp,*) 'IELDEB=',IELDEB
  107.          write(ioimp,*) 'IELFIN=',IELFIN
  108.          write(ioimp,*) 'NLCOU=',NLCOU
  109.          write(ioimp,*) 'NUM2=',NUMX(/2)
  110.          write(ioimp,*) 'MELEMX=',MELEMX
  111.          CALL ECMELX(MELEMX,0)
  112.          CALL ERREUR(5)
  113.          RETURN
  114.       ENDIF
  115. *     NARET=(NBNN-1)*NBNN/2
  116. c$$$      nbnn=2
  117. c$$$      NARET=(NBNN-1)*(NBNN/2)
  118. c$$$      WRITE(IOIMP,*) 'NBNN,NARET=',NBNN,NARET
  119. c$$$      nbnn=3
  120. c$$$      NARET=(NBNN-1)*(NBNN/2)
  121. c$$$      WRITE(IOIMP,*) 'NBNN,NARET=',NBNN,NARET
  122. c$$$      nbnn=4
  123. c$$$      NARET=(NBNN-1)*(NBNN/2)
  124. c$$$      WRITE(IOIMP,*) 'NBNN,NARET=',NBNN,NARET
  125. c$$$      stop 16
  126. *     2017/08/28 La ligne ci-dessus peut-elle induire des bugs ?
  127. *       Après test, non mais c'est chaud quand même...
  128.       NARET=((NBNN-1)*NBNN)/2
  129. * Coefficient de normalisation
  130.       I=IDIM
  131.       XCOF=DFACT(I)*(SQRT((DBLE(2**I))/(DBLE(I+1))))
  132. *
  133.       IF (IMET.EQ.1) XMETD=1.D0/DENSIT
  134.       IF (IMET.EQ.2) XMETD=1.D0/XDENS
  135.       IF (IMET.EQ.3) NFMET=1
  136.       IF (IMET.EQ.4) NFMET=IDIM*(IDIM+1)/2
  137.  
  138. *     DO 10 IBELEM=1,NUMX(/2)
  139.       DO 10 IBELEM=IELDEB,IELFIN
  140. *     WRITE(IOIMP,*) 'IBELEM=',IBELEM
  141. *     Calcul de la métrique moyenne et de son déterminant suivant les
  142. *     cas
  143.          if (imet.gt.0) then
  144.             IF (IMET.GE.1.AND.IMET.LE.3) THEN
  145.                IF (IMET.EQ.3) THEN
  146.                   YDENS=0.D0
  147.                   YDMIN=xgrand
  148.                   YDMAX=-xgrand
  149.                   DO I=1,IDIMP1
  150.                      INO=NUMX(I,IBELEM)
  151.                      IF (KPVIRT.NE.0) THEN
  152.                         IF (KPVIRT.EQ.INO) GOTO 10
  153.                      ENDIF
  154. *!
  155. *! MODIF !
  156. *!
  157.                      if (imomet.eq.1) then
  158.                         YDENS=YDENS+LOG(KCMETR.XIN(1,INO))
  159.                         YDMIN=MIN(YDMIN,LOG(KCMETR.XIN(1,INO)))
  160.                         YDMAX=MAX(YDMAX,LOG(KCMETR.XIN(1,INO)))
  161.                      else
  162.                         YDENS=YDENS+KCMETR.XIN(1,INO)
  163.                         YDMIN=MIN(YDMIN,KCMETR.XIN(1,INO))
  164.                         YDMAX=MAX(YDMAX,KCMETR.XIN(1,INO))
  165.                      endif
  166.                   ENDDO
  167. *!
  168. *! MODIF !
  169. *!
  170.                   if (imomet.eq.1) then
  171.                      YDENS=EXP(YDENS/IDIMP1)
  172.                      YDMIN=EXP(YDMIN)
  173.                      YDMAX=EXP(YDMAX)
  174.                   else
  175.                      YDENS=YDENS/IDIMP1
  176.                   endif
  177. *                  write(ioimp,*) 'ydens,ydmin,ydmax=',ydens,ydmin,ydmax
  178. *                  YDENS=YDMIN
  179.                   if (.false.) then
  180. *     ydens2=1.D0/SQRT(YDENS)
  181.                      ydens2=YDENS
  182.                      xxxx=abs(ydens2-densit)/densit
  183.                      if (xxxx.gt.xzprec) then
  184.                         write(ioimp,*) 'ydens2,densit,xxxx=',ydens2
  185.      $                       ,densit,xxxx
  186.                         call erreur(5)
  187.                         return
  188.                      endif
  189.                   endif
  190. *                  XMETD=1.D0/YDENS
  191.                   XMETD=SQRT(YDENS)
  192.                ENDIF
  193.                XDETMD=XMETD**IDIM
  194.             ELSEIF (IMET.EQ.4) THEN
  195.                DO J=1,NFMET
  196.                   XMET(J)=0.D0
  197.                ENDDO
  198.                DO I=1,IDIMP1
  199.                   INO=NUMX(I,IBELEM)
  200.                   IF (KPVIRT.NE.0) THEN
  201.                      IF (KPVIRT.EQ.INO) GOTO 10
  202.                   ENDIF
  203.                   DO J=1,NFMET
  204.                      XMET(J)=XMET(J)+KCMETR.XIN(J,INO)
  205.                   ENDDO
  206.                ENDDO
  207.                DO J=1,NFMET
  208.                   XMET(J)=XMET(J)/IDIMP1
  209.                ENDDO
  210.                if (imomet.eq.1) then
  211.                   DO J=1,IDIM
  212.                      DO I=1,IDIM
  213.                         A(I,J)=XMET(IDXSYM(I,J,IDIM))
  214.                      ENDDO
  215.                   ENDDO
  216. *     Exponentielle du tenseur symétrique
  217.                   IOTENS=8
  218.                   IKAS=3
  219.                   CALL TENS2(IOTENS,IKAS,A,D,R)
  220.                   IF (IERR.NE.0) RETURN
  221.                   DO IFMET=1,NFMET
  222.                      I=INVSYM(1,IFMET,IDIM)
  223.                      J=INVSYM(2,IFMET,IDIM)
  224.                      XMET(IFMET)=R(I,J)
  225.                   ENDDO
  226.                endif
  227.                IF (IDIM.EQ.2) THEN
  228.                   XMET2(1,1)=XMET(1)
  229.                   XMET2(2,1)=XMET(2)
  230.                   XMET2(1,2)=XMET(2)
  231.                   XMET2(2,2)=XMET(3)
  232.                   XDETMD=XMET2(1,1)*XMET2(2,2)-XMET2(2,1)*XMET2(1,2)
  233. c$$$                  if (xdetmd.lt.0.d0) then
  234. c$$$                     write(ioimp,*) 'xdetmd=',xdetmd
  235. c$$$                     write(ioimp,*) 'xmet=',(xmet(i),i=1,3)
  236. c$$$                     write(ioimp,*) 'a=',((a(i,j),i=1,2),j=1,2)
  237. c$$$                     write(ioimp,*) 'r=',((r(i,j),i=1,2),j=1,2)
  238. c$$$                  endif
  239.                   XDETMD=SQRT(XDETMD)
  240.                ELSEIF (IDIM.EQ.3) THEN
  241.                   XMET3(1,1)=XMET(1)
  242.                   XMET3(2,1)=XMET(2)
  243.                   XMET3(1,2)=XMET(2)
  244.                   XMET3(2,2)=XMET(3)
  245.                   XMET3(3,1)=XMET(4)
  246.                   XMET3(1,3)=XMET(4)
  247.                   XMET3(3,2)=XMET(5)
  248.                   XMET3(2,3)=XMET(5)
  249.                   XMET3(3,3)=XMET(6)
  250.                   XDETMD=DETTET(XMET3(1,1),XMET3(1,2),XMET3(1,3),XMET3(2
  251.      $                 ,1),XMET3(2,2),XMET3(2,3),XMET3(3,1),XMET3(3,2)
  252.      $                 ,XMET3(3,3))
  253.                   XDETMD=SQRT(XDETMD)
  254.                ELSE
  255.                   WRITE(IOIMP,*) 'quali6 idim=',IDIM
  256.                   INTERR(1)=IDIM
  257.                   CALL ERREUR(709)
  258.                   RETURN
  259.                ENDIF
  260.             ELSE
  261.                WRITE(IOIMP,*) 'quali6 imet=',IMET
  262.                CALL ERREUR(5)
  263.                RETURN
  264.             ENDIF
  265.          endif
  266. * Calcul des volumes
  267. * Volume d'un triangle
  268.          IF (IDIM.EQ.2) THEN
  269.             I0=NUMX(1,IBELEM)
  270.             I1=NUMX(2,IBELEM)
  271.             I2=NUMX(3,IBELEM)
  272.             IF (KPVIRT.NE.0) THEN
  273.                IF (KPVIRT.EQ.I0.OR.KPVIRT.EQ.I1.OR.KPVIRT.EQ.I2)
  274.      $              goto 10
  275.             ENDIF
  276.             IP0=(I0-1)*IDIMP1
  277.             IP1=(I1-1)*IDIMP1
  278.             IP2=(I2-1)*IDIMP1
  279.             X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  280.             Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  281.             X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  282.             Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  283. *     XVOL=ABS(DETTRI(X10,Y10,X20,Y20))/2.D0
  284. *     XVOL=ABS(X10*Y20-X20*Y10)/2.D0
  285.             XVOLO=(X10*Y20-X20*Y10)/2.D0
  286. *     sigvol=sign(1.d0,xvolo)
  287. *     Volume d'un tétraèdre
  288.          ELSEIF (IDIM.EQ.3) THEN
  289.             I0=NUMX(1,IBELEM)
  290.             I1=NUMX(2,IBELEM)
  291.             I2=NUMX(3,IBELEM)
  292.             I3=NUMX(4,IBELEM)
  293.             IF (KPVIRT.NE.0) THEN
  294.                IF (KPVIRT.EQ.I0.OR.KPVIRT.EQ.I1.OR.KPVIRT.EQ.I2
  295.      $              .OR.KPVIRT.EQ.I3) goto 10
  296.             ENDIF
  297.             IP0=(I0-1)*IDIMP1
  298.             IP1=(I1-1)*IDIMP1
  299.             IP2=(I2-1)*IDIMP1
  300.             IP3=(I3-1)*IDIMP1
  301.             X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  302.             Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  303.             Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
  304.             X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  305.             Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  306.             Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
  307.             X30=XCOOR(IP3+1)-XCOOR(IP0+1)
  308.             Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
  309.             Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
  310. *     XVOL=ABS(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  311. *     $              /6.D0
  312.             XVOLO=(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  313.      $           /6.D0
  314.          ENDIF
  315. * Déterminant de la métrique
  316.          if (imet.gt.0) then
  317.             XVOLO=XVOLO*XDETMD
  318.          endif
  319.          sigvol=sign(1.d0,xvolo+xvtol)
  320.          if (sigvol.lt.0.d0) then
  321.             istneg=istneg+1
  322.             xstmax=max(xstmax,xvolo)
  323.             xstmin=min(xstmin,xvolo)
  324.          endif
  325.          xvol=abs(xvolo)
  326. *     Calcul des longueurs
  327.          XLAR=0.D0
  328.          XLAR2=0.D0
  329.          IMEXP=IDIM
  330.             XLAMAX=XPETIT**(1.D0/IMEXP)
  331.             XLAMIN=XGRAND**(1.D0/IMEXP)
  332.          DO IBNN=1,NBNN
  333.             DO JBNN=IBNN+1,NBNN
  334.                I0=NUMX(IBNN,IBELEM)
  335.                I1=NUMX(JBNN,IBELEM)
  336.                IP0=(I0-1)*IDIMP1
  337.                IP1=(I1-1)*IDIMP1
  338.                DO IIDIM=1,IDIM
  339.                   DXI(IIDIM)=XCOOR(IP1+IIDIM)-XCOOR(IP0+IIDIM)
  340.                   IF (IMET.GE.1.AND.IMET.LE.3) THEN
  341. *!
  342. *! MODIF !
  343. *!
  344. *!                     DXI(IIDIM)=DXI(IIDIM)*(1.D0/XMETD)
  345.                      DXI(IIDIM)=DXI(IIDIM)*XMETD
  346.                   ENDIF
  347.                ENDDO
  348.                DXLAR2=0.D0
  349.                IF (IMET.EQ.4) THEN
  350.                   IF (IDIM.EQ.2) THEN
  351.                      DO J=1,IDIM
  352.                         DO I=1,IDIM
  353.                            DXLAR2=DXLAR2+XMET2(I,J)*DXI(I)*DXI(J)
  354.                         ENDDO
  355.                      ENDDO
  356.                   ELSEIF (IDIM.EQ.3) THEN
  357.                      DO J=1,IDIM
  358.                         DO I=1,IDIM
  359.                            DXLAR2=DXLAR2+XMET3(I,J)*DXI(I)*DXI(J)
  360.                         ENDDO
  361.                      ENDDO
  362.                   ELSE
  363.                      WRITE(IOIMP,*) 'quali6 idim=',IDIM
  364.  
  365.                      CALL ERREUR(5)
  366.                      RETURN
  367.                   ENDIF
  368.                ELSE
  369.                   DO I=1,IDIM
  370.                      DXLAR2=DXLAR2+DXI(I)*DXI(I)
  371.                   ENDDO
  372.                ENDIF
  373.                DXLAR=SQRT(DXLAR2)
  374.                XLAR=XLAR+DXLAR
  375.                XLAR2=XLAR2+DXLAR2
  376. *               XLAR2=XLAR2+(DXLAR**2)
  377.                XLAMAX=MAX(XLAMAX,DXLAR)
  378.                XLAMIN=MIN(XLAMIN,DXLAR)
  379.             ENDDO
  380.          ENDDO
  381.          XLAMOY=XLAR/NARET
  382.          if (imet.eq.0) then
  383.             XQUAL=XVOL/(XLAMOY**IDIM)
  384.             XQUALN=XQUAL*XCOF
  385.          else
  386.             XLAMO2=SQRT(XLAR2*(2.D0/(DBLE(IDIM)*DBLE(IDIM+1))))
  387. **tst 2017/08/29 XLAMO3 ne marche pas bien
  388. **     XLAMO3=SQRT(XLAMAX*XLAMIN)
  389.             XLAD=XLAMO2**IMEXP
  390. *            XLAD=XLAMOY**IDIM
  391. *
  392.             XLAMAD=XLAMAX**IMEXP
  393.             XLAMID=XLAMIN**IMEXP
  394.             XQUAL=XVOL/XLAD
  395.             XQUALN=XQUAL*XCOF
  396. *            PROG(IBELEM)=XQUALN
  397. *     PROG(IBELEM)=MIN(XQUALN,XLAD,1.D0/XLAD)
  398. *      Je n'aime pas trop ce critère
  399. * 2017/08/28            PROG(**)=MIN(XQUALN,XLAD,1.D0/XLAD)
  400. *     write(ioimp,*) 'xqualn,xqual2=',xqualn,xqual2
  401. *            write(ioimp,*) 'xlad=',xlad
  402. *      xqual2=min(XLAD,1.D0/XLAD)
  403.             xqualn=MIN(XQUALN,XLAD,1.D0/XLAD)
  404. *            xqualn=MIN(XQUALN,XQUAL2)
  405. *            write(ioimp,*) 'xqualn,xlamax,xlamin,xlamad,xlamid=',xqualn
  406. *     $           ,xlamax,xlamin,min(XLAMAD,1.D0/XLAMAD),min(XLAMID,1.D0
  407. *     $           /XLAMID)
  408. *            xqual2=min(XLAMAD,1.D0/XLAMAD,XLAMID,1.D0/XLAMID)
  409. *            xqual2=min(XLAMID,1.D0/XLAMID)
  410. *            write(ioimp,*) 'xqualn,xqual2=',xqualn,xqual2
  411. *!
  412. *! MODIF !
  413. *!
  414. *            xqualn=MIN(XQUALN,XLAMAD,1.D0/XLAMAD,XLAMID,1.D0/XLAMID)
  415. *            xqualn=MIN(XQUALN,XQUAL2)
  416.          endif
  417.          if (isign.eq.1) then
  418.             if (sigvol.lt.0.d0) xqualn=xqualn-1.d0
  419.          endif
  420. *     PROG(**)=XQUALN
  421. * faux            PROG(IBELEM)=XQUALN
  422.          PROG(IELDEB+NDQC)=XQUALN
  423.          NDQC=NDQC+1
  424.  10   CONTINUE
  425.  
  426. *      write(ioimp,*) 'quali6 : istneg=',istneg
  427. *      write(ioimp,*) 'xstmin',xstmin,'xstmax',xstmax
  428.  
  429.       RETURN
  430. *
  431. * formats
  432. *
  433.   188  FORMAT (2X,12(A6,'=',1PG12.5,2X))
  434. *
  435. * End of subroutine QUALI6
  436. *
  437.       END
  438.  
  439.  

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