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   25/07/24    21:15:05     12334          
  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=-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. *     Ici on fait la moyenne arithmétique
  155. *     mais kcmetr contient le log du tenseur si imomet=1
  156.                      YDENS=YDENS+KCMETR.XIN(1,INO)
  157.                      YDMIN=MIN(YDMIN,KCMETR.XIN(1,INO))
  158.                      YDMAX=MAX(YDMAX,KCMETR.XIN(1,INO))
  159.                   ENDDO
  160.                   if (imomet.eq.1) then
  161.                      YDENS=EXP(YDENS/IDIMP1)
  162.                      YDMIN=EXP(YDMIN)
  163.                      YDMAX=EXP(YDMAX)
  164.                   else
  165.                      YDENS=YDENS/IDIMP1
  166.                   endif
  167. *                  write(ioimp,*) 'ydens,ydmin,ydmax=',ydens,ydmin,ydmax
  168. *                  YDENS=YDMIN
  169.                   if (.false.) then
  170. *     ydens2=1.D0/SQRT(YDENS)
  171.                      ydens2=YDENS
  172.                      xxxx=abs(ydens2-densit)/densit
  173.                      if (xxxx.gt.xzprec) then
  174.                         write(ioimp,*) 'ydens2,densit,xxxx=',ydens2
  175.      $                       ,densit,xxxx
  176.                         call erreur(5)
  177.                         return
  178.                      endif
  179.                   endif
  180. *                  XMETD=1.D0/YDENS
  181.                   XMETD=SQRT(YDENS)
  182.                ENDIF
  183.                XDETMD=XMETD**IDIM
  184.             ELSEIF (IMET.EQ.4) THEN
  185.                DO J=1,NFMET
  186.                   XMET(J)=0.D0
  187.                ENDDO
  188.                DO I=1,IDIMP1
  189.                   INO=NUMX(I,IBELEM)
  190.                   IF (KPVIRT.NE.0) THEN
  191.                      IF (KPVIRT.EQ.INO) GOTO 10
  192.                   ENDIF
  193.                   DO J=1,NFMET
  194.                      XMET(J)=XMET(J)+KCMETR.XIN(J,INO)
  195.                   ENDDO
  196.                ENDDO
  197.                DO J=1,NFMET
  198.                   XMET(J)=XMET(J)/IDIMP1
  199.                ENDDO
  200.                if (imomet.eq.1) then
  201.                   DO J=1,IDIM
  202.                      DO I=1,IDIM
  203.                         A(I,J)=XMET(IDXSYM(I,J,IDIM))
  204.                      ENDDO
  205.                   ENDDO
  206. *     Exponentielle du tenseur symétrique
  207.                   IOTENS=8
  208.                   IKAS=3
  209.                   CALL TENS2(IOTENS,IKAS,A,D,R)
  210.                   IF (IERR.NE.0) RETURN
  211.                   DO IFMET=1,NFMET
  212.                      I=INVSYM(1,IFMET,IDIM)
  213.                      J=INVSYM(2,IFMET,IDIM)
  214.                      XMET(IFMET)=R(I,J)
  215.                   ENDDO
  216.                endif
  217.                IF (IDIM.EQ.2) THEN
  218.                   XMET2(1,1)=XMET(1)
  219.                   XMET2(2,1)=XMET(2)
  220.                   XMET2(1,2)=XMET(2)
  221.                   XMET2(2,2)=XMET(3)
  222.                   XDETMD=XMET2(1,1)*XMET2(2,2)-XMET2(2,1)*XMET2(1,2)
  223. c$$$                  if (xdetmd.lt.0.d0) then
  224. c$$$                     write(ioimp,*) 'xdetmd=',xdetmd
  225. c$$$                     write(ioimp,*) 'xmet=',(xmet(i),i=1,3)
  226. c$$$                     write(ioimp,*) 'a=',((a(i,j),i=1,2),j=1,2)
  227. c$$$                     write(ioimp,*) 'r=',((r(i,j),i=1,2),j=1,2)
  228. c$$$                  endif
  229.                   XDETMD=SQRT(XDETMD)
  230.                ELSEIF (IDIM.EQ.3) THEN
  231.                   XMET3(1,1)=XMET(1)
  232.                   XMET3(2,1)=XMET(2)
  233.                   XMET3(1,2)=XMET(2)
  234.                   XMET3(2,2)=XMET(3)
  235.                   XMET3(3,1)=XMET(4)
  236.                   XMET3(1,3)=XMET(4)
  237.                   XMET3(3,2)=XMET(5)
  238.                   XMET3(2,3)=XMET(5)
  239.                   XMET3(3,3)=XMET(6)
  240.                   XDETMD=DETTET(XMET3(1,1),XMET3(1,2),XMET3(1,3),XMET3(2
  241.      $                 ,1),XMET3(2,2),XMET3(2,3),XMET3(3,1),XMET3(3,2)
  242.      $                 ,XMET3(3,3))
  243.                   XDETMD=SQRT(XDETMD)
  244.                ELSE
  245.                   WRITE(IOIMP,*) 'quali6 idim=',IDIM
  246.                   INTERR(1)=IDIM
  247.                   CALL ERREUR(709)
  248.                   RETURN
  249.                ENDIF
  250.             ELSE
  251.                WRITE(IOIMP,*) 'quali6 imet=',IMET
  252.                CALL ERREUR(5)
  253.                RETURN
  254.             ENDIF
  255.          endif
  256. * Calcul des volumes
  257. * Volume d'un triangle
  258.          IF (IDIM.EQ.2) THEN
  259.             I0=NUMX(1,IBELEM)
  260.             I1=NUMX(2,IBELEM)
  261.             I2=NUMX(3,IBELEM)
  262.             IF (KPVIRT.NE.0) THEN
  263.                IF (KPVIRT.EQ.I0.OR.KPVIRT.EQ.I1.OR.KPVIRT.EQ.I2)
  264.      $              goto 10
  265.             ENDIF
  266.             IP0=(I0-1)*IDIMP1
  267.             IP1=(I1-1)*IDIMP1
  268.             IP2=(I2-1)*IDIMP1
  269.             X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  270.             Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  271.             X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  272.             Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  273. *     XVOL=ABS(DETTRI(X10,Y10,X20,Y20))/2.D0
  274. *     XVOL=ABS(X10*Y20-X20*Y10)/2.D0
  275.             XVOLO=(X10*Y20-X20*Y10)/2.D0
  276. *     sigvol=sign(1.d0,xvolo)
  277. *     Volume d'un tétraèdre
  278.          ELSEIF (IDIM.EQ.3) THEN
  279.             I0=NUMX(1,IBELEM)
  280.             I1=NUMX(2,IBELEM)
  281.             I2=NUMX(3,IBELEM)
  282.             I3=NUMX(4,IBELEM)
  283.             IF (KPVIRT.NE.0) THEN
  284.                IF (KPVIRT.EQ.I0.OR.KPVIRT.EQ.I1.OR.KPVIRT.EQ.I2
  285.      $              .OR.KPVIRT.EQ.I3) goto 10
  286.             ENDIF
  287.             IP0=(I0-1)*IDIMP1
  288.             IP1=(I1-1)*IDIMP1
  289.             IP2=(I2-1)*IDIMP1
  290.             IP3=(I3-1)*IDIMP1
  291.             X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  292.             Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  293.             Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
  294.             X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  295.             Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  296.             Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
  297.             X30=XCOOR(IP3+1)-XCOOR(IP0+1)
  298.             Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
  299.             Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
  300. *     XVOL=ABS(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  301. *     $              /6.D0
  302.             XVOLO=(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  303.      $           /6.D0
  304.          ENDIF
  305. * Déterminant de la métrique
  306.          if (imet.gt.0) then
  307.             XVOLO=XVOLO*XDETMD
  308.          endif
  309.          sigvol=sign(1.d0,xvolo+xvtol)
  310.          if (sigvol.lt.0.d0) then
  311.             istneg=istneg+1
  312.             xstmax=max(xstmax,xvolo)
  313.             xstmin=min(xstmin,xvolo)
  314.          endif
  315.          xvol=abs(xvolo)
  316. *     Calcul des longueurs
  317.          XLAR=0.D0
  318.          XLAR2=0.D0
  319.          IMEXP=IDIM
  320.             XLAMAX=XPETIT**(1.D0/IMEXP)
  321.             XLAMIN=XGRAND**(1.D0/IMEXP)
  322.          DO IBNN=1,NBNN
  323.             DO JBNN=IBNN+1,NBNN
  324.                I0=NUMX(IBNN,IBELEM)
  325.                I1=NUMX(JBNN,IBELEM)
  326.                IP0=(I0-1)*IDIMP1
  327.                IP1=(I1-1)*IDIMP1
  328.                DO IIDIM=1,IDIM
  329.                   DXI(IIDIM)=XCOOR(IP1+IIDIM)-XCOOR(IP0+IIDIM)
  330.                   IF (IMET.GE.1.AND.IMET.LE.3) THEN
  331. *!
  332. *! MODIF !
  333. *!
  334. *!                     DXI(IIDIM)=DXI(IIDIM)*(1.D0/XMETD)
  335.                      DXI(IIDIM)=DXI(IIDIM)*XMETD
  336.                   ENDIF
  337.                ENDDO
  338.                DXLAR2=0.D0
  339.                IF (IMET.EQ.4) THEN
  340.                   IF (IDIM.EQ.2) THEN
  341.                      DO J=1,IDIM
  342.                         DO I=1,IDIM
  343.                            DXLAR2=DXLAR2+XMET2(I,J)*DXI(I)*DXI(J)
  344.                         ENDDO
  345.                      ENDDO
  346.                   ELSEIF (IDIM.EQ.3) THEN
  347.                      DO J=1,IDIM
  348.                         DO I=1,IDIM
  349.                            DXLAR2=DXLAR2+XMET3(I,J)*DXI(I)*DXI(J)
  350.                         ENDDO
  351.                      ENDDO
  352.                   ELSE
  353.                      WRITE(IOIMP,*) 'quali6 idim=',IDIM
  354.  
  355.                      CALL ERREUR(5)
  356.                      RETURN
  357.                   ENDIF
  358.                ELSE
  359.                   DO I=1,IDIM
  360.                      DXLAR2=DXLAR2+DXI(I)*DXI(I)
  361.                   ENDDO
  362.                ENDIF
  363.                DXLAR=SQRT(DXLAR2)
  364.                XLAR=XLAR+DXLAR
  365.                XLAR2=XLAR2+DXLAR2
  366. *               XLAR2=XLAR2+(DXLAR**2)
  367.                XLAMAX=MAX(XLAMAX,DXLAR)
  368.                XLAMIN=MIN(XLAMIN,DXLAR)
  369.             ENDDO
  370.          ENDDO
  371.          XLAMOY=XLAR/NARET
  372.          if (imet.eq.0) then
  373.             XQUAL=XVOL/(XLAMOY**IDIM)
  374.             XQUALN=XQUAL*XCOF
  375.          else
  376.             XLAMO2=SQRT(XLAR2*(2.D0/(DBLE(IDIM)*DBLE(IDIM+1))))
  377. **tst 2017/08/29 XLAMO3 ne marche pas bien
  378. **     XLAMO3=SQRT(XLAMAX*XLAMIN)
  379.             XLAD=XLAMO2**IMEXP
  380. *            XLAD=XLAMOY**IDIM
  381. *
  382.             XLAMAD=XLAMAX**IMEXP
  383.             XLAMID=XLAMIN**IMEXP
  384.             XQUAL=XVOL/XLAD
  385.             XQUALN=XQUAL*XCOF
  386. *            PROG(IBELEM)=XQUALN
  387. *     PROG(IBELEM)=MIN(XQUALN,XLAD,1.D0/XLAD)
  388. *      Je n'aime pas trop ce critère
  389. * 2017/08/28            PROG(**)=MIN(XQUALN,XLAD,1.D0/XLAD)
  390. *     write(ioimp,*) 'xqualn,xqual2=',xqualn,xqual2
  391. *            write(ioimp,*) 'xlad=',xlad
  392. *      xqual2=min(XLAD,1.D0/XLAD)
  393.             xqualn=MIN(XQUALN,XLAD,1.D0/XLAD)
  394. *            xqualn=MIN(XQUALN,XQUAL2)
  395. *            write(ioimp,*) 'xqualn,xlamax,xlamin,xlamad,xlamid=',xqualn
  396. *     $           ,xlamax,xlamin,min(XLAMAD,1.D0/XLAMAD),min(XLAMID,1.D0
  397. *     $           /XLAMID)
  398. *            xqual2=min(XLAMAD,1.D0/XLAMAD,XLAMID,1.D0/XLAMID)
  399. *            xqual2=min(XLAMID,1.D0/XLAMID)
  400. *            write(ioimp,*) 'xqualn,xqual2=',xqualn,xqual2
  401. *!
  402. *! MODIF !
  403. *!
  404. *            xqualn=MIN(XQUALN,XLAMAD,1.D0/XLAMAD,XLAMID,1.D0/XLAMID)
  405. *            xqualn=MIN(XQUALN,XQUAL2)
  406.          endif
  407.          if (isign.eq.1) then
  408.             if (sigvol.lt.0.d0) xqualn=xqualn-1.d0
  409.          endif
  410. *     PROG(**)=XQUALN
  411. * faux            PROG(IBELEM)=XQUALN
  412.          PROG(IELDEB+NDQC)=XQUALN
  413.          NDQC=NDQC+1
  414.  10   CONTINUE
  415.  
  416. *      write(ioimp,*) 'quali6 : istneg=',istneg
  417. *      write(ioimp,*) 'xstmin',xstmin,'xstmax',xstmax
  418.  
  419.       RETURN
  420. *
  421. * formats
  422. *
  423.   188  FORMAT (2X,12(A6,'=',1PG12.5,2X))
  424. *
  425. * End of subroutine QUALI6
  426. *
  427.       END
  428.  
  429.  

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