Sources Esope | Cast3M

ella22
C ELLA22    SOURCE    CB215821  16/04/21    21:16:34     8920
      SUBROUTINE ELLA22(CORRES,FLAG,COOR,GAMA,ZA1,MASS,RMAS,
     &                  ZS,NP28,NMAS)
C
      IMPLICIT INTEGER(I-N)
      COMPLEX*16 ZA1,ZS
      REAL*8           COOR,GAMA,RMAS,
     &                 Z11,Z21,Z12,Z31,Z22,Z13,Z32,Z23,Z33,
     &                 X1,Y1,Z1,X2,Y2,Z2,XI1,XI2,XI3,
     &                 XJ1,XJ2,XJ3,XK1,XK2,XK3,
     &                 XL,GG,GX,GY,GZ,DET,DELTA
      INTEGER CORRES,FLAG,I,INP,IZA1,JZA1,MASS,NMAS,NP28,N1,N2
C
C
C   OPERATEUR ELFE LAPLACE ACOU
C
C   PRISE EN COMPTE DE LA MASSE PONCTUELLE
C
C   PARAMETRES :
C
C   COOR   : TABLEAU DES COORDONNEES DES NOEUDS (3,2*NP)
C   GAMA   : TABLEAU DES VECTEUR DEFINISSANT LE PLAN LOCALE OXY (3,2*NP)
C   ZA1    : MATRICE DE RIGIDITE
C   MASS   : LIEU DES MASSES PONCTUELLES
C   RASS   : VALEUR DES MASSES PONCTUELLES
C   ZS     : VALEUR DE LA VARIABLE DE LAPLACE
C   NP     : NOMBRE DE POUTRES
C   NP28   : --> NP*28
C   NMAS   : NOMBRE DE MASSES PONCTUELLES
C
C   AUTEURS : SAINT-DIZIER ET GORCY
C   DATE    : 04 FEVRIER 1991
C
      DIMENSION COOR(3,*),GAMA(3,*)
      DIMENSION ZA1(NP28,*),MASS(4,*)
      DIMENSION RMAS(4,*)
      DIMENSION FLAG(*),CORRES(*)
C
        IF (NMAS.GT.0) THEN
          DO 10 I = 1 , NMAS
            INP = MASS(2,I)
            N1  = 2*INP-1
            N2  = 2*INP
            X1  = COOR (1,N1)
            Y1  = COOR (2,N1)
            Z1  = COOR (3,N1)
            X2  = COOR (1,N2)
            Y2  = COOR (2,N2)
            Z2  = COOR (3,N2)
C
            XL = SQRT((X2-X1)**2 + (Y2-Y1)**2 + (Z2-Z1)**2)
C
C -- VECTEUR UNITAIRE OX REPERE LOCALE
C
            XI1 = (X2-X1)/XL
            XI2 = (Y2-Y1)/XL
            XI3 = (Z2-Z1)/XL
C
C -- VECTEUR UNITAIRE OY REPERE LOCALE
C
            GX = GAMA(1,INP)
            GY = GAMA(2,INP)
            GZ = GAMA(3,INP)
            GG = SQRT(GX*GX + GY*GY + GZ*GZ)
            GX = GX/GG
            GY = GY/GG
            GZ = GZ/GG
C
        DELTA =  SQRT (1.D0 - (XI1*GX + XI2*GY + XI3*GZ)**2)
C
        DET=-(GX*XI3-GZ*XI1)**2-(GY*XI1-GX*XI2)**2-(GY*XI3-GZ*XI2)**2
C
        IF (ABS(DET).LT.1.D-12) THEN
           XJ1 = -XI2
           XJ2 = 0.
           XJ3 = 0.
        ELSE
           XJ1 = (XI2*(GY*XI1-GX*XI2)-XI3*(GX*XI3-GZ*XI1))*DELTA/DET
           XJ2 = (XI3*(GZ*XI2-GY*XI3)-XI1*(GY*XI1-GX*XI2))*DELTA/DET
           XJ3 = (XI1*(GX*XI3-GZ*XI1)-XI2*(GZ*XI2-GY*XI3))*DELTA/DET
        ENDIF
