Sources Esope | Cast3M

coq3sh
C COQ3SH    SOURCE    CHAT      05/01/12    22:21:34     5004
      SUBROUTINE COQ3SH(QSI,ETA,XEL,VNM)
C=======================================================================
C
C     EVALUATION DES FONCTIONS DE FORME DU DKT OU COQ3
C        BASE SUR UNE INTERPOLATION CUBIQUE DES FCTNS DE FORME
C
C  ENTREES
C   QSI ETA COORDONNEES REDUITES DU POINT DE GAUSS
C   XEL(3,3)=COORDONNEES LOCALES DE L ELEMENT
C         NOEUD 1 2 SUR L AXE QSI  NOEUD 1 EN 0,0,0
C                                  NOEUD 3 DANS LE PLAN QSI ETA
C  SORTIES
C   VNM(3,9)  = FONCTION DE FORME POUR W
C   VNM(1,1:9)= VALEUR   DE W
C   VNM(2,1:9)= DERIVEE     W,X
C   VNM(3,1:9)= DERIVEE     W,Y
C
C   EBERSOLT SEPTEMBRE 85
C                 IL EXISTE UNE AUTRE ROUTINE QUI CALCULE QUE W
C                                         C EST MFDKT
C=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XEL(3,*),VNM(3,*),VNX(9),VNY(9)
      DATA UN,DEUX,TROIS,SIX,UNDEMI/1.D0,2.D0,3.D0,6.D0,.5D0/
      DATA XZER/0.D0/
C
C     MISE A ZERO  DE VNM
C
      DO 10 IA=1,9
      DO 10 IB=1,3
      VNM(IB,IA)=XZER
  10  CONTINUE
C
C         CALCUL DES FONCTIONS D INTERPOLATION
C
      AL=UN-QSI-ETA
      A=QSI*ETA*AL
C
      VNM(1,1)= AL*AL*(TROIS-DEUX*AL)     +A*DEUX
      VNM(1,2)= AL*AL*QSI                 +A*UNDEMI
      VNM(1,3)= AL*AL*ETA                 +A*UNDEMI
C
      VNM(1,4)=(TROIS-DEUX*QSI)*QSI*QSI   +A*DEUX
      VNM(1,5)=(QSI-UN)*QSI*QSI           -A
      VNM(1,6)= QSI*QSI*ETA               +A*UNDEMI
C
      VNM(1,7)=(TROIS-DEUX*ETA)*ETA*ETA   +A*DEUX
      VNM(1,8)= ETA*ETA*QSI               +A*UNDEMI
      VNM(1,9)=(ETA-UN)*ETA*ETA           -A
C
C     DERIVEES PAR RAPPORT A QSI
C
      B = ETA*(AL-QSI)
C
      VNX(1)= SIX*AL*(AL-UN)       +B*DEUX
      VNX(2)=(AL-DEUX*QSI)*AL      +B*UNDEMI
      VNX(3)=-DEUX*AL*ETA          +B*UNDEMI
C
      VNX(4)= SIX*QSI*(UN-QSI)     +B*DEUX
      VNX(5)= QSI*(TROIS*QSI-DEUX) -B
      VNX(6)= DEUX*QSI*ETA         +B*UNDEMI
C
      VNX(7)=                       B*DEUX
      VNX(8)= ETA*ETA              +B*UNDEMI
      VNX(9)=                      -B
C
C     DERIVEES PAR RAPPORT A ETA
C
      C = QSI*(AL-ETA)
C
      VNY(1)= SIX*AL*(AL-UN)        +C*DEUX
      VNY(2)=-DEUX*QSI*AL           +C*UNDEMI
      VNY(3)=(AL-DEUX*ETA)*AL       +C*UNDEMI
C
      VNY(4)=                        C*DEUX
      VNY(5)=                       -C
      VNY(6)= QSI*QSI               +C*DEUX
C
      VNY(7)= SIX*ETA*(UN-ETA)      +C*DEUX
      VNY(8)= DEUX*QSI*ETA          +C*UNDEMI
      VNY(9)= ETA*(TROIS*ETA-DEUX)  -C
C
C         VARIABLES DANS LE REPERE X,Y
C
      X21=XEL(1,2)-XEL(1,1)
      X13=XEL(1,1)-XEL(1,3)
      Y13=XEL(2,1)-XEL(2,3)
C
C    CHANGEMENT DE REPERE POUR W
C
      CC    =VNM(1,2)*X21-VNM(1,3)*X13
      VNM(1,2)=-VNM(1,3)*Y13
      VNM(1,3)=-CC
      CC    =VNM(1,5)*X21-VNM(1,6)*X13
      VNM(1,5)=-VNM(1,6)*Y13
      VNM(1,6)=-CC
      CC    =VNM(1,8)*X21-VNM(1,9)*X13
      VNM(1,8)=-VNM(1,9)*Y13
      VNM(1,9)=-CC
C
C     CHANGEMENT DE REPERE POUR W,X
C
      VNM(2,1)= VNX(1)/X21
      VNM(2,2)=-VNX(3)*Y13/X21
      VNM(2,3)=-VNX(2)         +VNX(3)*X13/X21
C
      VNM(2,4)= VNX(4)/X21
      VNM(2,5)=-VNX(6)*Y13/X21
      VNM(2,6)=-VNX(5)         +VNX(6)*X13/X21
C
      VNM(2,7)= VNX(7)/X21
      VNM(2,8)=-VNX(9)*Y13/X21
      VNM(2,9)=-VNX(8)         +VNX(9)*X13/X21
C
C     CHANGEMENT DE REPERE POUR W,Y
C
      VNM(3,1)=-VNX(1)*X13/(X21*Y13)  - VNY(1)/Y13
      VNM(3,2)= VNX(3)*X13/X21        + VNY(3)
      VNM(3,3)= VNX(2)*X13/Y13 -VNX(3)*X13*X13/(X21*Y13)
     1         +VNY(2)*X21/Y13 -VNY(3)*X13/Y13
C
      VNM(3,4)=-VNX(4)*X13/(X21*Y13)  - VNY(4)/Y13
      VNM(3,5)= VNX(6)*X13/X21        + VNY(6)
      VNM(3,6)= VNX(5)*X13/Y13 -VNX(6)*X13*X13/(X21*Y13)
     1         +VNY(5)*X21/Y13 -VNY(6)*X13/Y13
C
      VNM(3,7)=-VNX(7)*X13/(X21*Y13)  - VNY(7)/Y13
      VNM(3,8)= VNX(9)*X13/X21        + VNY(9)
      VNM(3,9)= VNX(8)*X13/Y13 -VNX(9)*X13*X13/(X21*Y13)
     1         +VNY(8)*X21/Y13 -VNY(9)*X13/Y13
C
      RETURN
      END