Sources Esope | Cast3M

cq4loc
C CQ4LOC    SOURCE    CB215821  26/08/24    21:15:50     12622          
      SUBROUTINE CQ4LOC (XE,  XEL,BPSS,NOQUAL,  IVRF)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
************************************************************************
*
*                             C Q 4 L O C
*                             -----------
*
* FONCTION:
* ---------
*
*     Calcul de caract{ristiques d'un {l{ment COQ4.
*
* PARAMETRES:   (E)=ENTREE   (S)=SORTIE   (+ = CONTENU DANS UN COMMUN)
* -----------
*
*     XE      (E)  Coordonn{es des 4 noeuds.
*     IVRF    (E)  = 1  si demande de v{rification de l'{l{ment,
*                  = 0  sinon.
*     XEL     (S)  Coordonn{es locales des 4 noeuds.
*     BPSS    (S)  Matrice de passage.
*     NOQUAL  (S)  Indice de non-qualit{:
*                  = 0  si OK,
*                  = 1  si noeuds trop voisins,
*                  = 3  SI NOEUDS NON COPLANAIRES.
*                  (fourni si demande de v{rification d'{l{ment)
*
      REAL*8 XE(3,4),XEL(3,4),BPSS(3,3)
      INTEGER NOQUAL,IVRF
*
* CONSTANTES:
* -----------
*
*     IND4   = indi\age circulaire de 1 @ 4.
*
      INTEGER IND4(0:5)
*
* VARIABLES:
* ----------
*
*     QSI1   = vecteur norm{ de la m{diane allant de 4-1 vers 2-3.
*     ETA1   = vecteur norm{ de la m{diane allant de 1-2 vers 3-4.
*     X1, Y1 = vecteurs du rep}re local, dans le plan moyen de
*              l'{l{ment.
*     Z1     = vecteur du rep}re local, normal au plan moyen de
*              l'{l{ment.
*
      REAL*8 QSI1(3),ETA1(3),X1(3),Y1(3),Z1(3)
      REAL*8 XD(3,4),U1(3),V1(3)
*
* MODE DE FONCTIONNEMENT:
* -----------------------
*
*     Pour le calcul du rep}re local et de la matrice de passage, on
*     fait une estimation du plan moyen.
*
* AUTEUR, DATE DE CREATION:
* -------------------------
*
*     PASCAL MANIGOT     09 JUILLET 1991
*
* LANGAGE:
* --------
*
*     FORTRAN77
*
************************************************************************
*
      DATA IND4/4,1,2,3,4,1/
*
      NOQUAL=0
C
C     VERIFICA DISTANZA MINIMA DEI PUNTI DELL ELEMENTO
C                                CALIBRE PAR RAPPORT AU PERIMETRE
*+*   A virer ?
      IF (IVRF.NE.1) GO TO 6
      PP=0.D0
      DO 2 I=1,4
      II=I+1
      IF(II.EQ.5) II=1
      C1 =  ABS(XE(1,I)-XE(1,II))
      C2 =  ABS(XE(2,I)-XE(2,II))
      C3 =  ABS(XE(3,I)-XE(3,II))
      C1 = C1*C1+C2*C2+C3*C3
      PP = PP + SQRT(C1)
   2  CONTINUE
      DMIN=PP/50.D0
      DO 103 I=1,3
        I1=I+1
        DO 3 N=I1,4
          IF(ABS(XE(1,I)-XE(1,N)).LE.DMIN.AND.
     $       ABS(XE(2,I)-XE(2,N)).LE.DMIN.AND.
     $       ABS(XE(3,I)-XE(3,N)).LE.DMIN)    THEN
          NOQUAL=1
          RETURN
          ENDIF
  3   CONTINUE
  103 CONTINUE
  6   CONTINUE
*
*     Calcul du rep}re local
*     ----------------------
*
*              Y
*         4    |          3
*          *---|---------*
*          |   |         |
*          |   |         |
*          |   |         |
*          |   +------------X
*          |             |
*          *-------------*
*         1               2
*
*
*     Calcul des m{dianes:
      QSI1(1) = XE(1,2)+XE(1,3) - XE(1,1)-XE(1,4)
      QSI1(2) = XE(2,2)+XE(2,3) - XE(2,1)-XE(2,4)
      QSI1(3) = XE(3,2)+XE(3,3) - XE(3,1)-XE(3,4)
      CALL NORMER (QSI1)
      ETA1(1) = XE(1,3)+XE(1,4) - XE(1,1)-XE(1,2)
      ETA1(2) = XE(2,3)+XE(2,4) - XE(2,1)-XE(2,2)
      ETA1(3) = XE(3,3)+XE(3,4) - XE(3,1)-XE(3,2)
      CALL NORMER (ETA1)
*
*     Normale = Normale aux 2 m{dianes.
      Z1(1)= QSI1(2)*ETA1(3) - QSI1(3)*ETA1(2)
      Z1(2)= QSI1(3)*ETA1(1) - QSI1(1)*ETA1(3)
      Z1(3)= QSI1(1)*ETA1(2) - QSI1(2)*ETA1(1)
      CALL NORMER (Z1)
*
*     Axes dans le Plan = bissectrices des bissectrices des m{dianes
*                       = m{dianes pour un {l{ment rectangulaire
      U1(1) = QSI1(1) - ETA1(1)
      U1(2) = QSI1(2) - ETA1(2)
      U1(3) = QSI1(3) - ETA1(3)
      CALL NORMER (U1)
      V1(1) = QSI1(1) + ETA1(1)
      V1(2) = QSI1(2) + ETA1(2)
      V1(3) = QSI1(3) + ETA1(3)
      CALL NORMER (V1)
*
      X1(1) = U1(1) + V1(1)
      X1(2) = U1(2) + V1(2)
      X1(3) = U1(3) + V1(3)
      CALL NORMER (X1)
*
      Y1(1)=X1(3)*Z1(2)-X1(2)*Z1(3)
      Y1(2)=X1(1)*Z1(3)-X1(3)*Z1(1)
      Y1(3)=X1(2)*Z1(1)-X1(1)*Z1(2)
      CALL NORMER (Y1)
*
*     Coordonn{es locales
*     -------------------
*
      DO 104 J=1,4
      DO 5 I=1,3
        XD(I,J)=XE(I,J)-XE(I,1)
   5  CONTINUE
  104 CONTINUE
*
      DO 10 J=1,4
        XEL(1,J) = XD(1,J)*X1(1)+XD(2,J)*X1(2)+XD(3,J)*X1(3)
        XEL(2,J) = XD(1,J)*Y1(1)+XD(2,J)*Y1(2)+XD(3,J)*Y1(3)
        XEL(3,J) = 0.D0
  10  CONTINUE
*
*     Matrice de passage
*     ------------------
*
      DO 15 I=1,3
        BPSS(1,I)=X1(I)
        BPSS(2,I)=Y1(I)
        BPSS(3,I)=Z1(I)
  15  CONTINUE
*
      IF(IVRF.NE.1) RETURN
*
*     Test de plan{it{
*     ----------------
*
*     Calcul des 4 "normales" locales:
      DO 102 K=1,4
         KP1 = IND4(K+1)
         KM1 = IND4(K-1)
         DO 100 J=1,3
            U1(J) = XE(J,KP1) - XE(J,K)
            V1(J) = XE(J,KM1) - XE(J,K)
  100    CONTINUE
         XD(1,K) = U1(2)*V1(3) - U1(3)*V1(2)
         XD(2,K) = U1(3)*V1(1) - U1(1)*V1(3)
         XD(3,K) = U1(1)*V1(2) - U1(2)*V1(1)
         XXD = (XD(1,K)**2) + (XD(2,K)**2) + (XD(3,K)**2)
         XXD = SQRT(XXD)
         DO 101 J=1,3
            XD(J,K) = XD(J,K) / XXD
  101    CONTINUE
  102 CONTINUE
*
*     Calcul de la non-plan{it{:
      COS13 = XD(1,3)*XD(1,1) + XD(2,3)*XD(2,1) + XD(3,3)*XD(3,1)
      COS24 = XD(1,4)*XD(1,2) + XD(2,4)*XD(2,2) + XD(3,4)*XD(3,2)
      IF (MIN(COS13,COS24) .LT. 0.99999) THEN
*        Non-plan{it{ de 0.25 degr{ ou plus:
         NOQUAL = 3
*        Rq: 0.9999 , qui {quivaut @ 1 degr{, est insuffisant.
      END IF
*
      END