Sources Esope | Cast3M

pb443
C PB443     SOURCE    CHAT      05/01/13    02:10:40     5004
      SUBROUTINE PB443(X,Y,Z,PG,FN,GR,FM,GM,ND,NP,MP,NGG,NPG,NOM2)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C************************************************************************
C
C     CALCULE LES FONCTIONS DE FORME D'UN : Iso-Q2 (iso P1/P0 nc) TE10
C
C
C
C    ^ zeta
C    |
C    |10
C    |
C    |
C    |  9
C    |  /\       /^ eta
C    | /  \     /
C    |/____\   /
C    |7     8 /
C    |       /
C    |      5
C    |     / \
C    |    /   \
C    |   6_____\4
C    |  / \    /\
C    | /   \  /  \
C    |/_____\/____\ ____________________>ksi
C    1      2      3
C
C
C************************************************************************
      CHARACTER*4 NOM2
      REAL*8 AL,BE
      PARAMETER (AL=.5854101966249684D0)
      PARAMETER (BE=.1381966011250105D0)
      REAL*8 X(NPG),Y(NPG),Z(NPG)
      REAL*8 FN(NP,NPG),GR(ND,NP,NPG),PG(NPG)
 
      REAL*8 FM(MP,NPG),GM(ND,MP,NPG)
      INTEGER CONEK(4,8),L,K,I,N,N1,N2,N3,N4,IGAU,IFN
      INTEGER ND,MP,NGG,NP,NPG
      REAL*8  KXYZ(3,10),AA(4)
      REAL*8 T1X,T1Y,T1Z,C1,P1C
      REAL*8 T2X,T2Y,T2Z,C2,P2C
      REAL*8 T3X,T3Y,T3Z,C3,P3C
      REAL*8 T4X,T4Y,T4Z,C4,P4C
C
C Donnees de la numerotation locale des 4 noeuds de chaque tetraedre
C
      DATA CONEK/1,2,6,7,
     &           2,8,6,7,
     &           2,3,4,8,
     &           2,4,6,8,
     &           4,5,6,9,
     &           8,9,4,6,
     &           6,7,8,9,
     &           7,8,9,10/
C
C Donnees des coordonnees des points du tetraedre TE10
C                X         Y         Z
      DATA KXYZ/0.0,      0.0,      0.0,
     &          0.5,      0.0,      0.0,
     &          1.0,      0.0,      0.0,
     &          0.5,      0.5,      0.0,
     &          0.0,      1.0,      0.0,
     &          0.0,      0.5,      0.0,
     &          0.0,      0.0,      0.5,
     &          0.5,      0.0,      0.5,
     &          0.0,      0.5,      0.5,
     &          0.0,      0.0,      1.0/
 
C
C On initialise les fonctions tests et les gradients a 0
C
      DO 10 L=1,NPG
C Poids de Gauss
         PG(L) = 0.1875D0
        DO 22 K=1,MP
        FM(K,L)=0.0D0
        DO 22 I=1,ND
        GM(I,K,L)=0.0D0
 22     CONTINUE
 
        DO 20 K=1,NP
        FN(K,L)=0.0D0
        DO 20 I=1,ND
        GR(I,K,L)=0.0D0
 20     CONTINUE
 
 10   CONTINUE
C
C On traite chaque tetraedre elementaire
C
      DO 40 K=1,8
C
C On travaille sur chacun des tetraedres composant le macro
C
      N1=CONEK(1,K)
      N2=CONEK(2,K)
      N3=CONEK(3,K)
      N4=CONEK(4,K)
C Calcul des equations de plans oppose a chaque point
C
C Plan 1
      T1X=(KXYZ(2,N3)-KXYZ(2,N2))*(KXYZ(3,N4)-KXYZ(3,N2))
     &   -(KXYZ(3,N3)-KXYZ(3,N2))*(KXYZ(2,N4)-KXYZ(2,N2))
      T1Y=(KXYZ(3,N3)-KXYZ(3,N2))*(KXYZ(1,N4)-KXYZ(1,N2))
     &   -(KXYZ(1,N3)-KXYZ(1,N2))*(KXYZ(3,N4)-KXYZ(3,N2))
      T1Z=(KXYZ(1,N3)-KXYZ(1,N2))*(KXYZ(2,N4)-KXYZ(2,N2))
     &   -(KXYZ(2,N3)-KXYZ(2,N2))*(KXYZ(1,N4)-KXYZ(1,N2))
      C1=-(T1X*KXYZ(1,N2)+T1Y*KXYZ(2,N2)+T1Z*KXYZ(3,N2))
      P1C=T1X*KXYZ(1,N1)+T1Y*KXYZ(2,N1)+T1Z*KXYZ(3,N1)+C1
      T1X=T1X/P1C
      T1Y=T1Y/P1C
      T1Z=T1Z/P1C
      C1=C1/P1C
