C*********************************************************
C*********************************************************
C*********************************************************
C****
C****
C****   E L E M E N T    F I N I    T E T R A E D R I Q U E
C****
C****   DE CLASSE C1,   DE DEGRE DEUX,   A SEIZE PARAMETRES.
C****
C****   VARIANTE I-C.
C****                                        MARCEL DELEZE 
C****                                        LE CHENE 21 
C****   FEVRIER 1978.                        1633 MARSENS
C****   NOVEMBRE 2008                        SWITZERLAND   
C****                           
C****	
C********************************************************
C********************************************************
C--------------------------------------------------------
C
       FUNCTION PT1 (J, WT,XT,YT,ZT)
C
C--------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DU POLYNOME DE DEGRE TROIS PT(1,J) AU POINT (XT,YT,ZT). 
C**** J PEUT PRENDRE LES VALEURS 0 OU 1.   WT=1.-XT-YT-ZT.
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (J) 10,10,11
   10 PT1=W*(2.-W-2.*(X*X+Y*Y+Z*Z+X*Y+Y*Z+Z*X))
      RETURN
   11 PT1=X*W*(W+0.5*(Y+Z))
      RETURN
      END
C--------------------------------------------------------- 
C
      FUNCTION WT11 (WT,XT,YT,ZT)
C 
C---------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DU POLYNOME PAR MORCEAUX WT(1,1) AU POINT (XT,YT,ZT). 
C**** WT EST INACTIF.
      X=XT
      Y=YT
      Z=ZT
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 20
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 30
C**** MORCEAU 1.
      WT11=X*( 2.*X*X/3.-X*Y-X*Z+2.*Y*Z )
      RETURN
C**** MORCEAU 2.
   20 WT11=Y*Y*( Z-Y/3. )
      RETURN
C**** MORCEAU 3.
   30 WT11=Z*Z*( Y-Z/3. )
      RETURN
      END
C-----------------------------------------------------------
C
      FUNCTION ST10 (WT,XT,YT,ZT)
C
C------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE RT(1,0), RESTREINTE AU
C**** MORCEAU 1, AU POINT (XT,YT,ZT).	WT=1.-XT-YT-ZT.
C**** NOUS SUPPOSONS QUE LE POINT (XT,YT,ZT) NE SOIT PAS SITUE SUR
C**** UNE DES CINQ ARETES DU MORCEAU 1 APPARTENANT AU BORD DU TETRAEDRE. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT
      Z=ZT 
      YMX=Y-X
      ZMX=Z-X
      IF (YMX.GT.EPS .OR. ZMX.GT.EPS) GO TO 10
C**** SUR LA DROITE X=Y=Z, ON A
      ST10=3.*X*X*W
      RETURN
   10 CONTINUE
C**** DANS LE RESTE DU MORCEAU 1, ON A
      VAR=-X*W*( 3./(1.-3.*X) + 2./(Y+ZMX) )
      VAR=-YMX*ZMX*( 6.*X-2.*W-3.+VAR )
      VAR=VAR*( ZMX/(Y*(W+YMX))+YMX/(Z*(W+ZMX)) )/(YMX+ZMX)
      VAR=W*( 2.*(YMX+ZMX)+VAR )/((1.-3.*X)*(Y+ZMX))
      ST10=X*X*W*(3.+VAR)
      RETURN
      END
C-----------------------------------------------------------------
C
      FUNCTION RT10 (WT,XT,YT,ZT)
C
C-----------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE PAR MORCEAUX RT(1,0) AU POINT 
C**** (XT,YT,ZT) DU TETRAEDRE DE REFERENCE.	WT=1.-XT-YT-ZT.
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (X.GT.EPS .AND. Y.GT.EPS .AND.
     1	  Z.GT.EPS .AND. W.GT.EPS ) GO TO 50
C**** SUR LE BORD DU TETRAEDRE, ON A
      RT10=0.0
      RETURN
   50 CONTINUE
C**** A L'INTERIEUR DU TETRAEDRE, ON A
      IF (X.GT.Y .AND. Z.GE.Y) GO TO 200
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 300
  100 CONTINUE
C**** MORCEAU 1.
      RT10=ST10(W,X,Y,Z)
      RETURN
  200 CONTINUE
C**** MORCEAU 2.
      CH=X
      X=Y
      Y=CH
      GO TO 100
  300 CONTINUE
C**** MORCEAU 3.
      CH=X
      X=Z
      Z=CH
      GO TO 100 
      END
C--------------------------------------------------------------
C
      FUNCTION ST13 (WT,XT,YT,ZT)
C
C--------------------------------------------------------------	
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE RT(1,3), RESTREINTE AU
C**** MORCEAU 1, AU POINT (XT,YT,ZT).	WT=1.-XT-YT-ZT.
C**** NOUS SUPPOSONS QUE LE POINT (XT,YT,ZT) NE SOIT PAS SITUE SUR
C**** UNE DES CINQ ARETES DU MORCEAU 1 APPARTENANT AU BORD DU TETRAEDRE. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT 
      YMX=Y-X
      ZMX=Z-X
      IF (YMX.GT.EPS .OR. ZMX.GT.EPS) GO TO 10
C**** SUR LA DROITE X=Y=Z, ON A
      ST13=0.0
      RETURN
   10 CONTINUE
C**** DANS LE RESTE DU MORCEAU 1, ON A
      VAR=1./(1.-3.*X) - 1./(W+Y)
      VAR=-4.*X + 2.*Z + X*ZMX*VAR
      VAR=W*YMX*VAR /( Y*(YMX+ZMX)*(W+YMX)*(W+Y)*(1.-3.*X) )
      VAR=VAR + 2./( (W+X)*(W+ZMX) )
      ST13=2.*X*X*ZMX*ZMX*W*VAR
      RETURN
      END
C------------------------------------------------------------------ 
C
      FUNCTION RT13 (WT,XT,YT,ZT)
C
C------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE PAR MORCEAUX RT(1,3) AU 
C**** POINT (XT,YT,ZT).	       WT=1.-XY-YT-ZT.
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 300
      IF (X.GT.EPS .AND. Y.GT.EPS .AND.
     1	  Z.GT.EPS .AND. W.GT.EPS ) GO TO 50
