Sources Esope | Cast3M

pb902
C PB902     SOURCE    MAGN      10/06/02    21:15:00     6681
      SUBROUTINE PB902
     &(XREF,X,Y,PG,FN,GR,FM,GM,ND,NP,MP,NG,NPG,NOM2,ITYPI)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C************************************************************************
C
C     CALCULE LES FONCTIONS DE FORME D'UN : QUA9
C        eta
C        ^
C        |
C        7      6      5
C 1.     x------x------x
C        |             |
C        |             |
C alfa  8x     9x     4x
C        |             |
C        |             |
C 0.     x------x------x--> ksi
C        1      2      3
C
C        0.     alfa   1.
C
C ITYPI=2 pts d'intégration de type Gauss Lobatto i.e. sommets élément
C
C REMARQUE : Pour l'élément quadratique.
C En coordonnées cylindriques si les points milieux (2,4,6,8,9) sont
C pile poil au milieu, la matrice masse diagonale s'annule sur l'axe.
C /1
C | 2 pi x N1 dx = 0.
C /0
C Pour remédier à cela on décale un peu les points mileux sur l'élément
C de référence, alfa=0.50005 au lieu de 0.5 .
C Si on rapprochait les points milieux de l'axe (alfa = 0.45 par exemple)
C la matrice masse diagonale serait négative sur l'axe, ce que l'on ne
C veut pas non plus. En conséquence le choix qu'on a fait doit être
C accompagné d'une numérotation dans le sens trigo de l'élément courant.
C************************************************************************
      REAL*8 XREF(ND,NP),X(NPG),Y(NPG)
      DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG)
      DIMENSION FM(MP,NPG),GM(ND,MP,NPG)
      CHARACTER*4 NOM2
      REAL*8 A,B,C,D,U(5),H(5),AUX1,AUX2
 
-INC PPARAM
-INC CCOPTIO
 
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      AUX1=1.D0/2.D0+3.D0**(1.D0/2.D0)/6.D0
      AUX2=1.D0/2.D0-3.D0**(1.D0/2.D0)/6.D0
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      IF(IFOMOD.EQ.0)THEN
      alfa=0.50005D0
      ELSE
      alfa=0.5D0
      ENDIF
 
      alfai=1.D0/alfa
      alfa2=1.D0/alfa/(alfa - 1.D0)
      alfa3=1.D0/(alfai - 1.D0)
 
      XREF(1,1)=0.D0
      XREF(2,1)=0.D0
 
      XREF(1,2)=alfa
      XREF(2,2)=0.D0
 
      XREF(1,3)=1.D0
      XREF(2,3)=0.D0
 
      XREF(1,4)=1.D0
      XREF(2,4)=alfa
 
      XREF(1,5)=1.D0
      XREF(2,5)=1.D0
 
      XREF(1,6)=alfa
      XREF(2,6)=1.D0
 
      XREF(1,7)=0.D0
      XREF(2,7)=1.D0
 
      XREF(1,8)=0.D0
      XREF(2,8)=alfa
 
      XREF(1,9)=alfa
      XREF(2,9)=alfa
 