C Plan 2
      T2X=(KXYZ(2,N4)-KXYZ(2,N3))*(KXYZ(3,N1)-KXYZ(3,N3))
     &   -(KXYZ(3,N4)-KXYZ(3,N3))*(KXYZ(2,N1)-KXYZ(2,N3))
      T2Y=(KXYZ(3,N4)-KXYZ(3,N3))*(KXYZ(1,N1)-KXYZ(1,N3))
     &   -(KXYZ(1,N4)-KXYZ(1,N3))*(KXYZ(3,N1)-KXYZ(3,N3))
      T2Z=(KXYZ(1,N4)-KXYZ(1,N3))*(KXYZ(2,N1)-KXYZ(2,N3))
     &   -(KXYZ(2,N4)-KXYZ(2,N3))*(KXYZ(1,N1)-KXYZ(1,N3))
      C2=-(T2X*KXYZ(1,N3)+T2Y*KXYZ(2,N3)+T2Z*KXYZ(3,N3))
      P2C=T2X*KXYZ(1,N2)+T2Y*KXYZ(2,N2)+T2Z*KXYZ(3,N2)+C2
      T2X=T2X/P2C
      T2Y=T2Y/P2C
      T2Z=T2Z/P2C
      C2=C2/P2C
C Plan 3
      T3X=(KXYZ(2,N1)-KXYZ(2,N4))*(KXYZ(3,N2)-KXYZ(3,N4))
     &   -(KXYZ(3,N1)-KXYZ(3,N4))*(KXYZ(2,N2)-KXYZ(2,N4))
      T3Y=(KXYZ(3,N1)-KXYZ(3,N4))*(KXYZ(1,N2)-KXYZ(1,N4))
     &   -(KXYZ(1,N1)-KXYZ(1,N4))*(KXYZ(3,N2)-KXYZ(3,N4))
      T3Z=(KXYZ(1,N1)-KXYZ(1,N4))*(KXYZ(2,N2)-KXYZ(2,N4))
     &   -(KXYZ(2,N1)-KXYZ(2,N4))*(KXYZ(1,N2)-KXYZ(1,N4))
      C3=-(T3X*KXYZ(1,N4)+T3Y*KXYZ(2,N4)+T3Z*KXYZ(3,N4))
      P3C=T3X*KXYZ(1,N3)+T3Y*KXYZ(2,N3)+T3Z*KXYZ(3,N3)+C3
      T3X=T3X/P3C
      T3Y=T3Y/P3C
      T3Z=T3Z/P3C
      C3=C3/P3C
C Plan 4
      T4X=(KXYZ(2,N2)-KXYZ(2,N1))*(KXYZ(3,N3)-KXYZ(3,N1))
     &   -(KXYZ(3,N2)-KXYZ(3,N1))*(KXYZ(2,N3)-KXYZ(2,N1))
      T4Y=(KXYZ(3,N2)-KXYZ(3,N1))*(KXYZ(1,N3)-KXYZ(1,N1))
     &   -(KXYZ(1,N2)-KXYZ(1,N1))*(KXYZ(3,N3)-KXYZ(3,N1))
      T4Z=(KXYZ(1,N2)-KXYZ(1,N1))*(KXYZ(2,N3)-KXYZ(2,N1))
     &   -(KXYZ(2,N2)-KXYZ(2,N1))*(KXYZ(1,N3)-KXYZ(1,N1))
      C4=-(T4X*KXYZ(1,N1)+T4Y*KXYZ(2,N1)+T4Z*KXYZ(3,N1))
      P4C=T4X*KXYZ(1,N4)+T4Y*KXYZ(2,N4)+T4Z*KXYZ(3,N4)+C4
      T4X=T4X/P4C
      T4Y=T4Y/P4C
      T4Z=T4Z/P4C
      C4=C4/P4C