C**** SUR LE BORD DU TETRAEDRE, ON A
      RT13=0.0
      RETURN
   50 CONTINUE
C**** A L'INTERIEUR DU TETRAEDRE, ON A 
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 200
  100 CONTINUE 
C**** MORCEAU 1. 
      RT13=ST13(W,X,Y,Z)
      RETURN
  200 CONTINUE
C**** MORCEAU 2.
      CH=X
      X=Y
      Y=CH
      GO TO 100
  300 CONTINUE
C**** MORCEAU 3. 
      RT13=0.0 
      RETURN
      END
C-----------------------------------------------------------------
C
      SUBROUTINE UT44 (NFCT, XT,YT,ZT, UT)
C
C-----------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE UT(1), ... , UT(NFCT) AU POINT (XT,YT,ZT) DU TETRAEDRE 
C**** DE REFERENCE.
C**** NFCT PEUT PRENDRE LES VALEURS 16, 28, 32 OU 44.
      DIMENSION UT(44)
      N=NFCT
      X=XT
      Y=YT
      Z=ZT
      W=1.-X-Y-Z
      DO 20 I=1,4
C****    BLOC 1.    FONCTIONS NUMERO 1 A 16.
         IND=4*I - 3
         UT(IND  )=PT1(0,W,X,Y,Z) 
         UT(IND+1)=PT1(1,W,X,Y,Z) 
         UT(IND+2)=PT1(1,W,Y,Z,X) 
         UT(IND+3)=PT1(1,W,Z,X,Y) 
         IF (N.LE.16) GO TO 10
C****    BLOC 2.    FONCTIONS NUMERO 17 A 28.
         IND=3*I+14
         UT(IND  )=WT11(W,X,Y,Z) 
         UT(IND+1)=WT11(W,Y,Z,X) 
         UT(IND+2)=WT11(W,Z,X,Y) 
         IF (N.LE.28) GO TO 10
C****    BLOC 3.    FONCTIONS NUMERO 29 A 32.
         UT(I+28 )=RT10(W,X,Y,Z) 
         IF (N.LE.32) GO TO 10
C****    BLOC 4.    FONCTIONS NUMERO 33 A 44.
         IND=3*I+30
         UT(IND  )=RT13(W,Y,Z,X) 
         UT(IND+1)=RT13(W,Z,X,Y) 
         UT(IND+2)=RT13(W,X,Y,Z)
   10    CH=W
         W=X
         X=Y
         Y=Z
         Z=CH
   20 CONTINUE
      RETURN
      END
C------------------------------------------------------------------
C
      FUNCTION DPT1 (J, TX,TY,TZ, WT,XT,YT,ZT)
C
C------------------------------------------------------------------	
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DU POLYNOME DE DEGRE TROIS PT(1,J) DANS LA DIRECTION 
C**** (TX,TY,TZ) AU POINT (XT,YT,ZT).
C**** J PEUT PRENDRE LES VALEURS 0 OU 1.       WT=1.-XT-YT-ZT.
      W=WT 
      X=XT
      Y=YT
      Z=ZT
      T1=TX 
      T2=TY
      T3=TZ 
      T0=-T1-T2-T3
      IF (J) 10,10,11
   10 P=2.-W - 2.*(X*X+Y*Y+Z*Z+X*Y+Y*Z+Z*X)
      DP=-T0 - 2.*( T1*(2.*X+Y+Z)+T2*(2.*Y+Z+X)+T3*(2.*Z+X+Y) ) 
      GO TO 20
   11 P=X*(W+0.5*(Y+Z))
      DP =T0*X + T1*(W+0.5*(Y+Z)) + (T2+T3)*0.5*X
   20 DPT1=T0*P + W*DP
      RETURN 
      END
C------------------------------------------------------------------
C
      FUNCTION DWT11 (TX,TY,TZ, WT,XT,YT,ZT)
C
C------------------------------------------------------------------	
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DU POLYNOME PAR MORCEAUX WT(1,1) DANS LA DIRECTION 
C**** (TX,TY,TZ) AU POINT (XT,YT,ZT).
C**** WT EST INACTIF.
      X=XT
      Y=YT
      Z=ZT
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 20
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 30
C**** MORCEAU 1.
      DWT11=2.*TX*(X*(X-Y-Z)+Y*Z) + X*(TY*(2.*Z-X)+TZ*(2.*Y-X)) 
      RETURN
C**** MORCEAU 2.
   20 DWT11=Y*(TY*(2.*Z-Y) + TZ*Y)
      RETURN
C**** MORCEAU 3.
   30 DWT11=Z*(TY*Z + TZ*(2.*Y-Z))
      RETURN
      END
C------------------------------------------------------------------ 
C
      FUNCTION DST10 (TX,TY,TZ, WT,XT,YT,ZT)
C 
C------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE RT(1,0), RESTREINTE AU 
C**** MORCEAU 1, DANS LA DIRECTION (TX,TY,TZ) AU POINT (XT,YT,ZT). 
C**** WT=1.-XT-YT-ZT.
C**** NOUS SUPPOSONS QUE LE POINT (XT,YT,ZT) NE SOIT PAS SITUE SUR
C**** UNE DES CINQ ARETES DU MORCEAU 1 APPARTENANT AU BORD DU TETRAEDRE. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT 
      YMX=Y-X
      ZMX=Z-X
      T1=TX 
      T2=TY 
      T3=TZ
      T0=-T1-T2-T3
      IF (YMX.GT.EPS .OR. ZMX.GT.EPS) GO TO 10
C**** SUR LA DROITE X=Y=Z, LA VALEUR ET LA DERIVEE DE C SONT NULLES. 
      V1=0.0
      DV1=0.0 
      GO TO 20
   10 CONTINUE