C Verification des coordonnées
c     IF(.TRUE.)THEN
      IF(.FALSE.)THEN
      DO 11 L=1,NP
      X(L)=XREF(1,L)
      Y(L)=XREF(2,L)
 11   CONTINUE
 
      DO 12 L=1,NP
      FN(1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
     &       *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
      FN(2,L)=alfa2*X(L)*(X(L)-1.D0)
     &       *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
      FN(3,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
     &       *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
      FN(4,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
     &       *alfa2*Y(L)*(Y(L)-1.D0)
      FN(5,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
     &       *alfa3*Y(L)*(alfai*Y(L)-1.D0)
      FN(6,L)=alfa2*X(L)*(X(L)-1.D0)
     &       *alfa3*Y(L)*(alfai*Y(L)-1.D0)
      FN(7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
     &       *alfa3*Y(L)*(alfai*Y(L)-1.D0)
      FN(8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
     &       *alfa2*Y(L)*(Y(L)-1.D0)
      FN(9,L)=alfa2*X(L)*(X(L)-1.D0)
     &       *alfa2*Y(L)*(Y(L)-1.D0)
 
      write(6,1033)L,FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L)
     &              ,FN(7,L),FN(8,L),FN(9,L)
 12   CONTINUE
 1033 FORMAT(1X,I4,' FN',10(1X,1PD11.4))
      ENDIF
C Fin Vérification
 
      CALL CALUHG(U,H,NG)
C
      A=0.D0
      B=1.D0
      C=0.D0
      D=1.D0
      IF(ITYPI.EQ.2)THEN
      X(1)=0.D0
      Y(1)=0.D0
      X(2)=alfa
      Y(2)=0.D0
      X(3)=1.D0
      Y(3)=0.D0
      X(4)=1.D0
      Y(4)=alfa
      X(5)=1.D0
      Y(5)=1.D0
      X(6)=alfa
      Y(6)=1.D0
      X(7)=0.D0
      Y(7)=1.D0
      X(8)=0.D0
      Y(8)=alfa
      X(9)=alfa
      Y(9)=alfa
      DO 2 L=1,NPG
      PG(L)=1.D0/9.D0
 2    CONTINUE
      ELSE
      CALL CALG2(A,B,C,D,NG,H,U,X,Y,PG)
      ENDIF
 
      DO 1 L=1,NPG
C
      IF(NOM2.EQ.'PRP0')THEN
      FM(1,L)= 1.D0
 
      GM(1,1,L)=0.D0
      GM(2,1,L)=0.D0
 
      ELSEIF(NOM2.EQ.'PRQ1')THEN
      FM(1,L)= 3.D0*(X(L)-AUX1)*(Y(L)-AUX1)
      FM(2,L)=-3.D0*(X(L)-AUX2)*(Y(L)-AUX1)
      FM(3,L)= 3.D0*(X(L)-AUX2)*(Y(L)-AUX2)
      FM(4,L)=-3.D0*(X(L)-AUX1)*(Y(L)-AUX2)
C
      GM(1,1,L)=3.D0*(Y(L)-AUX1)
      GM(2,1,L)=3.D0*(X(L)-AUX1)
      GM(1,2,L)=-3.D0*(Y(L)-AUX1)
      GM(2,2,L)=-3.D0*(X(L)-AUX2)
      GM(1,3,L)=3.D0*(Y(L)-AUX2)
      GM(2,3,L)=3.D0*(X(L)-AUX2)
      GM(1,4,L)=-3.D0*(Y(L)-AUX2)
      GM(2,4,L)=-3.D0*(X(L)-AUX1)
 
      ELSEIF(NOM2.EQ.'PRP1')THEN
C?    FM(1,L)= (3.D0-4.D0*X(L))/2.D0
C?    FM(2,L)= 2.D0*(X(L)-Y(L))
C?    FM(3,L)= (4.D0*Y(L)-1.D0)/2.D0
      FM(1,L)= (4.D0*X(L)+4.D0*Y(L)-3.D0)/3.D0
      FM(2,L)=-1.D0/3.D0*(8.D0*X(L)-4.D0*Y(L)-3.D0)
      FM(3,L)= 1.D0/3.D0*(-8.D0*Y(L)+4.D0*X(L)+3.D0)
C
      GM(1,1,L)=0.D0
      GM(2,1,L)=0.D0
      GM(1,2,L)=0.D0
      GM(2,2,L)=0.D0
      GM(1,3,L)=0.D0
      GM(2,3,L)=0.D0
      ELSEIF(NOM2.EQ.'PFP1')THEN
      FM(1,L)= (X(L)-1.D0)*(Y(L)-1.D0)
      FM(2,L)=-X(L)*(Y(L)-1.D0)
      FM(3,L)= X(L)*Y(L)
      FM(4,L)=-Y(L)*(X(L)-1.D0)
 
      GM(1,1,L)=Y(L)-1.D0
      GM(2,1,L)=X(L)-1.D0
      GM(1,2,L)=-(Y(L)-1.D0)
      GM(2,2,L)=-X(L)
      GM(1,3,L)=Y(L)
      GM(2,3,L)=X(L)
      GM(1,4,L)=-Y(L)
      GM(2,4,L)=-(X(L)-1.D0)
      ENDIF
C
c     N1X=(X(L)-1.D0)*(alfai*X(L)-1.D0)
c     N1Y=(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
c
c     N2X=alfa2*X(L)*(X(L)-1.D0)
c     N2Y=alfa2*Y(L)*(Y(L)-1.D0)
c
c     N3X=alfa3*X(L)*(alfai*X(L)-1.D0)
c     N3Y=alfa3*Y(L)*(alfai*Y(L)-1.D0)
C
 
c     FN(1,L)=N1X*N1Y
      FN(1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
     &       *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
c     FN(2,L)=N2X*N1Y
      FN(2,L)=alfa2*X(L)*(X(L)-1.D0)
     &       *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
c     FN(3,L)=N3X*N1Y
      FN(3,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
     &       *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
c     FN(4,L)=N3X*N2Y
      FN(4,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
     &       *alfa2*Y(L)*(Y(L)-1.D0)
c     FN(5,L)=N3X*N3Y
      FN(5,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
     &       *alfa3*Y(L)*(alfai*Y(L)-1.D0)
c     FN(6,L)=N2X*N3Y
      FN(6,L)=alfa2*X(L)*(X(L)-1.D0)
     &       *alfa3*Y(L)*(alfai*Y(L)-1.D0)
c     FN(7,L)=N1X*N3Y
      FN(7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
     &       *alfa3*Y(L)*(alfai*Y(L)-1.D0)
c     FN(8,L)=N1X*N2Y
      FN(8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
     &       *alfa2*Y(L)*(Y(L)-1.D0)
c     FN(9,L)=N2X*N2Y
      FN(9,L)=alfa2*X(L)*(X(L)-1.D0)
     &       *alfa2*Y(L)*(Y(L)-1.D0)
 
C
c     FN(1,L)=(X(L)-1.D0)*(Y(L)-1.D0)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
c     FN(2,L)=-4.D0*X(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)
c     FN(3,L)=X(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)
c     FN(4,L)=-4.D0*X(L)*Y(L)*(2.D0*X(L)-1.D0)*(Y(L)-1.D0)
c     FN(5,L)=X(L)*Y(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
c     FN(6,L)=-4.D0*X(L)*Y(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)
c     FN(7,L)=Y(L)*(2.D0*Y(L)-1.D0)*(2.D0*X(L)-1.D0)*(X(L)-1.D0)
c     FN(8,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(X(L)-1.D0)*(2.D0*X(L)-1.D0)
c     FN(9,L)=16.D0*X(L)*Y(L)*(X(L)-1.D0)*(Y(L)-1.D0)
C
c     N1X=(X(L)-1.D0)*(alfai*X(L)-1.D0)
c     N1Y=(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
c
c     N2X=alfa2*X(L)*(X(L)-1.D0)
c     N2Y=alfa2*Y(L)*(Y(L)-1.D0)
c
c     N3X=alfa3*X(L)*(alfai*X(L)-1.D0)
c     N3Y=alfa3*Y(L)*(alfai*Y(L)-1.D0)
c
c     G1X=2.D0*alfai*X(L)-(1.D0 + alfai)
c     G1Y=2.D0*alfai*Y(L)-(1.D0 + alfai)
c
c     G2X=alfa2*(2.D0*X(L)-1.D0)
c     G2Y=alfa2*(2.D0*Y(L)-1.D0)
c
c     G3X=alfa3*(2.D0*alfai*X(L)-1.D0)
c     G3Y=alfa3*(2.D0*alfai*Y(L)-1.D0)
 
C
c     GR(1,1,L)=G1X*N1Y
      GR(1,1,L)=(2.D0*alfai*X(L)-(1.D0 + alfai))
     &         *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
c     GR(2,1,L)=N1X*G1Y
      GR(2,1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
     &         *(2.D0*alfai*Y(L)-(1.D0 + alfai))
c     GR(1,2,L)=G2X*N1Y
      GR(1,2,L)=(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
     &         *alfa2*(2.D0*X(L)-1.D0)
c     GR(2,2,L)=N2X*G1Y
      GR(2,2,L)=alfa2*X(L)*(X(L)-1.D0)
     &         *(2.D0*alfai*Y(L)-(1.D0 + alfai))
c     GR(1,3,L)=G3X*N1Y
      GR(1,3,L)=alfa3*(2.D0*alfai*X(L)-1.D0)
     &         *(Y(L)-1.D0)*(alfai*Y(L)-1.D0)
c     GR(2,3,L)=N3X*G1Y
      GR(2,3,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
     &         *(2.D0*alfai*Y(L)-(1.D0 + alfai))
c     GR(1,4,L)=G3X*N2Y
      GR(1,4,L)=alfa3*(2.D0*alfai*X(L)-1.D0)
     &         *alfa2*Y(L)*(Y(L)-1.D0)
c     GR(2,4,L)=N3X*G2Y
      GR(2,4,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
     &         *alfa2*(2.D0*Y(L)-1.D0)
c     GR(1,5,L)=G3X*N3Y
      GR(1,5,L)=alfa3*(2.D0*alfai*X(L)-1.D0)
     &         *alfa3*Y(L)*(alfai*Y(L)-1.D0)
c     GR(2,5,L)=N3X*G3Y
      GR(2,5,L)=alfa3*X(L)*(alfai*X(L)-1.D0)
     &         *alfa3*(2.D0*alfai*Y(L)-1.D0)
c     GR(1,6,L)=G2X*N3Y
      GR(1,6,L)=alfa2*(2.D0*X(L)-1.D0)
     &         *alfa3*Y(L)*(alfai*Y(L)-1.D0)
c     GR(2,6,L)=N2X*G3Y
      GR(2,6,L)=alfa2*X(L)*(X(L)-1.D0)
     &         *alfa3*(2.D0*alfai*Y(L)-1.D0)
c     GR(1,7,L)=G1X*N3Y
      GR(1,7,L)=(2.D0*alfai*X(L)-(1.D0 + alfai))
     &         *alfa3*Y(L)*(alfai*Y(L)-1.D0)
c     GR(2,7,L)=N1X*G3Y
      GR(2,7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
     &         *alfa3*(2.D0*alfai*Y(L)-1.D0)
c     GR(1,8,L)=G1X*N2Y
      GR(1,8,L)=(2.D0*alfai*X(L)-(1.D0 + alfai))
     &         *alfa2*Y(L)*(Y(L)-1.D0)
c     GR(2,8,L)=N1X*G2Y
      GR(2,8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0)
     &         *alfa2*(2.D0*Y(L)-1.D0)
c     GR(1,9,L)=G2X*N2Y
      GR(1,9,L)=alfa2*Y(L)*(Y(L)-1.D0)
     &         *alfa2*(2.D0*X(L)-1.D0)
c     GR(2,9,L)=N2X*G2Y
      GR(2,9,L)=alfa2*X(L)*(X(L)-1.D0)
     &         *alfa2*(2.D0*Y(L)-1.D0)
C
c     GR(1,1,L)=(Y(L)-1.D0)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-3.D0)
c     GR(2,1,L)=(X(L)-1.D0)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-3.D0)
c     GR(1,2,L)=-4.D0*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)
c     GR(2,2,L)=-4.D0*X(L)*(X(L)-1.D0)*(4.D0*Y(L)-3.D0)
c     GR(1,3,L)=(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)*(4.D0*X(L)-1.D0)
c     GR(2,3,L)=X(L)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-3.D0)
c     GR(1,4,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(4.D0*X(L)-1.D0)
c     GR(2,4,L)=-4.D0*X(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
c     GR(1,5,L)=Y(L)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-1.D0)
c     GR(2,5,L)=X(L)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-1.D0)
c     GR(1,6,L)=-4.D0*Y(L)*(2.D0*Y(L)-1.D0)*(2.D0*X(L)-1.D0)
c     GR(2,6,L)=-4.D0*X(L)*(X(L)-1.D0)*(4.D0*Y(L)-1.D0)
c     GR(1,7,L)=Y(L)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-3.D0)
c     GR(2,7,L)=(2.D0*X(L)-1.D0)*(X(L)-1.D0)*(4.D0*Y(L)-1.D0)
c     GR(1,8,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(4.D0*X(L)-3.D0)
c     GR(2,8,L)=-4.D0*(X(L)-1.D0)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)
c     GR(1,9,L)=16.D0*Y(L)*(Y(L)-1.D0)*(2.D0*X(L)-1.D0)
c     GR(2,9,L)=16.D0*X(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)
C
 1    CONTINUE
C
C
C     WRITE(6,101)
C     WRITE(6,1002)FM
C     WRITE(6,1002)GM
C     WRITE(6,1002)FN
C     WRITE(6,1002)GR
C     WRITE(6,101)
 
C     write(6,*)' RET PB902 '
      RETURN
 1002 FORMAT(10(1X,1PD11.4))
 1001 FORMAT(20(1X,I5))
 101  FORMAT(1X,'... SUBPB902 ... FM,GM,FN,GR ',9(10H..........)/)
C
      END