C
C Boucle sur les points de Gauss
C
      DO 50 N=1,(NPG/8)
         IF (NPG.EQ.32) THEN
         DO 60 L=1,4
            AA(L)=BE
 60      CONTINUE
         AA(N)=AL
         ELSE
         DO 70 L=1,4
            AA(L)=0.5D1
 70      CONTINUE
         ENDIF
C
C Indice globale du point de Gauss
C
         IGAU=(K-1)*(NPG/8)+N
C
C Coordonnees du point de Gauss (barycentre des 4 points du tetraedre)
C
         X(IGAU)=AA(1)*KXYZ(1,N1)+AA(2)*KXYZ(1,N2)
     &          +AA(3)*KXYZ(1,N3)+AA(4)*KXYZ(1,N4)
         Y(IGAU)=AA(1)*KXYZ(2,N1)+AA(2)*KXYZ(2,N2)
     &          +AA(3)*KXYZ(2,N3)+AA(4)*KXYZ(2,N4)
         Z(IGAU)=AA(1)*KXYZ(3,N1)+AA(2)*KXYZ(3,N2)
     &          +AA(3)*KXYZ(3,N3)+AA(4)*KXYZ(3,N4)
C
C Boucle sur les fonctions tests
C
C
C Numerotation globale de la fonction test
C
         IFN=(K-1)*4
C
C Calcul des fonctions tests et du gradient au point IGAU
C
         FN((IFN+1),IGAU)=T1X*X(IGAU)+T1Y*Y(IGAU)
     &                   +T1Z*Z(IGAU)+C1
         GR(1,(IFN+1),IGAU)=T1X
         GR(2,(IFN+1),IGAU)=T1Y
         GR(3,(IFN+1),IGAU)=T1Z
         FN((IFN+2),IGAU)=T2X*X(IGAU)+T2Y*Y(IGAU)
     &                   +T2Z*Z(IGAU)+C2
         GR(1,(IFN+2),IGAU)=T2X
         GR(2,(IFN+2),IGAU)=T2Y
         GR(3,(IFN+2),IGAU)=T2Z
         FN((IFN+3),IGAU)=T3X*X(IGAU)+T3Y*Y(IGAU)
     &                   +T3Z*Z(IGAU)+C3
         GR(1,(IFN+3),IGAU)=T3X
         GR(2,(IFN+3),IGAU)=T3Y
         GR(3,(IFN+3),IGAU)=T3Z
         FN((IFN+4),IGAU)=T4X*X(IGAU)+T4Y*Y(IGAU)
     &                   +T4Z*Z(IGAU)+C4
         GR(1,(IFN+4),IGAU)=T4X
         GR(2,(IFN+4),IGAU)=T4Y
         GR(3,(IFN+4),IGAU)=T4Z
C
C Calcul des fonctions tests pour la pression
C
         FM(((K+1)/2),IGAU)=1.0D0
 50   CONTINUE
 40   CONTINUE
 
      IF(NOM2.EQ.'MCF1')THEN
      CO=6.D0 ** (1.D0/3.D0)
      UNSCO=1.D0/CO
      XXXX=-1.D0*UNSCO
      DO 51 L=1,NPG
      FM(1,L)=1.D0-(X(L)+Y(L)+Z(L))*UNSCO
      FM(2,L)=X(L)*UNSCO
      FM(3,L)=Y(L)*UNSCO
      FM(4,L)=Z(L)*UNSCO
C
         GM(1,1,L)=XXXX
         GM(2,1,L)=XXXX
         GM(3,1,L)=XXXX
C
         GM(1,2,L)=UNSCO
         GM(2,2,L)=0.D0
         GM(3,2,L)=0.D0
C
         GM(1,3,L)=0.D0
         GM(2,3,L)=UNSCO
         GM(3,3,L)=0.D0
C
         GM(1,4,L)=0.D0
         GM(2,4,L)=0.D0
         GM(3,4,L)=UNSCO
C
 51   CONTINUE
      ENDIF
 
C     write(6,*)'X '
C     write(6,1008)X
 1008 FORMAT(8(1X,1PE11.4))
 
      RETURN
      END