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   26/07/06    21:15:11     12592          
  2.       SUBROUTINE QUALI6(MELEMX,IELDEB,IELFIN,IMET,IMOMET,XDENS,KCMETR
  3.      $     ,NKPVIR,XVTOL,MLREEL,NDQC,jcritt,pcritq,qcritq)
  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. *
  31. *     En ce qui concerne jcritq, le critère est généralement le mini de
  32. *     : la qualité, la taille, l'inverse de la taille
  33. *     Suivant le chiffre des centaines c de jcritq :
  34. C     * c=0 : on renvoie le critère complet (mini des 3)
  35. C     * c=1 : on renvoie la qualité
  36. C     * c=2 : on renvoie la taille
  37. C     * c=3 : on renvoie l'inverse de la taille
  38. C
  39. C
  40. C
  41. C LANGAGE     : ESOPE
  42. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  43. C               mél : gounand@semt2.smts.cea.fr
  44. C***********************************************************************
  45. C VERSION    : v1, 27/11/2017, version initiale
  46. C HISTORIQUE : v1, 27/11/2017, création
  47. C HISTORIQUE : v2, 10/12/2025, on met des critères de qualité
  48. C     similaires à DEDUADAP
  49. C HISTORIQUE :
  50. C***********************************************************************
  51. -INC CCGEOME
  52. -INC PPARAM
  53. -INC CCOPTIO
  54. -INC CCREEL
  55. -INC SMCOORD
  56. -INC TMATOP1
  57. *-INC SMETRIQ
  58.       POINTEUR KCMETR.METRIQ
  59. *-INC SMELEMX
  60. -INC SMLREEL
  61.       PARAMETER(NMET=6)
  62.       DIMENSION XMET(NMET)
  63. *      DIMENSION XMET2(2,2)
  64. *      DIMENSION XMET3(3,3)
  65.       DIMENSION XMETV(3,3)
  66.       DIMENSION XJAC(3,3)
  67.       DIMENSION XJTMJ(3,3)
  68.       DIMENSION XMJ(3,3)
  69.       REAL*8 DXI(3)
  70.       REAL*8 DARET(6)
  71.       DIMENSION A(3,3),D(3)
  72. * Derivative of the affine barycentric map
  73.       DIMENSION DABM(4,3)
  74. * Simplex node coordinates
  75.       DIMENSION SNCO(3,4)
  76.  
  77.  
  78.       DIMENSION IDXSYM(3,3,3)
  79.       LOGICAL LROT
  80.       LOGICAL LQUAL0
  81. *
  82. * Statement functions
  83. *      DISTA(A,B,C,D)=SQRT((A-C)*(A-C)+(B-D)*(B-D))
  84. *      DISTB(A,B,C,D,E,F)=SQRT((A-D)*(A-D)+(B-E)*(B-E)+(C-F)*(C-F))
  85. *      DISTA(A,B)=SQRT(A*A+B*B)
  86. *      DISTB(A,B,C)=SQRT(A*A+B*B+C*C)
  87. *      DETTRI(A11,A12,A21,A22)=A11*A22-A21*A12
  88.       DETTET(A11,A12,A13,A21,A22,A23,A31,A32,A33)=
  89.      &  A11*(A22*A33-A23*A32)
  90.      &       +A12*(A23*A31-A21*A33)
  91.      &       +A13*(A21*A32-A22*A31)
  92. *
  93.       DATA ((XMETV(I,J),I=1,3),J=1,3) /9*0.D0/
  94.       DATA ((XJAC(I,J),I=1,3),J=1,3)  /9*0.D0/
  95.       DATA ((XJTMJ(I,J),I=1,3),J=1,3) /9*0.D0/
  96.       DATA ((XMJ(I,J),I=1,3),J=1,3) /9*0.D0/
  97.       DATA ((A(I,J),I=1,3),J=1,3) /9*0.D0/
  98.       DATA (D(I),I=1,3) /3*0.D0/
  99.       DATA ((DABM(I,J),I=1,4),J=1,3) /12*0.D0/
  100.       DATA ((SNCO(I,J),I=1,3),J=1,4) /12*0.D0/
  101.       DATA (((IDXSYM(I,J,K),I=1,1),J=1,1),K=1,1) /1/
  102.       DATA (((IDXSYM(I,J,K),I=1,2),J=1,2),K=2,2) /1,2,2,3/
  103.       DATA (((IDXSYM(I,J,K),I=1,3),J=1,3),K=3,3) /1,2,4,2,3,5,4,5,6/
  104. *
  105. *
  106. * Executable statements
  107. *
  108. *  Choix du critere de qualite sans metrique
  109. *  0 : Coupez
  110. *  1 : XALIN2
  111. *  2 : XALIN1
  112. *  3 : 2D : D(2)/D(1)
  113. *      3D : D(3)/D(1)
  114. *     write(ioimp,*) 'quali6: NKPVIR=',NKPVIR
  115.       JCRITC=JCRITT/100
  116.       JCRITQ=JCRITT-(JCRITC*100)
  117.       IF (IIMPI.EQ.666) THEN
  118.          WRITE(IOIMP,*) 'JCRITT,JCRITC,JCRITQ=',JCRITT,JCRITC,JCRITQ
  119.          WRITE(IOIMP,*) 'PCRITQ,QCRITQ=',PCRITQ,QCRITQ
  120.       ENDIF
  121. *
  122.       NDQC=0
  123.       IDIMP1=IDIM+1
  124. *      NBNN=NUMX(/1)
  125.       NBNN=NNCOU
  126. *      NBELEM=NUMX(/2)
  127.       IF (NBNN.NE.IDIMP1) THEN
  128.          CALL ERREUR(5)
  129.          RETURN
  130.       ENDIF
  131.       IF
  132.      $     (.NOT.(IELDEB.GE.1.AND.IELFIN.GE.IELDEB.AND.NLCOU.GE.IELFIN
  133.      $     .AND.NUMX(/2).GE.NLCOU))  THEN
  134.          WRITE(IOIMP,*) 'coucou quali6'
  135.          write(ioimp,*) 'IELDEB=',IELDEB
  136.          write(ioimp,*) 'IELFIN=',IELFIN
  137.          write(ioimp,*) 'NLCOU=',NLCOU
  138.          write(ioimp,*) 'NUM2=',NUMX(/2)
  139.          write(ioimp,*) 'MELEMX=',MELEMX
  140.          CALL ECMELX(MELEMX,0)
  141.          CALL ERREUR(5)
  142.          RETURN
  143.       ENDIF
  144. *         XPET=XPETIT*10.D0
  145.       XPET=sqrt(XPETIT)
  146.  
  147. *     DO 10 IBELEM=1,NUMX(/2)
  148.       DO 10 IBELEM=IELDEB,IELFIN
  149.          LQUAL0=.FALSE.
  150. *         WRITE(IOIMP,*) 'IBELEM=',IBELEM
  151. *     Calcul de la métrique moyenne M soit dans XMETD (scalaire)
  152. *                                     soit dans XMETV (tenseur SPD)
  153. * Derivative of the affine barycentric map M : lambda(x) = Mx +c
  154. *     Les coordonnees barycentriques sont definies par rapport au
  155. *     simplex regulier de cote 1, centre sur l'origine. Le noeud sommet
  156. *     a toutes ses coordonnees nulles sauf la derniere
  157. * Initialisations au premier pas
  158.          IF (IBELEM.EQ.IELDEB) THEN
  159.             IF (IDIM.GE.1) THEN
  160.                DABM(1,1)=-1.D0
  161.                DABM(2,1)=+1.D0
  162.                IF (IDIM.GE.2) THEN
  163.                   DABM(1,2)=-1.D0/SQRT(3.D0)
  164.                   DABM(2,2)=-1.D0/SQRT(3.D0)
  165.                   DABM(3,2)=+2.D0/SQRT(3.D0)
  166.                   IF (IDIM.GE.3) THEN
  167.                      DABM(1,3)=-1.D0/SQRT(6.D0)
  168.                      DABM(2,3)=-1.D0/SQRT(6.D0)
  169.                      DABM(3,3)=-1.D0/SQRT(6.D0)
  170.                      DABM(4,3)=+3.D0/SQRT(6.D0)
  171.                   ENDIF
  172.                ENDIF
  173.             ENDIF
  174. *
  175.             IF (IMET.EQ.1) XMETD=1.D0/DENSIT
  176.             IF (IMET.EQ.2) XMETD=1.D0/XDENS
  177.             IF (IMET.EQ.3) NFMET=1
  178.             IF (IMET.EQ.4) NFMET=IDIM*(IDIM+1)/2
  179.          ENDIF
  180. *
  181.          if (imet.gt.0) then
  182.             IF (IMET.GE.1.AND.IMET.LE.3) THEN
  183.                IF (IMET.EQ.3) THEN
  184.                   YDENS=0.D0
  185.                   DO I=1,IDIMP1
  186.                      INO=NUMX(I,IBELEM)
  187.                      IF (NKPVIR.NE.0) THEN
  188.                         IF (INO.LE.NKPVIR) GOTO 10
  189.                      ENDIF
  190. *     Ici on fait la moyenne arithmétique
  191. *     mais kcmetr contient le log du tenseur si imomet=1
  192.                      YDENS=YDENS+KCMETR.XIN(1,INO)
  193.                   ENDDO
  194.                   YDENS=YDENS/IDIMP1
  195.                   if (imomet.eq.1) then
  196.                      YDENS=EXP(YDENS)
  197.                   endif
  198.                   XMETD2=YDENS
  199.                   XMETD=SQRT(XMETD)
  200.                ELSE
  201.                   XMETD2=XMETD**2
  202.                ENDIF
  203.             ELSEIF (IMET.EQ.4) THEN
  204.                DO J=1,NFMET
  205.                   XMET(J)=0.D0
  206.                ENDDO
  207.                DO I=1,IDIMP1
  208.                   INO=NUMX(I,IBELEM)
  209.                   IF (NKPVIR.NE.0) THEN
  210.                      IF (INO.LE.NKPVIR) GOTO 10
  211.                   ENDIF
  212.                   DO J=1,NFMET
  213.                      XMET(J)=XMET(J)+KCMETR.XIN(J,INO)
  214.                   ENDDO
  215.                ENDDO
  216.                DO J=1,NFMET
  217.                   XMET(J)=XMET(J)/IDIMP1
  218.                ENDDO
  219. *
  220.                if (imomet.eq.1) then
  221.                   DO J=1,IDIM
  222.                      DO I=1,IDIM
  223.                         A(I,J)=XMET(IDXSYM(I,J,IDIM))
  224.                      ENDDO
  225.                   ENDDO
  226. *     Exponentielle du tenseur symétrique
  227.                   IOTENS=8
  228.                   IKAS=3
  229.                   CALL TENS2(IOTENS,IKAS,A,D,XMETV)
  230.                   IF (IERR.NE.0) RETURN
  231.                else
  232.                   DO J=1,IDIM
  233.                      DO I=1,IDIM
  234.                         XMETV(I,J)=XMET(IDXSYM(I,J,IDIM))
  235.                      ENDDO
  236.                   ENDDO
  237.                endif
  238.             ELSE
  239.                WRITE(IOIMP,*) 'quali6 imet=',IMET
  240.                CALL ERREUR(5)
  241.                RETURN
  242.             ENDIF
  243.          endif
  244. * Determinant de la metrique
  245.          if (jcritq.eq.0) then
  246.             if (imet.gt.0) then
  247.                IF (IMET.GE.1.AND.IMET.LE.3) THEN
  248.                   XDETMD=SQRT(XMETD2)**IDIM
  249.                ELSEIF (IMET.EQ.4) THEN
  250.                   IF (IDIM.EQ.1) THEN
  251.                      XDETMD=XMETV(1,1)
  252.                   ELSEIF (IDIM.EQ.2) THEN
  253.                      XDETMD=XMETV(1,1)*XMETV(2,2)-XMETV(2,1)*XMETV(1,2)
  254.                   ELSEIF (IDIM.EQ.3) THEN
  255.                      XDETMD=DETTET(XMETV(1,1),XMETV(1,2),XMETV(1,3)
  256.      $                    ,XMETV(2,1),XMETV(2,2),XMETV(2,3),XMETV(3,1)
  257.      $                    ,XMETV(3,2),XMETV(3,3))
  258.                   ELSE
  259.                      WRITE(IOIMP,*) 'quali6 jcrit=0 imet=4 idim=',IDIM
  260.                      INTERR(1)=IDIM
  261.                      CALL ERREUR(709)
  262.                      RETURN
  263.                   ENDIF
  264.                   XDETMD=SQRT(XDETMD)
  265.                ELSE
  266.                   WRITE(IOIMP,*) 'quali6 imet=',IMET
  267.                   CALL ERREUR(5)
  268.                   RETURN
  269.                ENDIF
  270.             endif
  271. * Volume du simplex
  272.             IF (IDIM.EQ.1) THEN
  273.                I0=NUMX(1,IBELEM)
  274.                I1=NUMX(2,IBELEM)
  275.                IF (NKPVIR.NE.0) THEN
  276.                   IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR) goto 10
  277.                ENDIF
  278.                IP0=(I0-1)*IDIMP1
  279.                IP1=(I1-1)*IDIMP1
  280.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  281.                XVOLO=X10
  282.             ELSEIF (IDIM.EQ.2) THEN
  283.                I0=NUMX(1,IBELEM)
  284.                I1=NUMX(2,IBELEM)
  285.                I2=NUMX(3,IBELEM)
  286.                IF (NKPVIR.NE.0) THEN
  287.                   IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR.OR.I2.LE.NKPVIR) goto
  288.      $                 10
  289.                ENDIF
  290.                IP0=(I0-1)*IDIMP1
  291.                IP1=(I1-1)*IDIMP1
  292.                IP2=(I2-1)*IDIMP1
  293.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  294.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  295.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  296.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  297.                XVOLO=(X10*Y20-X20*Y10)/2.D0
  298.             ELSEIF (IDIM.EQ.3) THEN
  299.                I0=NUMX(1,IBELEM)
  300.                I1=NUMX(2,IBELEM)
  301.                I2=NUMX(3,IBELEM)
  302.                I3=NUMX(4,IBELEM)
  303.                IF (NKPVIR.NE.0) THEN
  304.                   IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR.OR.I2.LE.NKPVIR.OR.
  305.      $                 I3.LE.NKPVIR) goto 10
  306.                ENDIF
  307.                IP0=(I0-1)*IDIMP1
  308.                IP1=(I1-1)*IDIMP1
  309.                IP2=(I2-1)*IDIMP1
  310.                IP3=(I3-1)*IDIMP1
  311.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  312.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  313.                Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
  314.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  315.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  316.                Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
  317.                X30=XCOOR(IP3+1)-XCOOR(IP0+1)
  318.                Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
  319.                Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
  320.                XVOLO=(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  321.      $              /6.D0
  322.             ENDIF
  323. * Déterminant de la métrique
  324.             if (imet.gt.0) then
  325.                XVOLO=XVOLO*XDETMD
  326. *               write(ioimp,*) 'XVOLO,XDETMD=',XVOLO,XDETMD
  327.             endif
  328.             IF (IIMPI.EQ.666) THEN
  329.                write(ioimp,*) 'Xvol,xvolo,xvtol=',Xvol,xvolo,xvtol
  330.             ENDIF
  331.             xvol=abs(xvolo)
  332.             if (xvol.LT.xvtol) then
  333.                LQUAL0=.TRUE.
  334.                goto 666
  335.             endif
  336. * Calcul de la longueur de reference XLAARI
  337.             XLAR=0.D0
  338.             IARET=0
  339.             DO IBNN=1,NBNN
  340.                DO JBNN=IBNN+1,NBNN
  341.                   I0=NUMX(IBNN,IBELEM)
  342.                   I1=NUMX(JBNN,IBELEM)
  343.                   IP0=(I0-1)*IDIMP1
  344.                   IP1=(I1-1)*IDIMP1
  345.                   DO IIDIM=1,IDIM
  346.                      DXI(IIDIM)=XCOOR(IP1+IIDIM)-XCOOR(IP0+IIDIM)
  347.                      IF (IMET.GE.1.AND.IMET.LE.3) THEN
  348.                         DXI(IIDIM)=DXI(IIDIM)*XMETD
  349.                      ENDIF
  350.                   ENDDO
  351.                   DXLAR2=0.D0
  352.                   IF (IMET.EQ.4) THEN
  353.                      DO J=1,IDIM
  354.                         DO I=1,IDIM
  355.                            DXLAR2=DXLAR2+XMETV(I,J)*DXI(I)*DXI(J)
  356.                         ENDDO
  357.                      ENDDO
  358.                   ELSE
  359.                      DO I=1,IDIM
  360.                         DXLAR2=DXLAR2+DXI(I)*DXI(I)
  361.                      ENDDO
  362.                   ENDIF
  363.                   DXLAR=SQRT(DXLAR2)
  364.                   IARET=IARET+1
  365.                   DARET(IARET)=DXLAR
  366.                   XLAR=XLAR+DXLAR
  367.                ENDDO
  368.             ENDDO
  369.             NARET=((NBNN-1)*NBNN)/2
  370.             XLAARI=XLAR/NARET
  371. * Coefficient de normalisation
  372.             I=IDIM
  373.             XCOQ=DFACT(I)*(SQRT((DBLE(2**I))/(DBLE(I+1))))
  374.             XQUALC=((XVOL*XCOQ)**(1.D0/IDIM))/XLAARI
  375.          else
  376. * Calcul du jacobien de la transformation geometrique entre
  377. * l'element regulier de coté 1 et l'element courant
  378. *   Coordonnees des noeuds
  379.             DO J=1,IDIMP1
  380.                INOD=NUMX(J,IBELEM)
  381.                IF (NKPVIR.NE.0) THEN
  382.                   IF (INOD.LE.NKPVIR) goto 10
  383.                ENDIF
  384.                IPNOD=(INOD-1)*IDIMP1
  385.                DO I=1,IDIM
  386.                   SNCO(I,J)=XCOOR(IPNOD+I)
  387.                ENDDO
  388.             ENDDO
  389. *         write(ioimp,*) 'SNCO,I,J=',IDIM,IDIMP1
  390. *         write(ioimp,*) ((SNCO(I,J),I=1,IDIM),J=1,IDIMP1)
  391. *         write(ioimp,*) 'DABM,I,J=',IDIMP1,IDIM
  392. *         write(ioimp,*) ((DABM(I,J),I=1,IDIMP1),J=1,IDIM)
  393. *   Matrice Jacobienne de la transformation J = SNCO*DABM
  394.             DO J=1,IDIM
  395.                DO I=1,IDIM
  396.                   XIJ=0.D0
  397.                   DO K=1,IDIMP1
  398.                      XIJ=XIJ+SNCO(I,K)*DABM(K,J)
  399.                   ENDDO
  400.                   XJAC(I,J)=XIJ
  401.                ENDDO
  402.             ENDDO
  403.             IF (IDIM.EQ.1) THEN
  404.                XDETJ=XJAC(1,1)
  405.             ELSEIF (IDIM.EQ.2) THEN
  406.                XDETJ=XJAC(1,1)*XJAC(2,2)-XJAC(2,1)*XJAC(1,2)
  407.             ELSEIF (IDIM.EQ.3) THEN
  408.                XDETJ=DETTET(XJAC(1,1),XJAC(1,2),XJAC(1,3),XJAC(2
  409.      $              ,1),XJAC(2,2),XJAC(2,3),XJAC(3,1),XJAC(3,2)
  410.      $             ,XJAC(3,3))
  411.             ELSE
  412.                INTERR(1)=IDIM
  413.                CALL ERREUR(709)
  414.                RETURN
  415.             ENDIF
  416. *            write(ioimp,*) 'XDETJ=',XDETJ
  417. *         write(ioimp,*) 'XJAC,I,J=',IDIM,IDIM
  418. *         write(ioimp,*) ((XJAC(I,J),I=1,IDIM),J=1,IDIM)
  419. * Matrice JtMJ
  420.             IF (IMET.LT.4) THEN
  421.                DO K=1,IDIM
  422.                   DO I=1,IDIM
  423.                      XIK=0.D0
  424.                      DO J=1,IDIM
  425.                         XIK=XIK+XJAC(J,I)*XJAC(J,K)
  426.                      ENDDO
  427.                      IF (IMET.EQ.0) THEN
  428.                         XJTMJ(I,K)=XIK
  429.                      ELSE
  430.                         XJTMJ(I,K)=XIK*XMETD2
  431.                      ENDIF
  432.                   ENDDO
  433.                ENDDO
  434.                IF (IMET.EQ.0) THEN
  435.                   XDETM=1.D0
  436.                ELSE
  437.                   XDETM=XMETD2**IDIM
  438.                ENDIF
  439.             ELSE
  440.                DO L=1,IDIM
  441.                   DO J=1,IDIM
  442.                      XJL=0.D0
  443.                      DO K=1,IDIM
  444. * Utilisons la symetrie de XMETV
  445.                         XJL=XJL+XMETV(K,J)*XJAC(K,L)
  446.                      ENDDO
  447.                      XMJ(J,L)=XJL
  448.                   ENDDO
  449.                ENDDO
  450.                DO L=1,IDIM
  451.                   DO I=1,IDIM
  452.                      XIL=0.D0
  453.                      DO J=1,IDIM
  454.                         XIL=XIL+XJAC(J,I)*XMJ(J,L)
  455.                      ENDDO
  456.                      XJTMJ(I,L)=XIL
  457.                   ENDDO
  458.                ENDDO
  459.                IF (IDIM.EQ.1) THEN
  460.                   XDETM=XMETV(1,1)
  461.                ELSEIF (IDIM.EQ.2) THEN
  462.                   XDETM=XMETV(1,1)*XMETV(2,2)-XMETV(2,1)*XMETV(1,2)
  463.                ELSEIF (IDIM.EQ.3) THEN
  464.                   XDETM=DETTET(XMETV(1,1),XMETV(1,2),XMETV(1,3),XMETV(2
  465.      $                 ,1),XMETV(2,2),XMETV(2,3),XMETV(3,1),XMETV(3,2)
  466.      $                 ,XMETV(3,3))
  467.                ELSE
  468.                   INTERR(1)=IDIM
  469.                   CALL ERREUR(709)
  470.                   RETURN
  471.                ENDIF
  472.             ENDIF
  473. *      WRITE(IOIMP,*) 'XDETM=',XDETM
  474. *         write(ioimp,*) 'XJTMJ,I,J=',IDIM,IDIM
  475. *     write(ioimp,*) ((XJTMJ(I,J),I=1,IDIM),J=1,IDIM)
  476.             XDETJM=XDETJ*SQRT(XDETM)
  477. *     WRITE(IOIMP,*) 'XDETJM=',XDETJM
  478.             IF (IIMPI.EQ.666) THEN
  479.                write(ioimp,*) 'XDETJ,XDETM,Xdetjm,xvtol=',xdetj,xdetm
  480.      $              ,Xdetjm,xvtol
  481.             ENDIF
  482. *            if (XDETJM.LT.xvtol) then
  483.             if (ABS(XDETJM).LT.xvtol) then
  484.                LQUAL0=.TRUE.
  485.                goto 666
  486.             endif
  487. * Determinant et trace de JTMJ
  488.             IF (IDIM.EQ.1) THEN
  489.                D(1)=XDETJM**2
  490.             ELSEIF (IDIM.EQ.2) THEN
  491. *     XDET=XJTMJ(1,1)*XJTMJ(2,2)-XJTMJ(2,1)*XJTMJ(1,2)
  492.                CALL JACOD2(XJTMJ,D)
  493.             ELSEIF (IDIM.EQ.3) THEN
  494. *     XDET=DETTET(XJTMJ(1,1),XJTMJ(1,2),XJTMJ(1,3),XJTMJ(2
  495. *     $           ,1),XJTMJ(2,2),XJTMJ(2,3),XJTMJ(3,1),XJTMJ(3,2)
  496. *     $           ,XJTMJ(3,3))
  497.                CALL JACOD3(XJTMJ,3,D)
  498.             ELSE
  499.                WRITE(IOIMP,*) 'quali6 idim=',IDIM
  500.                INTERR(1)=IDIM
  501.                CALL ERREUR(709)
  502.                RETURN
  503.             ENDIF
  504. * On stocke les racines des valeurs propres (longueurs) dans DARET
  505. *         WRITE(IOIMP,*) '1',(D(II),II=1,IDIM)
  506.             DO I=1,IDIM
  507. *     D(I)=ABS(D(I))
  508.                DARET(I)=SQRT(MAX(D(I),XZERO))
  509.             ENDDO
  510.             NARET=IDIM
  511. *
  512. * Pour le calcul de XALIN2 et XALIN1, on stocke D(I) ou DARET(I) dans DXI(I)
  513. *
  514. * Calcul de XALIN2 = inverse de QALI dans Huang
  515. *
  516. *     Calcul de XALIN1 pareil que XALIN2 mais avec les valeurs propres
  517. *     au lieu de leur carré
  518.             IF (JCRITQ.EQ.1.OR.JCRITQ.EQ.2) THEN
  519.                IF (IDIM.EQ.1) THEN
  520.                   XALIN=1
  521.                ELSE
  522.                   if (jcritq.eq.1) then
  523.                      DO I=1,IDIM
  524.                         DXI(I)=D(I)
  525.                      ENDDO
  526.                   else
  527.                      DO I=1,IDIM
  528.                         DXI(I)=DARET(I)
  529.                      ENDDO
  530.                   endif
  531. * Trace
  532.                   XTR=DXI(1)
  533.                   DO I=2,IDIM
  534.                      XTR=XTR+DXI(I)
  535.                   ENDDO
  536.                   XLT=XTR/IDIM
  537.                   if (jcritq.eq.1) then
  538.                      IF (IDIM.EQ.2) THEN
  539.                         XLTD=XLT
  540.                      ELSEIF (IDIM.EQ.3) THEN
  541.                         XLTD=XLT*SQRT(XLT)
  542.                      ELSE
  543.                         INTERR(1)=IDIM
  544.                         CALL ERREUR(709)
  545.                         RETURN
  546.                      ENDIF
  547.                   else
  548.                      XLTD=XLT**IDIM
  549.                   endif
  550.                   XALIN=ABS(XDETJM)/(MAX(XLTD,XPET))
  551.                   IF (IDIM.EQ.3) XALIN=SQRT(XALIN)
  552.                ENDIF
  553.             ENDIF
  554.          endif
  555. * Par rapport au livre de Huang p.205, XALIN2 vaut 1 / (Qali^(n-1)) (n>=2)
  556. * Comme Qali est minore par le rapport d'aspect d'un element, il
  557. * faut comparer XALIN2 a des (rapports de longueur)^(n-1)
  558. * Il y a donc un carre en dimension 3 cf. les expressions de XQUALN
  559. * plus bas (voir aussi deadutil.procedur qui exprime des indicateurs
  560. * en rapports de longueur)
  561. * On a choisi d'exprimer XALIN2 en fonction de XDETJM directement
  562. * car la presque nullite de XDETJ permet de detecter les elements
  563. * plats.
  564. * Si on l'eleve a la puissance (1/3), ca ne marche pas.
  565. * On devra sans doute faire mieux pour etre vraiment robuste....
  566. *         write(ioimp,*) 'XALIN2=',XALIN2
  567. * Les valeurs propres sont censees etre positives mais pas garanti
  568. *     donc on prend la valeur absolue et on reclasse
  569. * Raccourci pour les elements de qualite nulle
  570.  666     CONTINUE
  571.          JELDEB=IELDEB+NDQC
  572. *     write(ioimp,*) 'jeldeb,prog=',jeldeb,prog(jeldeb)
  573.          IF (LQUAL0) THEN
  574.             IF (IIMPI.EQ.666) write(ioimp,*) 'Cas LQUAL0 IELEM=',JELDEB
  575.             PROG(JELDEB)=0.D0
  576.          ELSE
  577.             IF (IDIM.EQ.1) THEN
  578.                XQUALN=1.D0
  579.             ELSE
  580.                IF (JCRITQ.EQ.0) THEN
  581.                   XQUALN=XQUALC
  582.                ELSEIF (JCRITQ.EQ.1.OR.JCRITQ.EQ.2) THEN
  583.                   XQUALN=XALIN
  584.                ELSEIF (JCRITQ.EQ.3) THEN
  585.                   IF (DARET(1).NE.XZERO) THEN
  586.                      XQUALN=DARET(IDIM)/DARET(1)
  587.                   ELSE
  588.                      XQUALN=XZERO
  589.                   ENDIF
  590.                ELSE
  591.                   WRITE(IOIMP,*) 'quali6 : jcritq=',JCRITQ
  592.                   CALL ERREUR(5)
  593.                   RETURN
  594.                ENDIF
  595.             ENDIF
  596.             IF (IMET.EQ.0) THEN
  597.                PROG(JELDEB)=XQUALN
  598.             ELSE
  599.                XQUALQ=XQUALN**(1.d0/QCRITQ)
  600. *     WRITE(IOIMP,*) '3',(DARET(II),II=1,NARET)
  601.                IF (JCRITQ.EQ.0) THEN
  602.                   DMIN=XGRAND
  603.                   DMAX=XZERO
  604.                   DO K=1,NARET
  605.                      DMIN=MIN(DMIN,DARET(K))
  606.                      DMAX=MAX(DMAX,DARET(K))
  607.                   ENDDO
  608.                ELSEIF (JCRITQ.GE.1.AND.JCRITQ.LE.3) THEN
  609.                   DMAX=DARET(1)
  610.                   DMIN=DARET(NARET)
  611.                ELSE
  612.                   WRITE(IOIMP,*) 'quali6 : jcritq=',JCRITQ
  613.                   CALL ERREUR(5)
  614.                   RETURN
  615.                ENDIF
  616.                dinf=abs(DMAX)
  617.                if (dinf.lt.xpet) then
  618.                   XL1=XZERO
  619.                else
  620.                   if (pcritq.eq.0.d0) then
  621.                      DMOYP=1.D0
  622.                      DO K=1,NARET
  623.                         DMOYP=DMOYP*DARET(K)
  624.                      ENDDO
  625.                      XP=1.D0/NARET
  626.                      DMOYP=DMOYP**XP
  627.                   elseif (pcritq.gt.100.d0) then
  628.                      dmoyp=abs(dmax)
  629.                   elseif (pcritq.lt.-100.d0) then
  630.                      dmoyp=abs(dmin)
  631.                   else
  632.                      xsumn=(abs(daret(NARET))/dinf)**pcritq
  633.                      DO K=NARET-1,1,-1
  634.                         xsumn=xsumn+(abs(daret(k))/dinf)**pcritq
  635.                      ENDDO
  636.                      xsumn=xsumn/NARET
  637.                      dmoyp=dinf*xsumn**(1.d0/pcritq)
  638.                   endif
  639.                endif
  640.                XL1=DMOYP
  641.                XL1I=1.D0/MAX(XL1,XPET)
  642.                XL1T=MIN(XL1,XL1I)
  643. *               WRITE(IOIMP,'(A,2X,I3,2X,10(2(f9.3),1X))')
  644. *     $          'critq,d1,d2,dmoyp,xl1t,xalin1,xqual=',
  645. *     $           jcritq,pcritq,qcritq,daret(1),daret(2),dmoyp,xl1t
  646. *     ,xalin1,xqual
  647.                IF (JCRITC.EQ.0) THEN
  648.                   PROG(JELDEB)=MIN(XQUALQ,XL1T)
  649.                ELSEIF (JCRITC.EQ.1) THEN
  650.                   PROG(JELDEB)=XQUALQ
  651.                ELSEIF (JCRITC.EQ.2) THEN
  652.                   PROG(JELDEB)=XL1
  653.                ELSEIF (JCRITC.EQ.3) THEN
  654.                   PROG(JELDEB)=XL1I
  655.                ELSE
  656.                   CALL ERREUR(5)
  657.                   RETURN
  658.                ENDIF
  659.             ENDIF
  660.          ENDIF
  661.  
  662. * Scaling ? mmmmm, petit doute
  663.          NDQC=NDQC+1
  664. *         Write(ioimp,*) 'jeldeb,prog=',jeldeb,prog(jeldeb)
  665.  10   CONTINUE
  666.       RETURN
  667. *
  668. * formats
  669. *
  670.   188  FORMAT (2X,12(A6,'=',1PG12.5,2X))
  671. *
  672. * End of subroutine QUALI6
  673. *
  674.       END
  675.  
  676.  

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