Sources Esope | Cast3M

pb503
C PB503     SOURCE    MAGN      10/05/19    21:15:11     6676
      SUBROUTINE PB503(XREF,X,Y,Z,PG,FN,GR,FM,GM,ND,NP,MP,NPG,NOM2)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C*********************************************************************
C                           PYRAMIDE
C   CALCUL DES FONCTIONS DE FORME ET DIVERS POINTS D'INTEGRATION
C                VALEURS F(P)
C                          Z      Y                            XREF
C                          |     /                     P1 | 1.  0.  0.
C                          1.   1.                        |
C                          |   /                       P2 | 0.  1.  0.
C                          |  /                           |
C                          | /                         P3 |-1.  0.  0.
C                          |/                             |
C   ----------(-1.)--------|-------------1.-------->X  P4 | 0. -1.  0.
C                         /                               |
C                        /                             P5 | 0.  0.  1.
C                       /
C                      /
C                    -1.
C*********************************************************************
      CHARACTER*4 NOM2
      REAL*8 XREF(ND,NP),X(NPG),Y(NPG),Z(NPG)
      DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG)
      DIMENSION FM(MP,NPG),GM(ND,MP,NPG)
      REAL*8 A,B,C,D,U(5),H(5)
      REAL*8 CO,AL,BE,DIV
 
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      XZERO=0.0D0
 
      XREF(1,1)=1.D0
      XREF(2,1)=0.D0
      XREF(3,1)=0.D0
 
      XREF(1,2)=0.D0
      XREF(2,2)=1.D0
      XREF(3,2)=0.D0
 
      XREF(1,3)=-1.D0
      XREF(2,3)=0.D0
      XREF(3,3)=0.D0
 
      XREF(1,4)=0.D0
      XREF(2,4)=-1.D0
      XREF(3,4)=0.D0
 
      XREF(1,5)=0.D0
      XREF(2,5)=0.D0
      XREF(3,5)=1.D0
 
C Verification des coordonnées
      IF(.FALSE.)THEN
C      IF(.TRUE.)THEN
      DO 11 L=1,NP
      X(L)=XREF(1,L)
      Y(L)=XREF(2,L)
      Z(L)=XREF(3,L)
 11   CONTINUE
 
      DO 12 L=1,NP
      DZEM = 1.D0 - Z(L) + 1.D-10
      DZEM4= (1.0D0 - Z(L))*4.0D0 + 1.D-10
      AUX  = X(L)+Y(L)+Z(L)-1.0D0
      AUX1 =-X(L)+Y(L)+Z(L)-1.0D0
      AUX2 =-X(L)-Y(L)+Z(L)-1.0D0
      AUX3 = X(L)-Y(L)+Z(L)-1.0D0
C
      FN(1,L)=AUX1*AUX2/DZEM4
      FN(2,L)=AUX2*AUX3/DZEM4
      FN(3,L)=AUX*AUX3/DZEM4
      FN(4,L)=AUX*AUX1/DZEM4
      FN(5,L)=Z(L)
      write(6,1033)L,FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L)
 12   CONTINUE
 1033 FORMAT(1X,I4,' FN',10(1X,1PD11.4))
      ENDIF
C Fin Vérification
 
      IF(NPG.EQ.1)THEN
      X(1)=0.D0
      Y(1)=0.D0
      Z(1)=0.25D0
      PG(1)=2.D0/3.D0
      ENDIF
      IF(NPG.EQ.5)THEN
      AUX=0.5D0
      H1 = 0.1531754163448146D0
      H2 = 0.6372983346207416D0
      PO=2.D0 / 15.D0
C
      X(1)=AUX
      Y(1)=XZERO
      Z(1)=H1
      PG(1)=PO
C
      X(2)=XZERO
      Y(2)=-AUX
      Z(2)=H1
      PG(2)=PO
C
      X(3)=AUX
      Y(3)=XZERO
      Z(3)=H1
      PG(3)=PO
C
      X(4)=XZERO
      Y(4)=-AUX
      Z(4)=H1
      PG(4)=PO
C
      X(5)=XZERO
      Y(5)=XZERO
      Z(5)=H2
      PG(5)=PO