C**** DANS LE RESTE DU MORCEAU 1, ON A
C**** VALEUR ET DERIVEE DE A (PLACEES DANS V1, DV1).
      P  =X*W
      DP =T1*W + T0*X
      V1 =-P*( 3./(1.-3.*X) + 2./(Y+ZMX) )
      DV1=-3.*(DP - P*( -3.*T1/(1.-3.*X) ))/(1.-3.*X)
     1    -2.*(DP - P*( (T2+T3-T1)/(Y+ZMX) ))/(Y+ZMX)
C**** VALEUR ET DERIVEE DE B (PLACEES DANS V2,DV2).
      P  =6.*X-2.*W-3.+V1
      V2 =-YMX*ZMX*P
      DV2=-((T2-T1)*ZMX+YMX*(T3-T1))*P - YMX*ZMX*(6.*T1-2.*T0+DV1)
C**** VALEUR ET DERIVEE DE C (PLACEES DANS V1, DV1).
      P =ZMX*V2
      DP =(T3-T1)*V2 + ZMX*DV2
      Q =(YMX+ZMX)*Y*(W+YMX)
      V1 =P/Q
      DV1=(DP - P*( (T2+T3-2.*T1)/(YMX+ZMX)+T2/Y+(T0+T2-T1)/(W+YMX) ))/Q 
      P =YMX*V2
      DP =(T2-T1)*V2 + YMX*DV2
      Q  =(YMX+ZMX)*Z*(W+ZMX)
      V1 =V1 + P/Q
      DV1=DV1 + (DP - P*( (T2+T3-2.*T1)/(YMX+ZMX)+T3/Z 
     1                   +(T0+T3-T1)/(W+ZMX) ))/Q 
   20 CONTINUE
C**** VALEUR ET DERIVEE DE D (PLACEES DANS V2, DV2). 
      P =W*( 2.*(YMX+ZMX)+V1 )
      DP =T0*( 2.*(YMX+ZMX)+V1)+ W*( 2.*(T2+T3-2.*T1)+DV1 ) 
      Q =(1.-3.*X)*(Y+ZMX)
      V2 =P/Q
      DV2=(DP - P*( -3.*T1/(1.-3.*X)+(T2+T3-T1)/(Y+ZMX) ))/Q 
C**** DERIVEE DE RT(1,0).
      DST10=X*( (T1*2.*W+X*T0)*(3.+V2) + X*W*DV2 ) 
      RETURN
      END
C------------------------------------------------------------------
C
      FUNCTION DRT10 (TX,TY,TZ, WT,XT,YT,ZT)
C 
C------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE PAR MORCEAUX RT(1,0) DANS 
C**** LA DIRECTION (TX,TY,TZ) AU POINT (XT,YT,ZT).  WT=1.-XT-YT-ZT.
      DATA EPS/0.5E-06/
      W=WT 
      X=XT
      Y=YT
      Z=ZT
      IF (X.GT.EPS .AND. Y.GT.EPS .AND. Z.GT.EPS) GO TO 50
C**** DANS LES FACES 2, 3 ET 4, ON A
      DRT10=0.0
      RETURN
   50 CONTINUE
      T1=TX 
      T2=TY 
      T3=TZ
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 200
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 300
  100 CONTINUE
C**** MORCEAU 1.
      IF (W.GT.EPS) GO TO 150
C**** DANS L'INTERSECTION DU MORCEAU 1 ET DE LA FACE 1, ON A
      DRT10= (-T1-T2-T3)*3.*X*X
      RETURN 
  150 CONTINUE
C**** DANS L'INTERSECTION DU MORCEAU 1 ET DE L'INTERIEUR DU TETRAEDRE, 
C**** ON A
      DRT10=DST10(T1,T2,T3,W,X,Y,Z)
      RETURN 
  200 CONTINUE
C**** MORCEAU 2.
      CH=X
      X=Y
      Y=CH 
      CH=T1 
      T1=T2 
      T2=CH
      GO TO 100
  300 CONTINUE
C**** MORCEAU 3.
      CH=X
      X=Z
      Z=CH 
      CH=T1 
      T1=T3 
      T3=CH
      GO TO 100
      END
C------------------------------------------------------------- 
C
      FUNCTION DST13 (TX,TY,TZ, WT,XT,YT,ZT)
C
C-------------------------------------------------------------	
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE RT(1,3), RESTREINTE AU
C**** MORCEAU 1, DANS LA DIRECTION (TX,TY,TZ) AU POINT (XT,YT,ZT).
C**** WT=1.-XT-YT-ZT.
C**** NOUS SUPPOSONS QUE LE POINT (XT,YT,ZT) NE SOIT PAS SITUE SUR
C**** UNE DES CINQ ARETES DU MORCEAU 1 APPARTENANT AU BORD DU TETRAEDRE. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT
      Z=ZT
      YMX=Y-X 
      ZMX=Z-X
      IF (YMX.GT.EPS .OR. ZMX.GT.EPS) GO TO 10
C**** SUR LA DROITE X=Y=Z, ON A
      DST13=0.0
      RETURN
   10 CONTINUE
C**** DANS LE RESTE DU MORCEAU 1, ON A
      T1=TX 
      T2=TY 
      T3=TZ 
      T0=-T1-T2-T3
C**** VALEUR ET DERIVEE DE E.
      P =X*ZMX
      DP=T1*ZMX + X*(T3-T1)
      Q =1.-3.*X
      E =P/Q
      DE=(DP - E*(-3.*T1))/Q
      Q =W+Y
      E =-4.*X+2.*Z + E - P/Q
      DE=-4.*T1+2.*T3 + DE - ( DP - P*(T0+T2)/Q )/Q
C**** VALEUR ET DERIVEE DE F.
      P =W*YMX*E
      DP=T0*YMX*E + W*(T2-T1)*E + W*YMX*DE
      Q =Y*(YMX+ZMX)*(W+YMX)*(W+Y)*(1.-3.*X)
      F =P/Q
      DF=(DP - P*( T2/Y+(T2+T3-2.*T1)/(YMX+ZMX)+(T0+T2-T1)/(W+YMX)
     1          +(T0+T2)/(W+Y)-3.*T1/(1.-3.*X) ))/Q
      Q =(W+X)*(W+ZMX)
      F =F + 2./Q
      DF=DF - 2.*( (T0+T1)/(W+X)+(T0+T3-T1)/(W+ZMX) )/Q
