Sources Esope | Cast3M

pb301
C PB301     SOURCE    MAGN      10/05/18    21:16:31     6675
      SUBROUTINE PB301(XREF,X,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 : SEG3 Element quadratique Q2
C                                         : SEG3 Element macro     2xQ1
C                                           (pression Q1 dans les 2 cas)
C                                           Cf VNIMP pression continue
C ITYPI=0 1 pts d'integration de type Gauss (sous integration)
C ITYPI=1 pts d'integration de type Gauss (en nombre suffisant)
C ITYPI=2 pts d'integration de type Gauss Lobatto i.e. sommets element
C
C REMARQUE :
C En coordonnées cylindriques si le point milieux (2) est
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 le point mileux sur l'élément
C de référence, alfa=0.50000000000005 au lieu de 0.5 .
C Voir aussi pb902.eso et pb702.eso.
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 orientation correcte de l'élément courant.
C************************************************************************
 
      DIMENSION X(NPG),XREF(ND,NP)
      DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG)
      REAL*8 FM(MP,NPG),GM(ND,MP,NPG)
      DIMENSION U(5),H(5)
      CHARACTER*4 NOM2
 
-INC PPARAM
-INC CCOPTIO
 
      XREF(1,1)=0.D0
      XREF(1,2)=0.5D0
      XREF(1,3)=1.D0
 
      IF(IFOMOD.EQ.0)THEN
      alfa=0.50000000000005D0
      ELSE
      alfa=0.5D0
      ENDIF
 
      h1=1.D0/alfa
      h2=-1.D0/alfa/(alfa - 1.D0)
      h3=1.D0/(alfa - 1.D0)
      CALL CALUHG(U,H,NG)
      A=0.D0
      B=1.D0
      IF(ITYPI.EQ.2)THEN
      X(1)=0.D0
      X(2)=alfa
      X(3)=1.D0
      DO 4 L=1,3
      PG(L)=1.D0/3.D0
 4    CONTINUE
      ELSE
      CALL CALG1(A,B,NG,H,U,X,PG)
      ENDIF
 
C     write(6,*)'Xg=',(X(I),I=1,NPG)
C     write(6,*)'Pg=',(Pg(I),I=1,NPG)
 
      IF(NOM2.NE.'MCF1')THEN
      DO 1 L=1,NPG
      FN(1,L)=h1*(X(L)-1.D0)*(X(L)-alfa)
      FN(2,L)=h2*X(L)*(1.D0-X(L))
      FN(3,L)=h3*X(L)*(alfa-X(L))
      GR(1,1,L)=h1*(2.D0*X(L)- 1.D0 - alfa)
      GR(1,2,L)=h2*(1.D0-2.D0*X(L))
      GR(1,3,L)=h3*(alfa-2.D0*X(L))
 1    CONTINUE
      ELSE
C Cas MACRO ELEMENT
      DO 10 L=1,NPG
        IF(X(L).LE.0.5D0)THEN
        FN(1,L)=1.D0-2.D0*X(L)
        FN(2,L)=2.D0*X(L)
        FN(3,L)=0.D0
        GR(1,1,L)=-2.D0
        GR(1,2,L)=2.D0
        GR(1,3,L)=0.D0
        ELSE
        FN(1,L)=0.D0
        FN(2,L)=2.D0-2.D0*X(L)
        FN(3,L)=2.D0*X(L)-1.D0
        GR(1,1,L)=0.D0
        GR(1,2,L)=-2.D0
        GR(1,3,L)=2.D0
        ENDIF
 10   CONTINUE
      ENDIF
 
      IF(MP.EQ.1)THEN
      DO 2 L=1,NPG
      FM(1,L)=1.D0
      GM(1,1,L)=0.D0
 2    CONTINUE
      ELSEIF(MP.EQ.2)THEN
 
      FM(1,1)= 1.D0
      FM(1,2)= 0.D0
 
      FM(2,1)= 0.D0
      FM(2,2)= 1.D0
 
      CALL INITD(GM,(ND*MP*NPG),0.D0)
 
      ENDIF
      IF(NOM2.EQ.'P1P1'.OR.NOM2.EQ.'PFP1'.OR.NOM2.EQ.'MCF1')THEN
      DO 3 L=1,NPG
C
      FM(1,L)=1.D0-X(L)
      FM(2,L)=X(L)
      GM(1,1,L)=-1.D0
      GM(1,2,L)=1.D0
 3    CONTINUE
 
      ENDIF
 
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C     WRITE(6,101)
C     DO 104 L=1,NPG
C     write(6,*)'FN npg=',L
C     WRITE(6,1002)(FN(I,L),I=1,NP)
C     write(6,*)'GR(1 npg=',L
C     WRITE(6,1002)(GR(1,I,L),I=1,NP)
C     write(6,*)'GR(2 npg=',L
C     WRITE(6,1002)(GR(2,I,L),I=1,NP)
C     write(6,*)'FM npg=',L
C     WRITE(6,1002)(FM(I,L),I=1,MP)
C     WRITE(6,1002)FM
C     WRITE(6,1002)GM
C104  CONTINUE
C     WRITE(6,101)
      RETURN
 1002 FORMAT(10(1X,1PD11.4))
 1001 FORMAT(20(1X,I5))
 101  FORMAT(1X,'... SUB PB301 ... FN,GR ',9(10H..........)/)
C
      END