C
C -- VECTEUR UNITAIRE OZ REPERE LOCALE
C
        XK1 = XI2*XJ3 - XI3*XJ2
        XK2 = XI3*XJ1 - XI1*XJ3
        XK3 = XI1*XJ2 - XI2*XJ1
C
        Z11 = XJ2*XK3 - XJ3*XK2
        Z12 = XJ3*XK1 - XJ1*XK3
        Z13 = XJ1*XK2 - XJ2*XK1
        Z21 = XI3*XK2 - XI2*XK3
        Z22 = XI1*XK3 - XI3*XK1
        Z23 = XI2*XK1 - XI1*XK2
        Z31 = XI2*XJ3 - XI3*XJ2
        Z32 = XI3*XJ1 - XI1*XJ3
        Z33 = XI1*XJ2 - XI2*XJ1
C
        IZA1 = MASS(4,I)
        JZA1 = MASS(3,I)
C
           ZA1(IZA1  ,JZA1)   = CMPLX(-1.D0)*RMAS(1,I)*ZS*ZS
           ZA1(IZA1+1,JZA1+1) = CMPLX(-1.D0)*RMAS(1,I)*ZS*ZS
           ZA1(IZA1+2,JZA1+2) = CMPLX(-1.D0)*RMAS(1,I)*ZS*ZS
C
           ZA1(IZA1+3,JZA1+3) = CMPLX(-1.D0)*ZS*ZS*( RMAS(2,I)*Z11*XI1
     *                        + RMAS(3,I)*Z21*XJ1 + RMAS(4,I)*Z31*XK1 )
           ZA1(IZA1+3,JZA1+4) = CMPLX(-1.D0)*ZS*ZS*( RMAS(2,I)*Z12*XI1
     *                        + RMAS(3,I)*Z22*XJ1 + RMAS(4,I)*Z32*XK1 )
           ZA1(IZA1+3,JZA1+5) = CMPLX(-1.D0)*ZS*ZS*( RMAS(2,I)*Z13*XI1
     *                        + RMAS(3,I)*Z23*XJ1 + RMAS(4,I)*Z33*XK1 )
C
           ZA1(IZA1+4,JZA1+3) = CMPLX(-1.D0)*ZS*ZS*( RMAS(2,I)*Z11*XI2
     *                        + RMAS(3,I)*Z21*XJ2 + RMAS(4,I)*Z31*XK2 )
           ZA1(IZA1+4,JZA1+4) = CMPLX(-1.D0)*ZS*ZS*( RMAS(2,I)*Z12*XI2
     *                        + RMAS(3,I)*Z22*XJ2 + RMAS(4,I)*Z32*XK2 )
           ZA1(IZA1+4,JZA1+5) = CMPLX(-1.D0)*ZS*ZS*( RMAS(2,I)*Z13*XI2
     *                        + RMAS(3,I)*Z23*XJ2 + RMAS(4,I)*Z33*XK2 )
C
           ZA1(IZA1+5,JZA1+3) = CMPLX(-1.D0)*ZS*ZS*( RMAS(2,I)*Z11*XI3
     *                        + RMAS(3,I)*Z21*XJ3 + RMAS(4,I)*Z31*XK3 )
           ZA1(IZA1+5,JZA1+4) = CMPLX(-1.D0)*ZS*ZS*( RMAS(2,I)*Z12*XI3
     *                        + RMAS(3,I)*Z22*XJ3 + RMAS(4,I)*Z32*XK3 )
           ZA1(IZA1+5,JZA1+5) = CMPLX(-1.D0)*ZS*ZS*( RMAS(2,I)*Z13*XI3
     *                        + RMAS(3,I)*Z23*XJ3 + RMAS(4,I)*Z33*XK3 )
C
C
10    CONTINUE
C
      END IF
C
      END