C**** DERIVEE DE DST13.
      DST13=2.*X*ZMX*( 2.*(T1*ZMX+X*(T3-T1))*W*F + X*ZMX*(T0*F+W*DF) ) 
      RETURN
      END
C------------------------------------------------------------------
C
      FUNCTION DRT13 (TX,TY,TZ, WT,XT,YT,ZT)
C 
C-------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE PAR MORCEAUX RT(1,3) DANS 
C**** LA DIRECTION (TX,TY,TZ) AU POINT (XT,YT,ZT).   WT=1.-XT-YT-ZT. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 300
      IF (X.GT.EPS .AND. Y.GT.EPS .AND. Z.GT.EPS) GO TO 50
C**** DANS LES FACES 2, 3 ET 4, ON A
      DRT13=0.0
      RETURN
   50 CONTINUE
      T1=TX 
      T2=TY 
      T3=TZ
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 200
  100 CONTINUE
C**** MORCEAU 1.
      IF (W.GT.EPS) GO TO 150
C**** DANS L'INTERSECTION DU MORCEAU 1 ET DE LA FACE 1, ON A 
      DRT13=4.*(-T1-T2-T3)*X*(Z-X)
      RETURN
  150 CONTINUE
C**** DANS L'INTERSECTION DU MORCEAU 1 ET DE L'INTERIEUR DU TETRAEDRE, 
C**** ON A
      DRT13=DST13(T1,T2,T3,W,X,Y,Z)
      RETURN
  200 CONTINUE
C**** MORCEAU 2.
      CH=X
      X=Y
      Y=CH
      CH=T1
      T1=T2
      T2=CH
      GO TO 100
  300 CONTINUE
C**** MORCEAU 3. 
      DRT13=0.0 
      RETURN
      END
C-------------------------------------------------------------
C
      SUBROUTINE DUT44 (NFCT, TX,TY,TZ, XT,YT,ZT, DUT)
C 
C-------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE UT(1), ..., UT(NFCT) DANS LA DIRECTION (TX,TY,TZ) 
C**** AU POINT (XT,YT,ZT).
C**** NFCT PEUT PRENDRE LES VALEURS 16, 28, 32 OU 44. 
      DIMENSION DUT(44)
      N=NFCT
      T1=TX
      T2=TY
      T3=TZ
      T0=-T1-T2-T3
      X=XT
      Y=YT
      Z=ZT
      W=1.-X-Y-Z
      DO 20 I=1,4
C****    BLOC 1.   FONCTIONS NUMERO 1 A 16. 
         IND=4*I-3
         DUT(IND  )=DPT1(0,T1,T2,T3,W,X,Y,Z) 
         DUT(IND+1)=DPT1(1,T1,T2,T3,W,X,Y,Z) 
         DUT(IND+2)=DPT1(1,T2,T3,T1,W,Y,Z,X) 
         DUT(IND+3)=DPT1(1,T3,T1,T2,W,Z,X,Y) 
         IF (N.LE.16) GO TO 10
C****    BLOC 2.   FONCTIONS NUMERO 17 A 28. 
         IND=3*I+14
         DUT(IND  )=DWT11(T1,T2,T3,W,X,Y,Z) 
         DUT(IND+1)=DWT11(T2,T3,T1,W,Y,Z,X) 
         DUT(IND+2)=DWT11(T3,T1,T2,W,Z,X,Y) 
         IF (N.LE.28) GO TO 10
C****    BLOC 3.   FONCTIONS NUMERO 29 A 32.
         DUT(I+28 )=DRT10(T1,T2,T3,W,X,Y,Z) 
         IF (N.LE.32) GO TO 10
C****    BLOC 4.   FONCTIONS NUMERO 33 A 44. 
         IND=3*I+30
         DUT(IND  )=DRT13(T2,T3,T1,W,Y,Z,X) 
         DUT(IND+1)=DRT13(T3,T1,T2,W,Z,X,Y) 
         DUT(IND+2)=DRT13(T1,T2,T3,W,X,Y,Z)
   10    CH=W
         W=X 
         X=Y 
         Y=Z
         Z=CH
         CH=T0 
         T0=T1 
         T1=T2 
         T2=T3 
         T3=CH
   20 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------
C
      SUBROUTINE ELEM0 (X1,Y1,Z1, X2,Y2,Z2, X3,Y3,Z3, X4,Y4,Z4) 
C
C-----------------------------------------------------------------	
C**** ELEMENT GENERIQUE ELARGI.
C**** INITIALISATIONS RELATIVES AU TETRAEDRE GENERIQUE.
C**** ENTREES  :  (X1,Y1,Z1), (X2,Y2,Z2), (X3,Y3,Z3), (X4,Y4,Z4) SOMMETS
C****                    DU TETRAEDRE GENERIQUE.
C**** SORTIES  :  TETRAS BLOC COMMUN PAR LEQUEL SE FONT LES SORTIES.
C****             LL     PARTIE LINEAIRE DE L'APPLICATION AFFINE L,
C****                    FACTORISEE SELON LA METHODE DE GAUSS.
C****             IPERM  VECTEUR CONTENANT LA NOUVELLE NUMEROTATION DES
C****                    LIGNES DE LL, CONSECUTIVE AU CHOIX DES PIVOTS.
C****             MX(I), MY(I), MZ(I), VECTEUR NORMAL A LA I-EME FACE DU
C****                    TETRAEDRE, I=1(1)4.
C****             NX(I,J), NY(I,J), NZ(I,J), J=1,2,3, VECTEURS PARAL-
C****                    LELES A LA I-EME FACE ET NORMAUX AUX ARETES DU
C****                    TETRAEDRE, I=1(1)4.
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION AUX(3)
      X(1)=X1 
      Y(1)=Y1 
      Z(1)=Z1 
      X(2)=X2 
      Y(2)=Y2 
      Z(2)=Z2 
      X(3)=X3 
      Y(3)=Y3 
      Z(3)=Z3 
      X(4)=X4 
      Y(4)=Y4 
      Z(4)=Z4
