C FFFACE    SOURCE    CHAT      05/01/12    23:58:45     5004
      SUBROUTINE FFFACE(XA,NPF)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C************************************************************************
C
C
C
C************************************************************************
      DIMENSION FN(20),XA(3,21),COOR(3,7)
C
      IF(NPF.EQ.8)THEN
      X=0.5D0
      Y=0.5D0

      FN(1)=-2.D0*(1.D0-X)*(1.D0-Y)*(X+Y-0.5D0)
      FN(2)= 4.D0*X*(1.D0-X)*(1.D0-Y)
      FN(3)=-2.D0*X*(1.D0-Y)*(Y-X+0.5D0)
      FN(4)= 4.D0*X*Y*(1.D0-Y)
      FN(5)=-2.D0*X*Y*(1.5D0-X-Y)
      FN(6)= 4.D0*X*Y*(1.D0-X)
      FN(7)=-2.D0*Y*(1.D0-X)*(X-Y+0.5D0)
      FN(8)= 4.D0*Y*(1.D0-Y)*(1.D0-X)

      DO 2 N=1,3
      XA(N,NPF+1)=
     &         FN(1)*XA(N,1)+FN(2)*XA(N,2)+FN(3)*XA(N,3)+FN(4)*XA(N,4)
     &+        FN(5)*XA(N,5)+FN(6)*XA(N,6)+FN(7)*XA(N,7)+FN(8)*XA(N,8)
 2    CONTINUE

      ELSEIF(NPF.EQ.6)THEN
      A=SQRT(2.D0)
      AS2=A/2.D0
      X=AS2
      Y=AS2
      FN(1)=(X+Y-AS2)*(X+Y-A)
      FN(2)=-2.D0*X*(X+Y-A)
      FN(3)=-X*(X-AS2)
      FN(4)=2.D0*X*Y
      FN(5)=Y*(Y-AS2)
      FN(6)=-2.D0*Y*(X+Y-A)

      DO 4 N=1,3
      XA(N,NPF+1)=(XA(N,1)+XA(N,3)+XA(N,5))/3.D0
C?    XA(N,NPF+1)=
C?   &         FN(1)*XA(N,1)+FN(2)*XA(N,2)+FN(3)*XA(N,3)+FN(4)*XA(N,4)
C?   &+        FN(5)*XA(N,5)+FN(6)*XA(N,6)
 4    CONTINUE

      ELSEIF(NPF.EQ.20)THEN

      X=0.5D0
      Y=0.5D0
      Z=0.5D0

      FN(1)=-2.D0*(1.D0-X)*(1.D0-Y)*(1.D0-Z)*(X+Y+Z-0.5D0)
      FN(2)= 4.D0*X*(1.D0-X)*(1.D0-Y)*(1.D0-Z)
      FN(3)=-2.D0*X*(1.D0-Y)*(1.D0-Z)*(Z+Y-X+0.5D0)
      FN(4)= 4.D0*X*Y*(1.D0-Y)*(1.D0-Z)
      FN(5)=-2.D0*X*Y*(1.D0-Z)*(1.5D0-X-Y+Z)
      FN(6)= 4.D0*X*Y*(1.D0-X)*(1.D0-Z)
      FN(7)=-2.D0*Y*(1.D0-X)*(1.D0-Z)*(X-Y+Z+0.5D0)
      FN(8)= 4.D0*Y*(1.D0-Y)*(1.D0-X)*(1.D0-Z)
      FN(9) = 4.D0*Z*(1.D0-Z)*(1.D0-X)*(1.D0-Y)
      FN(10)= 4.D0*X*Z*(1.D0-Z)*(1.D0-Y)
      FN(11)= 4.D0*X*Y*Z*(1.D0-Z)
      FN(12)= 4.D0*Y*Z*(1.D0-Z)*(1.D0-X)
      FN(13)=-2.D0*(1.D0-X)*(1.D0-Y)*Z*(X+Y-Z+0.5D0)
      FN(14)= 4.D0*X*(1.D0-X)*(1.D0-Y)*Z
      FN(15)=-2.D0*X*(1.D0-Y)*Z*(-Z+Y-X+1.5D0)
      FN(16)= 4.D0*X*Y*(1.D0-Y)*Z
      FN(17)=-2.D0*X*Y*Z*(2.5D0-X-Y-Z)
      FN(18)= 4.D0*X*Y*(1.D0-X)*Z
      FN(19)=-2.D0*Y*(1.D0-X)*Z*(X-Y-Z+1.5D0)
      FN(20)= 4.D0*Y*(1.D0-Y)*(1.D0-X)*Z

      DO 20 N=1,3
      XA(N,NPF+1)=
     &         FN(1)*XA(N,1)+FN(2)*XA(N,2)+FN(3)*XA(N,3)+FN(4)*XA(N,4)
     &+        FN(5)*XA(N,5)+FN(6)*XA(N,6)+FN(7)*XA(N,7)+FN(8)*XA(N,8)
     &+  FN(9)*XA(N,9)+FN(10)*XA(N,10)+FN(11)*XA(N,11)+FN(12)*XA(N,12)
     &+FN(13)*XA(N,13)+FN(14)*XA(N,14)+FN(15)*XA(N,15)+FN(16)*XA(N,16)
     &+FN(17)*XA(N,17)+FN(18)*XA(N,18)+FN(19)*XA(N,19)+FN(20)*XA(N,20)
 20   CONTINUE

      ELSEIF(NPF.EQ.15)THEN

      A=SQRT(2.D0)
      AS2=A*0.5D0
      CA2=1.D0
      CA4=2.D0
      A4=4.D0/A

