Sources Esope | Cast3M

ffq
C FFQ       SOURCE    CHAT      05/01/12    23:59:09     5004
      SUBROUTINE FFQ(NOMS,XA,COOR,XG,IDIM,NP)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C************************************************************************
C
C
C
C************************************************************************
      DIMENSION FN(20),XA(3,27),COOR(3,NP),XG(3,NP)
      PARAMETER (NBE=18)
      CHARACTER*8 NOMS,LISTE(NBE)
      DATA LISTE/'SEG3    ','TRI6    ','QUA8    ','TRI7    ','QUA9    ',
     &'CU27    ','PR21    ','TE15    ','PY19    ',
     &'PR18    ','PY14    ',
     &'SEG4    ','TR12    ','QU16    ','CU64    ','PR48    ','TE32    ',
     &'PY57    '/
C
      R2=SQRT(2.D0)
 
      CALL OPTLI(IP,LISTE,NOMS,NBE)
C     write(6,*)' FFQ : IP=',ip,' NOMS=',noms,np
      IF(IP.EQ.0)RETURN
      GO TO (301,602,802,602,802,2703,2103,1503),IP
 
 301  CONTINUE
 
C SEG3
      DO 311 L=1,NP
      X=COOR(1,L)
      Y=COOR(2,L)
 
C
      FN(1)=(X-0.5D0)*(X-1.D0)
      FN(2)=-2.D0*X*(X+Y-R2)
      FN(3)= X*(X-R2/2.D0)
C
      DO 321 N=1,IDIM
      XG(N,L)=0.D0
      DO 331 I=1,3
      XG(N,L)=XG(N,L)+FN(I)*XA(N,I)
 331  CONTINUE
 321  CONTINUE
 311  CONTINUE
 
      RETURN
 
 
 
 602  CONTINUE
C TRI6
      DO 612 L=1,NP
      X=COOR(1,L)
      Y=COOR(2,L)
 
C
      FN(1)=(X+Y-R2/2.D0)*(X+Y-R2)
      FN(2)=-2.D0*X*(X+Y-R2)
      FN(3)= X*(X-R2/2.D0)
      FN(4)=2.D0*X*Y
      FN(5)= Y*(Y-R2/2.D0)
      FN(6)=-2.D0*Y*(X+Y-R2)
C
      DO 622 N=1,IDIM
      XG(N,L)=0.D0
      DO 632 I=1,6
      XG(N,L)=XG(N,L)+FN(I)*XA(N,I)
 632  CONTINUE
 622  CONTINUE
 612  CONTINUE
 
      RETURN
 
 
 
 802  CONTINUE
C QUA8
      DO 812 L=1,NP
      X=COOR(1,L)
      Y=COOR(2,L)
C
      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*(X+Y-1.5D0)
      FN(6)= 4.D0*X*Y*(1.D0-X)
      FN(7)= 2.D0*Y*(1.D0-X)*(Y-X-0.5D0)
      FN(8)= 4.D0*Y*(1.D0-X)*(1.D0-Y)
C
      DO 822 N=1,IDIM
      XG(N,L)=0.D0
      DO 832 I=1,8
      XG(N,L)=XG(N,L)+FN(I)*XA(N,I)
 832  CONTINUE
 822  CONTINUE
 812  CONTINUE
      RETURN
 
 2703 CONTINUE
C CU27 -> On ne prend que les 20 premiers points correspondant au CU20
      DO 2712 L=1,NP
      X=COOR(1,L)
      Y=COOR(2,L)
      Z=COOR(3,L)
C
      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)*(Y-X+Z+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)*(0.5D0+X-Y+Z)
      FN(8)= 4.D0*Y*(1.D0-X)*(1.D0-Y)*(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-X)*(1.D0-Z)
      FN(13)=-2.D0*Z*(1.D0-X)*(1.D0-Y)*(X+Y-Z+0.5D0)
      FN(14)= 4.D0*X*Z*(1.D0-X)*(1.D0-Y)
      FN(15)=-2.D0*X*Z*(1.D0-Y)*(-X+Y-Z+1.5D0)
      FN(16)= 4.D0*X*Y*Z*(1.D0-Y)
      FN(17)=-2.D0*X*Y*Z*(-X-Y-Z+2.5D0)
      FN(18)= 4.D0*X*Y*Z*(1.D0-X)
      FN(19)=-2.D0*Y*Z*(1.D0-X)*(X-Y-Z+1.5D0)
      FN(20)= 4.D0*Y*Z*(1.D0-X)*(1.D0-Y)