C**** PARTIE LINEAIRE DE L'APPLICATION AFFINE L.
      LL(1,1)=X(2)-X(1)
      LL(2,1)=Y(2)-Y(1)
      LL(3,1)=Z(2)-Z(1)
      LL(1,2)=X(3)-X(1)
      LL(2,2)=Y(3)-Y(1) 
      LL(3,2)=Z(3)-Z(1) 
      LL(1,3)=X(4)-X(1) 
      LL(2,3)=Y(4)-Y(1)
      LL(3,3)=Z(4)-Z(1)
C**** FACTORISATION DE LL.
      CALL DECOMP(3,3,LL,AUX,IPERM,DET)
C**** VOL=ABS(DET)/6. EST LE VOLUME DU TETRAEDRE.
C**** SI LE TETRAEDRE EST DEGENERE (VOL=0.), LA SUBROUTINE DECOMP 
C**** IMPRIME UN MESSAGE D'ERREUR.
      DO 10 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1 
         L=MOD(K,4)+1
C****    INITIALISATION DU VECTEUR (MX(I), MY(I), MZ(I)). 
         MX(I)=(Y(K)-Y(J))*(Z(L)-Z(J)) - (Z(K)-Z(J))*(Y(L)-Y(J)) 
         MY(I)=(Z(K)-Z(J))*(X(L)-X(J)) - (X(K)-X(J))*(Z(L)-Z(J)) 
         MZ(I)=(X(K)-X(J))*(Y(L)-Y(J)) - (Y(K)-Y(J))*(X(L)-X(J)) 
         S=SQRT(MX(I)**2 + MY(I)**2 + MZ(I)**2)
         MX(I)=MX(I)/S
         MY(I)=MY(I)/S
         MZ(I)=MZ(I)/S
C****    VECTEURS   (NX(I,J), NY(I,J), NZ(I,J)),  J=1,2,3. 
         NX(I,1)=MY(I)*(Z(L)-Z(K)) - MZ(I)*(Y(L)-Y(K)) 
         NY(I,1)=MZ(I)*(X(L)-X(K)) - MX(I)*(Z(L)-Z(K)) 
         NZ(I,1)=MX(I)*(Y(L)-Y(K)) - MY(I)*(X(L)-X(K)) 
         S=SQRT(NX(I,1)**2+NY(I,1)**2+NZ(I,1)**2) 
         NX(I,1)=NX(I,1)/S
         NY(I,1)=NY(I,1)/S
         NZ(I,1)=NZ(I,1)/S
         NX(I,2)=MY(I)*(Z(J)-Z(L)) - MZ(I)*(Y(J)-Y(L)) 
         NY(I,2)=MZ(I)*(X(J)-X(L)) - MX(I)*(Z(J)-Z(L)) 
         NZ(I,2)=MX(I)*(Y(J)-Y(L)) - MY(I)*(X(J)-X(L)) 
         S=SQRT(NX(I,2)**2+NY(I,2)**2+NZ(I,2)**2)
         NX(I,2)=NX(I,2)/S
         NY(I,2)=NY(I,2)/S
         NZ(I,2)=NZ(I,2)/S
         NX(I,3)=MY(I)*(Z(K)-Z(J)) - MZ(I)*(Y(K)-Y(J)) 
         NY(I,3)=MZ(I)*(X(K)-X(J)) - MX(I)*(Z(K)-Z(J)) 
         NZ(I,3)=MX(I)*(Y(K)-Y(J)) - MY(I)*(X(K)-X(J)) 
         S=SQRT(NX(I,3)**2+NY(I,3)**2+NZ(I,3)**2)
         NX(I,3)=NX(I,3)/S
         NY(I,3)=NY(I,3)/S
         NZ(I,3)=NZ(I,3)/S
   10 CONTINUE
      RETURN
      END
C----------------------------------------------------------------
C
      SUBROUTINE U44 (NFCT, PX,PY,PZ, U)