C
      ENDIF
      DO 1 L=1,NPG
      DZEM = 1.D0 - Z(L)
      DZEM4= (1.0D0 - Z(L))*4.0D0
      AUX  = X(L)+Y(L)+Z(L)-1.0D0
      AUX1 =-X(L)+Y(L)+Z(L)-1.0D0
      AUX2 =-X(L)-Y(L)+Z(L)-1.0D0
      AUX3 = X(L)-Y(L)+Z(L)-1.0D0
C
      FN(1,L)=AUX1*AUX2/DZEM4
      FN(2,L)=AUX2*AUX3/DZEM4
      FN(3,L)=AUX*AUX3/DZEM4
      FN(4,L)=AUX*AUX1/DZEM4
      FN(5,L)=Z(L)
C
      GR(1,1,L)=(-AUX1-AUX2)/DZEM4
      GR(2,1,L)=( AUX2-AUX1)/DZEM4
      GR(3,1,L)=(AUX1+AUX2+AUX1*AUX2/DZEM)/DZEM4
C
      GR(1,2,L)=( AUX2-AUX3)/DZEM4
      GR(2,2,L)=(-AUX2-AUX3)/DZEM4
      GR(3,2,L)=(AUX2+AUX3+AUX2*AUX3/DZEM)/DZEM4
C
      GR(1,3,L)=( AUX +AUX3)/DZEM4
      GR(2,3,L)=( AUX3-AUX )/DZEM4
      GR(3,3,L)=(AUX3+AUX +AUX3*AUX /DZEM)/DZEM4
C
      GR(1,4,L)=( AUX1-AUX )/DZEM4
      GR(2,4,L)=( AUX +AUX1)/DZEM4
      GR(3,4,L)=( AUX+AUX1+ AUX*AUX1/DZEM)/DZEM4
C
      GR(1,5,L)=XZERO
      GR(2,5,L)=XZERO
      GR(3,5,L)=1.0D0
C
 
      IF(NOM2.EQ.'P1P1')THEN
      FM(1,L)=FN(1,L)
      FM(2,L)=FN(2,L)
      FM(3,L)=FN(3,L)
      FM(4,L)=FN(4,L)
      FM(5,L)=FN(5,L)
C
      GM(1,1,L)=(-AUX1-AUX2)/DZEM4
      GM(2,1,L)=( AUX2-AUX1)/DZEM4
      GM(3,1,L)=(AUX1+AUX2+AUX1*AUX2/DZEM)/DZEM4
C
      GM(1,2,L)=( AUX2-AUX3)/DZEM4
      GM(2,2,L)=(-AUX2-AUX3)/DZEM4
      GM(3,2,L)=(AUX2+AUX3+AUX2*AUX3/DZEM)/DZEM4
C
      GM(1,3,L)=( AUX +AUX3)/DZEM4
      GM(2,3,L)=( AUX3-AUX )/DZEM4
      GM(3,3,L)=(AUX3+AUX +AUX3*AUX /DZEM)/DZEM4
C
      GM(1,4,L)=( AUX1-AUX )/DZEM4
      GM(2,4,L)=( AUX +AUX1)/DZEM4
      GM(3,4,L)=( AUX+AUX1+ AUX*AUX1/DZEM)/DZEM4
C
      GM(1,5,L)=XZERO
      GM(2,5,L)=XZERO
      GM(3,5,L)=1.0D0
C
      ELSE
      FM(1,L)=1.D0
      GM(1,1,L)=0.D0
      GM(2,1,L)=0.D0
      GM(3,1,L)=0.D0
      ENDIF
 1    CONTINUE
 
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C     WRITE(6,100)
C     WRITE(6,101)
C     WRITE(6,1002)FN
C     WRITE(6,1002)GR
C     WRITE(6,101)
      RETURN
 1003 FORMAT(1X,' FN',10(1X,1PD11.4))
 1002 FORMAT(1X,'  GR',3(1X,1PD11.4))
 1001 FORMAT(20(1X,I5))
 100  FORMAT(1H1)
 101  FORMAT(1X,'... SUB PB503 ... FN,GR,FOM,GM ',9(10H..........)/)
      END