C     COOR(1,1)=AS2
C     COOR(2,1)=0.D0
C     COOR(3,1)=0.5D0

      COOR(1,2)=AS2
      COOR(2,2)=AS2
      COOR(3,2)=0.5D0

C     COOR(1,3)=0.D0
C     COOR(2,3)=AS2
C     COOR(3,3)=0.5D0

      X=COOR(1,2)
      Y=COOR(2,2)
      Z=COOR(3,2)
C     write(6,*)' X,Y,Z=',X,Y,Z

      FN(1)= CA2*(X+Y-A)*(1.D0-Z)*(X+Y+A*Z-AS2)
      FN(2)=-CA4*X*(X+Y-A)*(1.D0-Z)
      FN(3)=-CA2*X*(1.D0-Z)*(-X+A*Z+AS2)
      FN(4)= CA4*X*Y*(1.D0-Z)
      FN(5)=-CA2*Y*(1.D0-Z)*(-Y+A*Z+AS2)
      FN(6)=-CA4*Y*(X+Y-A)*(1.D0-Z)
      FN(7)=-A4*Z*(1.D0-Z)*(X+Y-A)
      FN(8)= A4*X*Z*(1.D0-Z)
      FN(9) = A4*Y*Z*(1.D0-Z)
      FN(10)= CA2*(X+Y-A)*Z*(X+Y+A*(1.D0-Z)-AS2)
      FN(11)=-CA4*X*(X+Y-A)*Z
      FN(12)=-CA2*X*Z*(-X+A*(1.D0-Z)+AS2)
      FN(13)= CA4*X*Y*Z
      FN(14)=-CA2*Y*Z*(-Y+A*(1.D0-Z)+AS2)
      FN(15)=-CA4*Y*Z*(X+Y-A)

      DO 15 N=1,3
      XA(N,NPF+1)=(XA(N,7)+XA(N,8)+XA(N,9))/3.D0
C?    XA(N,NPF+1)=
C?   &         FN(1)*XA(N,1)+FN(2)*XA(N,2)+FN(3)*XA(N,3)+FN(4)*XA(N,4)
C?   &+        FN(5)*XA(N,5)+FN(6)*XA(N,6)+FN(7)*XA(N,7)+FN(8)*XA(N,8)
C?   &+  FN(9)*XA(N,9)+FN(10)*XA(N,10)+FN(11)*XA(N,11)+FN(12)*XA(N,12)
C?   &+FN(13)*XA(N,13)+FN(14)*XA(N,14)+FN(15)*XA(N,15)
 15   CONTINUE
      ELSE

      DO 30 N=1,3
      U=0.D0
      DO 31 I=1,NPF
      U=U+XA(N,I)
 31   CONTINUE
      XA(N,NPF+1)=U/FLOAT(NPF)
 30   CONTINUE

      ENDIF



1002  format(10(1x,1pe11.4))
      RETURN
      END