C
      DO 2722 N=1,IDIM
      XG(N,L)=0.D0
      DO 2732 I=1,20
      XG(N,L)=XG(N,L)+FN(I)*XA(N,I)
 2732 CONTINUE
 2722 CONTINUE
 2712 CONTINUE
      RETURN
 
 2103 CONTINUE
C PR21 -> On ne prend que les 15 premiers points correspondant au PR15
      DO 2112 L=1,NP
      X=COOR(1,L)
      Y=COOR(2,L)
      Z=COOR(3,L)
C
      FN(1)=      (X+Y-R2)*(1.D0-Z)*(X+Y+R2*Z-R2/2.D0)
      FN(2)=-2.D0*X*(X+Y-R2)*(1.D0-Z)
      FN(3)=-X*(1.D0-Z)*(-X+R2*Z+R2/2.D0)
      FN(4)= 2.D0*X*Y*(1.D0-Z)
      FN(5)=-Y*(1.D0-Z)*(-Y+R2*Z+R2/2.D0)
      FN(6)=-2.D0*Y*(1.D0-Z)*(X+Y-R2)
      FN(7)=-4.D0*Z*(1.D0-Z)*(X+Y-R2)/R2
      FN(8)= 4.D0*X*Z*(1.D0-Z)/R2
      FN(9)= 4.D0*Y*Z*(1.D0-Z)/R2
      FN(10)= (X+Y-R2)*Z*(X+Y+R2*(1.D0-Z)-R2/2.D0)
      FN(11)=-2.D0*X*Z*(X+Y-R2)
      FN(12)=-X*Z*(-X+R2*(1.D0-Z)+R2/2.D0)
      FN(13)= 2.D0*X*Y*Z
      FN(14)=-Y*Z*(-Y+R2*(1.D0-Z)+R2/2.D0)
      FN(15)=-2.D0*Y*Z*(X+Y-R2)
C
      DO 2122 N=1,IDIM
      XG(N,L)=0.D0
      DO 2132 I=1,15
      XG(N,L)=XG(N,L)+FN(I)*XA(N,I)
 2132 CONTINUE
 2122 CONTINUE
 2112 CONTINUE
      RETURN
 
 1503 CONTINUE
C TE15 -> On ne prend que les 10 premiers points correspondant au TE10
      A=SQRT(6.D0)
      A2=A*A
      DO 1512 L=1,NP
      X=COOR(1,L)
      Y=COOR(2,L)
      Z=COOR(3,L)
C
      FN(1)= 2.D0*(X+Y+A*Z-A/2.D0)*(X+Y+A*Z-A)/A2
      FN(2)=-4.D0*(X+Y+A*Z-A)*X/A2
      FN(3)= 2.D0*(X-A/2.D0)*X/A2
      FN(4)= 4.D0*X*Y/A2
      FN(5)= 2.D0*(Y-A/2.D0)*Y/A2
      FN(6)=-4.D0*(X+Y+A*Z-A)*Y/A2
      FN(7)=-4.D0*(X+Y+A*Z-A)*Z/A
      FN(8)= 4.D0*X*Z/A
      FN(9)= 4.D0*Y*Z/A
      FN(10)= 2.D0*Z*(Z-0.5D0)
C
      DO 1522 N=1,IDIM
      XG(N,L)=0.D0
      DO 1532 I=1,10
      XG(N,L)=XG(N,L)+FN(I)*XA(N,I)
 1532 CONTINUE
 1522 CONTINUE
 1512 CONTINUE
      RETURN
 
1002  format(10(1x,1pe11.4))
      END