Sources Esope | Cast3M

pb702
C PB702     SOURCE    MAGN      10/05/31    21:15:16     6679
      SUBROUTINE PB702(XREF,X,Y,PG,FN,GR,FM,GM,ND,NP,MP,NPG,NOM2,ITYPI)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C
C     CALCULE LES FONCTIONS DE FORME D'UN : TRI7
C
C Cf Crouzeix Raviard. Conforming and nonconforming finite element
C methods for solving the stationary Stokes equations I
C R.A.I.R.O. (7e année, décembre 1973,R-3,p.33 à 76)
C
C
C       ^ eta
C       |
C       |n3
C  r2   5                              1 -> A1
C       |\                             2 -> A12  -
C       | \    r2=sqrt(2)              3 -> A2
C       |  \                           4 -> A23
C       |   \                          5 -> A3
C       |    \                         6 -> A13
C       |     \                        7 -> A123
C r2/2  6------4
C       |      |\                       D1=1-x/r2-y/r2
C       |      | \                      D2=x/r2
C       |   7  |  \                     D3=y/r2
C       |      |   \
C       |      |    \
C       1------2-----3--->ksi
C       0     r2/2   r2
C
C  Pi  = Di(aDi-c)+FD1xD2xD3       1 &lt;= i &lt;= 3
C  Pij = bxDi(aDi-c)-GxD1xD2xD3    1 &lt;= i &lt; j &lt;= 3
C  P123= QxD1xD2xD3
C  Pour l'élément de référence on a a=2 b=4 c=1 F=3 G=12
C  Si on veut perturber un peu la position des points il faut néanmoins
C  vérifier les propriétés essentielles des fonctions de forme qui outre
C  Pi=1 en X=i et Pi=0 ailleur on doit avoir Somme des Pi = 1
C (intégration exacte d'une fonction constante). Ceci conduit à des
C  relations entre a,b,c,F,G et Q
C  b=2a c=a-1 Q=3(G-F) (Somme Pi = 1) a ne 2 points Aij non au milieux
C des arètes
C
C ITYPI=2 pts d'integration de type Gauss Lobatto i.e. sommets element
C
C REMARQUE :
C En coordonnées cylindriques,pour une distribution régulière des noeuds
C conforme à l'article de Crouzeix Raviard, la matrice masse diagonale
C s'annule sur l'axe, comme pour le QUA9 (voir pb902.eso).
C /1
C | 2 pi x N1 dx = 0.
C /0
C Pour remédier à cela on se sert de la bulle en la décalant un peu par
C rapport à l'élément de référence en jouant sur les constantes F et G.
C On a pas à toucher "a" position des points 2,4 et 6.
C Le choix que l'on a fait conduit à une matrice masse diagonale >0. En
C conséquence ce choix doit être accompagné d'une numérotation dans le
C sens trigo de l'élément courant.
C
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)
      REAL*8 A,B,C,D,U(5),H(5)
      CHARACTER*4 NOM2
 
-INC PPARAM
-INC CCOPTIO
C***
 
      IF(IFOMOD.EQ.0)THEN
      F=3.0001
      G=11.9999
      F=3.
      G=12.
      ELSE
      F=3.
      G=12.
      ENDIF
      Q=3.*(G-F)
      a=2.
      b=2.*a
      c=a-1.
 
 
      R2=SQRT(2.D0)
      IPOSP=2
 
 
      XREF(1,1)=0.D0
      XREF(2,1)=0.D0
 
      XREF(1,2)=R2/2.D0
      XREF(2,2)=0.D0
 
      XREF(1,3)=R2
      XREF(2,3)=0.D0
 
      XREF(1,4)=R2/2.D0
      XREF(2,4)=R2/2.D0
 
      XREF(1,5)=0.D0
      XREF(2,5)=R2
 
      XREF(1,6)=0.
      XREF(2,6)=R2/2.D0
 
      XREF(1,7)=R2/3.D0
      XREF(2,7)=R2/3.D0
 
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
 
      D1=1.D0-X(L)/R2-Y(L)/R2
      D2=X(L)/R2
      D3=Y(L)/R2
 
      FN(1,L)=D1*(a*D1-c)+F*D1*D2*D3
      FN(2,L)=b*D1*D2-G*D1*D2*D3
      FN(3,L)=D2*(a*D2-c)+F*D1*D2*D3
      FN(4,L)=b*D2*D3-G*D1*D2*D3
      FN(5,L)=D3*(a*D3-c)+F*D1*D2*D3
      FN(6,L)=b*D1*D3-G*D1*D2*D3
      FN(7,L)=Q*D1*D2*D3
      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)
 12   CONTINUE
 1033 FORMAT(1X,I4,' FN',10(1X,1PD11.4))
      ENDIF