C 
C----------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT 
C****                    INITIALISE A L'AIDE DE LA SUBROUTINE ELEM0.
C****             NFCT   PEUT PRENDRE LES VALEURS 16, 28, 32 OU 44.
C****             (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE   :  U(1),	..., U(NFCT) VALEURS DES NFCT PREMIERES FONC-
C****                    TIONS DE BASE AU POINT (PX, PY, PZ). 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION U(44), GEN(3), REF(3)
      DATA NERR/0/, EPS/-1.E-03/
      GEN(1)=PX-X(1) 
      GEN(2)=PY-Y(1) 
      GEN(3)=PZ-Z(1)
      CALL SOLVE(3,3,LL,GEN,REF,IPERM)
      XT=REF(1)
      YT=REF(2)
      ZT=REF(3)
      IF (XT.GT.EPS .AND. YT.GT.EPS .AND. ZT.GT.EPS .AND.
     1    1.-XT-YT-ZT.GT.EPS) GO TO 20
C**** LE POINT (PX,PY,PZ) EST A L'EXTERIEUR DU TETRAEDRE GENERIQUE.
      IF (NERR.LT.10) WRITE(UNIT=*, FMT=*)
     1                ' *** ERROR MESSAGE CALLED FROM U44 :',
     2                ' THE POINT IS OUTSIDE THE ELEMENT.'
      NERR=NERR+1
   20 CONTINUE
C**** VALEURS DES FONCTIONS DE BASE.
      CALL UT44(NFCT,XT,YT,ZT,U)
      RETURN
      END
C------------------------------------------------------------
C
      SUBROUTINE DU44 (NFCT, VX,VY,VZ, PX,PY,PZ, DU) 
C
C------------------------------------------------------------ 	
C**** BASE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT 
C****                    INITIALISE A L'AIDE DE LA SUBROUTINE ELEM0.
C****             NFCT   PEUT PRENDRE LES VALEURS 16, 28, 32 OU 44.
C****	          (VX, VY, VZ) UN VECTEUR.
C****             (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE   :  DU(1), ..., DU(NFCT) DERIVEES DES NFCT PREMIERES
C****                    FONCTIONS DE BASE DANS LA DIRECTION (VX,VY,VZ)
C****                    AU POINT (PX, PY, PZ). 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION DU(44), GEN(3), REF(3)
      DATA NERR/0/, EPS/-1.E-03/
      GEN(1)=VX 
      GEN(2)=VY 
      GEN(3)=VZ
      CALL SOLVE(3,3,LL,GEN,REF,IPERM)
      TX=REF(1) 
      TY=REF(2) 
      TZ=REF(3) 
      GEN(1)=PX-X(1)
      GEN(2)=PY-Y(1)
      GEN(3)=PZ-Z(1)
      CALL SOLVE(3,3,LL,GEN,REF,IPERM)
      XT=REF(1) 
      YT=REF(2) 
      ZT=REF(3)
      IF (XT.GT.EPS .AND. YT.GT.EPS .AND. ZT.GT.EPS .AND. 
     1    1.-XT-YT-ZT.GT.EPS) GO TO 20
C**** LE POINT (PX,PY,PZ) EST A L'EXTERIEUR DU TETRAEDRE GENERIQUE. 
      IF (NERR.LT.10) WRITE(UNIT=*, FMT=*)
     1                ' *** ERROR MESSAGE CALLED FROM DU44 :',
     2                ' THE POINT IS OUTSIDE THE ELEMENT.'
      NERR=NERR+1
   20 CONTINUE
C**** DERIVEES DES FONCTIONS DE BASE.
      CALL DUT44(NFCT,TX,TY,TZ,XT,YT,ZT,DU)
      RETURN
      END
C-------------------------------------------------------------
C
      SUBROUTINE ELEM1
C
C-------------------------------------------------------------
C**** MATRICE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** INITIALISATION DE L'INVERSE DES BLOCS DIAGONAUX DE LA MATRICE Q. 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DO 10 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1 
         L=MOD(K,4)+1
C****    INITIALISATION DE D1, D2, D3 ET D4.
         D(I,1,1)=X(J)-X(I)
         D(I,1,2)=Y(J)-Y(I)
         D(I,1,3)=Z(J)-Z(I) 
         D(I,2,1)=X(K)-X(I) 
         D(I,2,2)=Y(K)-Y(I) 
         D(I,2,3)=Z(K)-Z(I) 
         D(I,3,1)=X(L)-X(I) 
         D(I,3,2)=Y(L)-Y(I) 
         D(I,3,3)=Z(L)-Z(I) 
C****    INITIALISATION DE D5. 
         IND=3*I-2
         D5(IND  )=NX(I,1)*(2.*X(J)-X(K)-X(L))
     1            +NY(I,1)*(2.*Y(J)-Y(K)-Y(L))
     2            +NZ(I,1)*(2.*Z(J)-Z(K)-Z(L))
         D5(IND+1)=NX(I,2)*(2.*X(K)-X(L)-X(J))
     1            +NY(I,2)*(2.*Y(K)-Y(L)-Y(J))
     2            +NZ(I,2)*(2.*Z(K)-Z(L)-Z(J))
         D5(IND+2)=NX(I,3)*(2.*X(L)-X(J)-X(K))
     1            +NY(I,3)*(2.*Y(L)-Y(J)-Y(K))
     2            +NZ(I,3)*(2.*Z(L)-Z(J)-Z(K))
C****    INITIALISATION DE D6.
         D6(I)=MX(I)*(3.*X(I)-X(J)-X(K)-X(L))
     1        +MY(I)*(3.*Y(I)-Y(J)-Y(K)-Y(L))
     2        +MZ(I)*(3.*Z(I)-Z(J)-Z(K)-Z(L))
C****	 INITIALISATION DE D7.
         IND=3*I-2
         D7(IND  )=D6(I)
         D7(IND+1)=D6(I) 
         D7(IND+2)=D6(I) 
   10 CONTINUE
      RETURN
      END
C-------------------------------------------------------------
C
      SUBROUTINE ELEM2
C
C-------------------------------------------------------------	
C**** MATRICE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI. 
C**** INITIALISATION DES BLOCS SOUS-DIAGONAUX DE LA MATRICE Q. 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	     NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION DU(44)
      DO 80 I=1,4
         J=MOD(I,4)+1
         K=MOD(J,4)+1
         L=MOD(K,4)+1
C****    INITIALISATION DU BLOC Q1.
         PX=0.5*(X(K)+X(L)) 
         PY=0.5*(Y(K)+Y(L)) 
         PZ=0.5*(Z(K)+Z(L))
         CALL DU44(16,NX(I,1),NY(I,1),NZ(I,1),PX,PY,PZ,DU)
         IND=3*I-2
         DO 10 M=1,16
            Q1(IND,M)=DU(M)
   10    CONTINUE
         PX=0.5*(X(J)+X(L)) 
         PY=0.5*(Y(J)+Y(L)) 
         PZ=0.5*(Z(J)+Z(L))
         CALL DU44(16,NX(I,2),NY(I,2),NZ(I,2),PX,PY,PZ,DU)
         IND=IND+1
         DO 20 M=1,16
            Q1(IND,M)=DU(M)
   20    CONTINUE
         PX=0.5*(X(K)+X(J)) 
         PY=0.5*(Y(K)+Y(J)) 
         PZ=0.5*(Z(K)+Z(J))
         CALL DU44(16,NX(I,3),NY(I,3),NZ(I,3),PX,PY,PZ,DU)
         IND=IND+1
         DO 30 M=1,16
            Q1(IND,M)=DU(M)
   30    CONTINUE
C****    INITIALISATION DU BLOC Q2. 
         PX=(X(J)+X(K)+X(L))/3. 
         PY=(Y(J)+Y(K)+Y(L))/3. 
         PZ=(Z(J)+Z(K)+Z(L))/3.
         CALL DU44(28,MX(I),MY(I),MZ(I),PX,PY,PZ,DU)
         DO 40 M=1,28
            Q2(I,M)=DU(M)
   40    CONTINUE
C****    INITIALISATION DU BLOC Q3. 
         PX=(4.*X(J)+X(K)+X(L))/6. 
         PY=(4.*Y(J)+Y(K)+Y(L))/6. 
         PZ=(4.*Z(J)+Z(K)+Z(L))/6.
         CALL DU44(32,MX(I),MY(I),MZ(I),PX,PY,PZ,DU)
         IND=3*I-2
         DO 50 M=1,32
            Q3(IND,M)=DU(M)
   50    CONTINUE
         PX=(X(J)+4.*X(K)+X(L))/6. 
         PY=(Y(J)+4.*Y(K)+Y(L))/6.
         PZ=(Z(J)+4.*Z(K)+Z(L))/6.
         CALL DU44(32,MX(I),MY(I),MZ(I),PX,PY,PZ,DU) 
         IND=IND+1
         DO 60 M=1,32
            Q3(IND,M)=DU(M) 
   60    CONTINUE
         PX=(X(J)+X(K)+4.*X(L))/6.
         PY=(Y(J)+Y(K)+4.*Y(L))/6. 
         PZ=(Z(J)+Z(K)+4.*Z(L))/6.
         CALL DU44(32,MX(I),MY(I),MZ(I),PX,PY,PZ,DU) 
         IND=IND+1
         DO 70 M=1,32
            Q3(IND,M)=DU(M) 
   70    CONTINUE
   80 CONTINUE
      RETURN
      END
C--------------------------------------------------------------
C
      SUBROUTINE FE2P44 (C, S)
C
C--------------------------------------------------------------
C**** ELEMENT GENERIQUE ELARGI.
C**** RESOLUTION DU PROBLEME D'INTERPOLATION.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT
C****                    INITIALISE A L'AIDE DES SUBROUTINES ELEM0,
C****                    ELEM1 ET ELEM2.
C****             C(1), ..., C(44) VALEURS DES PARAMETRES DE L'ELEMENT.
C**** SORTIE   :  S(1), ..., S(44) COORDONNEES DE LA FONCTION D'INTER-
C****                    POLATION DANS LA BASE D'INTERPOLATION DE
C****                    L'ELEMENT GENERIQUE ELARGI. 
C**** METHODE  :  RESOLUTION DU SYSTEME LINEAIRE Q*S=C. 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION C(44), S(44)
C**** CALCUL DE S(1),	..., S(16).
      DO 10 I=1,4
         K=4*I-3
         S(K)=C(K)
         S(K+1)=D(I,1,1)*C(K+1)+D(I,1,2)*C(K+2)+D(I,1,3)*C(K+3) 
         S(K+2)=D(I,2,1)*C(K+1)+D(I,2,2)*C(K+2)+D(I,2,3)*C(K+3) 
         S(K+3)=D(I,3,1)*C(K+1)+D(I,3,2)*C(K+2)+D(I,3,3)*C(K+3)
   10 CONTINUE
C**** CALCUL DE S(17), ..., S(28).
      DO 30 I=1,12
         SI=0.0
         DO 20 J=1,16
            SI=SI+Q1(I,J)*S(J)
   20    CONTINUE 
         K=I+16 
         S(K)=D5(I)*(C(K)-SI)
   30 CONTINUE
C**** CALCUL DE S(29), ..., S(32).
      DO 50 I=1,4
      SI=0.0
         DO 40 J=1,28 
            SI=SI+Q2(I,J)*S(J)
   40    CONTINUE
         K=28+I
         S(K)=D6(I)*(C(K)-SI)
   50 CONTINUE
C**** CALCUL DE S(33),	..., S(44).
      DO 70 I=1,12
         SI=0.0
         DO 60 J=1,32 
            SI=SI+Q3(I,J)*S(J)
   60    CONTINUE
         K=I+32 
         S(K)=D7(I)*(C(K)-SI) 
   70 CONTINUE
      RETURN
      END
C--------------------------------------------------------
C
      FUNCTION D1NOD (VX,VY,VZ, M, C)
C
C--------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX,
C**** A 16 PARAMETRES.
C**** DONNEES  :  M      PEUT PRENDRE LES VALEURS 1, 2, 3 OU 4.
C****             C(4*M-2), C(4*M-1) ET C(4*M) DERIVEES PARTIELLES DX,
C****                    DY ET DZ AU SOMMET A(M).
C****             (VX, VY, VZ) UN VECTEUR.
C**** SORTIE   :  D1NOD  DERIVEE DANS LA DIRECTION (VX, VY, VZ)
C****                    AU SOMMET A(M).
      DIMENSION C(44)
      IND=4*M-2
      D1NOD=VX*C(IND) + VY*C(IND+1) + VZ*C(IND+2)
      RETURN
      END
C--------------------------------------------------------
C
      SUBROUTINE FE2P16 (C, S)
C
C--------------------------------------------------------	
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX, 
C**** A 16 PARAMETRES.
C**** RESOLUTION DU PROBLEME D'INTERPOLATION.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT
C****                    INITIALISE A L'AIDE DES SUBROUTINES ELEM0,
C****                    ELEM1 ET ELEM2.
C****             C(1), ..., C(16) VALEURS DES PARAMETRES DE L'ELEMENT.
C**** SORTIE   :  S(1), ..., S(44) COORDONNEES DE LA FONCTION D'INTER-
C****                    POLATION DANS LA BASE D'INTERPOLATION DE
C****                    L'ELEMENT GENERIQUE ELARGI.
C**** LES MEMOIRES C(17), ..., C(44) SONT DES PLACES DE TRAVAIL. 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION C(44), S(44)
      DO 10 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1
         L=MOD(K,4)+1
C****    CALCUL DES PARAMETRES C(17) A C(28).
         IND=3*I+14
         C(IND  )=0.5*( D1NOD(NX(I,1),NY(I,1),NZ(I,1),K,C)
     1                 +D1NOD(NX(I,1),NY(I,1),NZ(I,1),L,C) ) 
         C(IND+1)=0.5*( D1NOD(NX(I,2),NY(I,2),NZ(I,2),L,C)
     1                 +D1NOD(NX(I,2),NY(I,2),NZ(I,2),J,C) )
         C(IND+2)=0.5*( D1NOD(NX(I,3),NY(I,3),NZ(I,3),J,C)
     1                 +D1NOD(NX(I,3),NY(I,3),NZ(I,3),K,C) ) 
C****    CALCUL DES PARAMETRES C(29) A C(32). 
         DJ=D1NOD(MX(I),MY(I),MZ(I),J,C)
         DK=D1NOD(MX(I),MY(I),MZ(I),K,C)
         DL=D1NOD(MX(I),MY(I),MZ(I),L,C)
         C(I+28)=(DJ+DK+DL)/3.
C****    CALCUL DES PARAMETRES C(33) A C(44).
         IND=3*I+30
         C(IND  )=(4.*DJ+DK+DL)/6.
         C(IND+1)=(DJ+4.*DK+DL)/6.
         C(IND+2)=(DJ+DK+4.*DL)/6.
   10 CONTINUE
C**** CALCUL DES COORDONNEES S(1), ... , S(44).
      CALL FE2P44(C,S)
      RETURN 
      END
C------------------------------------------------------------ 
C
      FUNCTION FINT (PX,PY,PZ, S)
C 
C------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX,
C**** A 16 PARAMETRES.
C**** DONNEES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT 
C****	                 INITIALISE A L'AIDE DE LA SUBROUTINE ELEM0.
C****             S(1), ..., S(44) COORDONNEES DE LA FONCTION D'INTER-
C****                    POLATION DANS LA BASE D'INTERPOLATION DE
C****                    L'ELEMENT GENERIQUE ELARGI.
C****             (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE   :  FINT	 VALEUR DE LA FONCTION D'INTERPOLATION DE
C****                    COORDONNEES S AU POINT (PX, PY, PZ).
      DIMENSION U(44), S(44)
      CALL U44(44,PX,PY,PZ,U)
      V=0.0
      DO 10 I=1,44
         V=V+S(I)*U(I)
   10 CONTINUE 
      FINT=V 
      RETURN 
      END
C--------------------------------------------------------------
C
      FUNCTION DFINT (VX,VY,VZ, PX,PY,PZ, S)
C 
C--------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX,
C**** A 16 PARAMETRES.
C**** DONNEES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT 
C****                    INITIALISE A L'AIDE DE LA SUBROUTINE ELEM0.
C****             S(1), ..., S(44) COORDONNEES DE LA FONCTION D'INTER-
C****                    POLATION DANS LA BASE D'INTERPOLATION DE
C****                    L'ELEMENT GENERIQUE ELARGI.
C****             (VX, VY, VZ) UN VECTEUR.
C****             (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE   :  DFINT DERIVEE DE LA FONCTION D'INTERPOLATION DE
C****                   COORDONNEES S DANS LA DIRECTION (VX, VY, VZ)
C****                   AU POINT (PX,PY,PZ).
      DIMENSION DU(44), S(44)
      CALL DU44(44,VX,VY,VZ,PX,PY,PZ,DU)
      DV=0.0
      DO 10 I=1,44
         DV=DV+S(I)*DU(I)
   10 CONTINUE 
      DFINT=DV 
      RETURN
      END
C-------------------------------------------------------------
C
      SUBROUTINE ELEM3
C
C-------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX,
C**** A 16 PARAMETRES.
C**** BASE D'HERMITE.   INITIALISATION DE LA MATRICE T.
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION C(44), S(44) 
      DATA C/44*0.0/
      DO 20 I=1,16
C****    CALCUL DES COORDONNEES T(1,I), ..., T(44,I) DU I-EME ELEMENT
C****	 DE LA BASE D'HERMITE PAR RAPPORT A LA BASE D'INTERPOLATION
C****	 DE L'ELEMENT GENERIQUE ELARGI.
         C(I)=1.0
         CALL FE2P16(C,S)
         DO 10 J=1,44
            T(J,I)=S(J)
   10    CONTINUE 
         C(I)=0.0 
   20 CONTINUE
      RETURN
      END
C-------------------------------------------------------------- 
C
      SUBROUTINE HERM (PX,PY,PZ, V)
C
C--------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX, 
C**** A 16 PARAMETRES.
C**** BASE D'HERMITE.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT
C****                    INITIALISE A L'AIDE DES SUBROUTINES ELEM0,
C****                    ELEM1, ELEM2 ET ELEM3.
C****             (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE   :  V(1), ..., V(16) VALEURS DES 16 FONCTIONS DE LA BASE
C****                    D'HERMITE AU POINT (PX, PY, PZ).
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3) 
      DIMENSION V(16), U(44)
      CALL U44(44,PX,PY,PZ, U)
      DO 20 I=1,16
         SI=0.0
         DO 10 J=1,44
            SI=SI + T(J,I)*U(J)
   10    CONTINUE
         V(I)=SI
   20 CONTINUE
      RETURN
      END
C----------------------------------------------------------
C
      SUBROUTINE DHERM (VX,VY,VZ, PX,PY,PZ, DV)
C
C----------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX, 
C**** A 16 PARAMETRES.
C**** BASE D'HERMITE.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT 
C****                    INITIALISE A L'AIDE DES SUBROUTINES ELEM0, 
C****                    ELEM1, ELEM2 ET ELEM3.
C****             (VX, VY, VZ) UN VECTEUR.
C****             (PX, PY,PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE   :  DV(1), ..., DV(16) DERIVEES DES 16 FONCTIONS DE LA
C****                    BASE D'HERMITE DANS LA DIRECTION (VX, VY, VZ)
C****                    AU POINT (PX, PY, PZ). 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3) 
      DIMENSION DV(16), DU(44)
      CALL DU44(44,VX,VY,VZ,PX,PY,PZ,DU) 
      DO 20 I=1,16
         SI=0.0
         DO 10 J=1,44
            SI=SI + T(J,I)*DU(J)
   10    CONTINUE 
         DV(I)=SI 
   20 CONTINUE
      RETURN
      END