C Fin Vérification
 
      IF(ITYPI.EQ.2)THEN
      X(1)=0.D0
      Y(1)=0.D0
      X(2)=R2/2.D0
      Y(2)=0.D0
      X(3)=R2
      Y(3)=0.D0
      X(4)=R2/2.D0
      Y(4)=R2/2.D0
      X(5)=0.D0
      Y(5)=R2
      X(6)=0.D0
      Y(6)=R2/2.D0
      X(7)=(X(1)+X(3)+X(5))/3.D0
      Y(7)=(Y(1)+Y(3)+Y(5))/3.D0
      DO 2 L=1,7
      PG(L)=1.D0/7.D0
 2    CONTINUE
      ELSE
      CALL CALUHH(X,Y,PG,NPG)
      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.'PRP1')THEN
C?    FM(1,L)=R2*(X(L)+Y(L)-R2/2.D0)
C?    FM(2,L)=-R2*(X(L)-R2/2.D0)
C?    FM(3,L)=-R2*(Y(L)-R2/2.D0)
      A1=6.D0/5.D0/R2
      FM(1,L)=5.D0/3.D0*(1.D0-A1*X(L)-A1*Y(L))
      FM(2,L)=-1.D0/3.D0*(1.D0-6.D0/R2*X(L))
      FM(3,L)=-1.D0/3.D0*(1.D0-6.D0/R2*Y(L))
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)=-R2/2.D0*(X(L)+Y(L)-R2)
      FM(2,L)=R2/2.D0*X(L)
      FM(3,L)=R2/2.D0*Y(L)
 
      GM(1,1,L)=-R2/2.D0
      GM(2,1,L)=-R2/2.D0
      GM(1,2,L)=R2/2.D0
      GM(2,2,L)=0.D0
      GM(1,3,L)=0.D0
      GM(2,3,L)=R2/2.D0
      ENDIF
 
C
c     IF(.FALSE.)THEN
      IF(.TRUE.)THEN
C
 
      D1=1.D0-X(L)/R2-Y(L)/R2
      D2=X(L)/R2
      D3=Y(L)/R2
 
      FN(1,L)=D1*(a*D1-c)+F*D1*D2*D3
      FN(2,L)=b*D1*D2-G*D1*D2*D3
      FN(3,L)=D2*(a*D2-c)+F*D1*D2*D3
      FN(4,L)=b*D2*D3-G*D1*D2*D3
      FN(5,L)=D3*(a*D3-c)+F*D1*D2*D3
      FN(6,L)=b*D1*D3-G*D1*D2*D3
      FN(7,L)=Q*D1*D2*D3
c     write(6,*)FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L),FN(7,L)
 
      D1x=-1./R2
      D1y=-1./R2
 
      D2x=1./R2
      D2y=0.
 
      D3x=0.
      D3y=1/R2
 
      GR(1,1,L)=D1x*(a*D1-c)+D1*a*D1x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
      GR(2,1,L)=D1y*(a*D1-c)+D1*a*D1y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
 
      GR(1,2,L)=b*(D1x*D2+D1*D2x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
      GR(2,2,L)=b*(D1y*D2+D1*D2y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
 
      GR(1,3,L)=D2x*(a*D2-c)+D2*a*D2x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
      GR(2,3,L)=D2y*(a*D2-c)+D2*a*D2y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
 
      GR(1,4,L)=b*(D2x*D3+D2*D3x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
      GR(2,4,L)=b*(D2y*D3+D2*D3y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
 
      GR(1,5,L)=D3x*(a*D3-c)+D3*a*D3x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
      GR(2,5,L)=D3y*(a*D3-c)+D3*a*D3y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
 
      GR(1,6,L)=b*(D1x*D3+D1*D3x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
      GR(2,6,L)=b*(D1y*D3+D1*D3y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
 
      GR(1,7,L)=Q*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x)
      GR(2,7,L)=Q*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y)
 
      ELSE
 
c     write(6,*)'La base'
 
C
C
      FN(1,L)=(1.D0-X(L)/R2-Y(L)/R2)*(1.D0-2.D0*X(L)/R2-2.D0*Y(L)/R2+
     & 3.D0/2.D0*X(L)*Y(L))
      FN(2,L)=(1.D0-X(L)/R2-Y(L)/R2)*(4.D0/R2-6.D0*Y(L))*X(L)
      FN(3,L)=X(L)/R2*(2.D0*X(L)/R2-1.D0+3.D0*(1.D0-X(L)/R2-Y(L)/R2)
     & *Y(L)/R2)
      FN(4,L)=2.D0*X(L)*Y(L)*(1.D0-3.D0*(1.D0-X(L)/R2-Y(L)/R2))
      FN(5,L)=Y(L)*Y(L)-Y(L)*R2/2.D0+3.D0/2.D0*X(L)*Y(L)-
     & 3.D0*R2/4.D0*X(L)*X(L)*Y(L)-3.D0*R2/4.D0*X(L)*Y(L)*Y(L)
      FN(6,L)=2.D0*R2*(-R2/2.D0*Y(L)*Y(L)+3.D0/2.D0*(Y(L)*Y(L)*X(L)+
     & X(L)*X(L)*Y(L))-2.D0*R2*X(L)*Y(L)+Y(L))
      FN(7,L)=27.D0/2.D0*(X(L)*Y(L)-R2/2.D0*X(L)*X(L)*Y(L)-
     &R2/2.D0*X(L)*Y(L)*Y(L))
C
C
C
      GR(1,1,L)=-3.D0*R2/4.D0*Y(L)*Y(L)+2.D0*X(L)+7./2.D0*Y(L)-
     * 3.D0*R2/2.D0*X(L)*Y(L)-3.D0*R2/2.D0
      GR(2,1,L)=-3.D0*R2/4.D0*X(L)*X(L)+2.D0*Y(L)+7./2.D0*X(L)-
     * 3.D0*R2/2.D0*X(L)*Y(L)-3.D0*R2/2.D0
      GR(1,2,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-4.D0*X(L)-8.D0*Y(L)+
     * 2.D0*R2
      GR(2,2,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-8.D0*X(L)
      GR(1,3,L)=-3.D0*R2/4.D0*Y(L)*Y(L)-3.D0*R2/2.D0*X(L)*Y(L)+
     * 2.D0*X(L)+3.D0/2.D0*Y(L)-R2/2.D0
      GR(2,3,L)=R2/2.D0*(-3.D0/2.D0*X(L)*X(L)-3.D0*X(L)*Y(L)+
     * 3.D0*R2/2.D0*X(L))
      GR(1,4,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-4.D0*Y(L)
      GR(2,4,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-4.D0*X(L)
      GR(1,5,L)=-3.D0*R2/4.D0*Y(L)*Y(L)-3.D0*R2/2.D0*X(L)*Y(L)
     & +3.D0/2.D0*Y(L)
      GR(2,5,L)=-3.D0*R2/4.D0*X(L)*X(L)-3.D0*R2/2.D0*X(L)*Y(L)
     & +2.D0*Y(L)+
     * 3.D0/2.D0*X(L)-R2/2.D0
      GR(1,6,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-8.D0*Y(L)
      GR(2,6,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-8.D0*X(L)
     & -4.D0*Y(L)+2.D0*R2
      GR(1,7,L)=27.D0/2.D0*(-R2/2*Y(L)*Y(L)-R2*X(L)*Y(L)+Y(L))
      GR(2,7,L)=27.D0/2.D0*(-R2/2*X(L)*X(L)-R2*X(L)*Y(L)+X(L))
C
      ENDIF
 1    CONTINUE
 
      IF(.FALSE.)THEN
c     IF(.TRUE.)THEN
      DO 121 L=1,NPG
      F1=0.
      DO 131 I=1,NP
      F1=F1+FN(I,L)
 131  CONTINUE
      Write(6,*)' L=',L,F1
 121  CONTINUE
 
      DO 122 L=1,NPG
      G1=0.
      G2=0.
      DO 141 I=1,NP
      G1=G1+GR(1,I,L)
      G2=G2+GR(2,I,L)
 141  CONTINUE
      Write(6,*)' L=',L,G1,G2
 122  CONTINUE
      ENDIF
 
c     stop
 
 
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)
 
      RETURN
 1002 FORMAT(10(1X,1PD11.4))
 1001 FORMAT(20(1X,I5))
 101  FORMAT(1X,'... SUBPB702 ... FM,GM,FN,GR ',9(10H..........)/)
C
      END