C------------------------------------------------------
C
      PROGRAM VARII
C
C**** PROGRAMME PRINCIPAL.
C**** TEST DE L'ELEMENT FINI DE DEGRE 2 A 16 PARAMETRES. 
C**** VARIANTE II.
C
C**** INITIALISATIONS.
      DIMENSION X(4), Y(4), Z(4)
C
      DO 30 NEL=1,2
C****    CHOIX D'UN TETRAEDRE, AU HASARD.
         WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*) ' *** NOUVEL ELEMENT : X, Y, Z'
         DO 20 I=1,4
            CALL RANDOM(YFL) 
            X(I)=4.*YFL-2. 
            CALL RANDOM(YFL) 
            Y(I)=4.*YFL-2. 
            CALL RANDOM(YFL) 
            Z(I)=4.*YFL-2.
   20    CONTINUE
         WRITE(UNIT=*, FMT=110)  (X(I), Y(I), Z(I), I=1,4)
C****    INITIALISATIONS RELATIVES AU TETRAEDRE CHOISI.
         CALL ELEM0 (X(1),Y(1),Z(1), X(2),Y(2),Z(2), 
     1               X(3),Y(3),Z(3), X(4),Y(4),Z(4))
         CALL ELEM1
         CALL ELEM2
         CALL ELEM3
C****    VERIFICATIONS DE L'ELEMENT.
         WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*)
     1   '---*** THESE P. 119, TEST 7.4-1 ***---'
         WRITE(UNIT=*, FMT=*) 'CONTINUITE'
         WRITE(UNIT=*, FMT=120)
         CALL CONT
         WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*)
     1   '---*** THESE P. 119-120, TEST 7.4-2 ***---'
         WRITE(UNIT=*, FMT=*) 'BASE DE LAGRANGE (HERMITE)'
         WRITE(UNIT=*, FMT=120)
         CALL LAGR
         WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*)
     1   '---*** THESE P. 120, TEST 7.4-3 ***---'
         WRITE(UNIT=*, FMT=*) 'POLYNOMES DE DEGRE DEUX'
         WRITE(UNIT=*, FMT=120)
         CALL POL2
         WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*)
     1   '---*** THESE P. 120, TEST 7.4-4 ***---'
         WRITE(UNIT=*, FMT=*) 'DERIVEE NORMALE DE DEGRE .LE. 1'
         WRITE(UNIT=*, FMT=120)
         CALL NORLIN
   30 CONTINUE
       WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*)
     1   '---*** THESE P. 121, TEST 7.4-5 ***---'
         WRITE(UNIT=*, FMT=*)
     1   'MOSAIQUE DE DEUX TETRAEDRES (RACCORDEMENT)'
         WRITE(UNIT=*, FMT=120)
      CALL MOSAIQ
       WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*)
     1   '---*** THESE P. 121, TEST 7.4-6 ***---'
         WRITE(UNIT=*, FMT=*) 'PERMUTATION DES SOMMETS'
         WRITE(UNIT=*, FMT=120)
      CALL PERM (X, Y, Z) 
      STOP
  110 FORMAT (3F18.12) 
  120 FORMAT (/)
      END
C-------------------------------------------------------------
C**** TEST.
      SUBROUTINE CONT
C**** ON VERIFIE SI LES FONCTIONS FINT ET DFINT SONT CONTINUES 
C**** AU VOISINAGE DU BORD.
C**** PASSAGE DE L'EXPRESSION LIMITE SUR LE BORD A L'EXPRESSION 
C**** GENERALE A L'INTERIEUR.
C
      COMMON /TETRAS/ X(4),Y(4),Z(4), TETR(745),IPERM(3)
      DIMENSION S(28)
      DATA S/16*0.0, 12*1.0/
      DATA EPS, RN, M /0.5E-06, 1., 20/ 
      PX=0.2*X(2)+0.3*X(3)+0.5*X(4) 
      PY=0.2*Y(2)+0.3*Y(3)+0.5*Y(4) 
      PZ=0.2*Z(2)+0.3*Z(3)+0.5*Z(4)
      H=2.*RN*EPS/M
      AN=(1.-RN)*EPS
      MP1=M+1
      WRITE(UNIT=*,FMT=*)' *** CONTINU : T, V0, V1, V2, V3'
      DO 30 I=1,MP1
         T=AN+(I-1)*H
         XT=(1.-T)*PX+T*X(1) 
         YT=(1.-T)*PY+T*Y(1) 
         ZT=(1.-T)*PZ+T*Z(1) 
         V0=FINT(XT,YT,ZT,S) 
         V1=DFINT(1.,0.,0.,XT,YT,ZT,S) 
         V2=DFINT(0.,1.,0.,XT,YT,ZT,S) 
         V3=DFINT(0.,0.,1.,XT,YT,ZT,S)
         WRITE(UNIT=*, FMT=20)  T,V0,V1,V2,V3 
   20    FORMAT(E8.2, 4E18.9) 
   30 CONTINUE
      RETURN
      END
C----------------------------------------------------
C**** TEST.
      SUBROUTINE LAGR
C**** NOUS VERIFIONS SI LA BASE D'HERMITE POSSEDE LES PARAMETRES VOULUS.
C
      COMMON /TETRAS/ X(4),Y(4),Z(4), TETR(745), IPERM(3)
      DIMENSION E(16) 
      WRITE(UNIT=*, FMT=*) ' *** LAGR : ERR REL SUR HERM ET DHERM'
      DO 10 I=1,4
C****    J = 0.
         IND=4*I-3
         CALL HERM(X(I), Y(I), Z(I), E)
         E(IND)=E(IND)-1.0
         WRITE(UNIT=*, FMT=110)   (E(K), K=1,16)
C****    J = 1.
         IND=IND+1
         CALL DHERM(1., 0., 0., X(I), Y(I), Z(I), E)
         E(IND)=E(IND)-1.0
         WRITE(UNIT=*, FMT=110)  (E(K), K=1,16)
C****	 J = 2. 
         IND=IND+1
         CALL DHERM(0., 1., 0., X(I), Y(I), Z(I), E)
         E(IND)=E(IND)-1.0
         WRITE(UNIT=*, FMT=110) (E(K), K=1,16)
C****	 J = 3. 
         IND=IND+1
         CALL DHERM(0., 0., 1., X(I), Y(I), Z(I), E)
         E(IND)=E(IND)-1.0
         WRITE(UNIT=*, FMT=110)  (E(K), K=1,16)
         WRITE(UNIT=*, FMT=120)
   10 CONTINUE
      RETURN
  110 FORMAT(16E10.1)
  120 FORMAT(/) 
      END
C----------------------------------------------------------------
C**** TEST.
      SUBROUTINE POL2
C**** NOUS VERIFIONS SI TOUS LES POLYNOMES DE DEGRE DEUX APPARTIENNENT 
C**** A L'ESPACE ELEMENT FINI.
C
      COMMON /TETRAS/ X(4),Y(4),Z(4), TETR(745), IPERM(3)
      DIMENSION C(28), S(28), Q(10), E(20) 
      DATA Q/10*0.0/, H/1.0/
      WRITE(UNIT=*, FMT=*) ' *** POL2 : ERR SUR FINT ET DFINT'
      DO 60 IPI=1,10
C****    CHOIX D'UN POLYNOME DE DEGRE DEUX. 
         Q(IPI)=1.0
C****    CALCUL DES PARAMETRES DE CE POLYNOME.
         DO 20 I=1,4
            IND=4*I-3
            C(IND  )=P(Q,X(I),Y(I),Z(I))
            C(IND+1)=0.5*( P(Q,X(I)+H,Y(I),Z(I))
     1                    -P(Q,X(I)-H,Y(I),Z(I)) )/H
            C(IND+2)=0.5*( P(Q,X(I),Y(I)+H,Z(I))
     1                    -P(Q,X(I),Y(I)-H,Z(I)) )/H
            C(IND+3)=0.5*( P(Q,X(I),Y(I),Z(I)+H)
     1                    -P(Q,X(I),Y(I),Z(I)-H) )/H
   20    CONTINUE
C****    COORDONNEES DE LA FONCTION D'INTERPOLATION CORRESPONDANTE. 
         CALL FE2P16(C, S)
C****    COMPARAISON EN CERTAINS POINTS.
   40       CALL RANDOM(B2)
            CALL RANDOM(B3)
            CALL RANDOM(B4)
            B1=1.-B2-B3-B4
            IF (B1.LT.0.) GO TO 40
         PX=B1*X(1)+B2*X(2)+B3*X(3)+B4*X(4) 
         PY=B1*Y(1)+B2*Y(2)+B3*Y(3)+B4*Y(4) 
         PZ=B1*Z(1)+B2*Z(2)+B3*Z(3)+B4*Z(4)
         IND=2*IPI-1 
         VALP=P(Q,PX,PY,PZ)
         E(IND)=FINT(PX,PY,PZ,S) - VALP
         IND=IND+1
         CALL RANDOM(B1)
         B1=4.*B1-2.
         CALL RANDOM(B2)
         B2=4.*B2-2.
         CALL RANDOM(B3)
         B3=4.*B3-2.
         VALP=0.5*( P(Q,PX+B1*H,PY+B2*H,PZ+B3*H)
     1             -P(Q,PX-B1*H,PY-B2*H,PZ-B3*H) )/H
         E(IND)=DFINT(B1,B2,B3,PX,PY,PZ,S) - VALP
   60 CONTINUE
      WRITE(UNIT=*, FMT=110)  (E(K), K=1,20)
      RETURN
  110 FORMAT(20E10.1)
      END
C----------------------------------------------------------------
C**** TEST.
      FUNCTION P(Q,PX,PY,PZ)
C**** VALEUR DU POLYNOME DE DEGRE 2, DE COEFFICIENTS Q(1), ..., Q(10) 
C**** AU POINT (PX,PY,PZ).
C
      DIMENSION Q(10)
      X=PX 
      Y=PY 
      Z=PZ
      P= Q(1) + Q(2)*X + Q(3)*Y + Q(4)*Z
     1   + Q(5)*X*X + Q(6)*Y*Y + Q(7)*Z*Z
     2   + Q(8)*Y*Z + Q(9)*X*Z + Q(10)*X*Y 
      RETURN
      END
C----------------------------------------------------
C**** TEST.
      SUBROUTINE NORLIN
C**** ON VERIFIE SI LES DERIVEES NORMALES SONT LINEAIRES
C**** DANS LES FACES DU TETRAEDRE.
C 
      REAL MIX,MIY,MIZ
      COMMON /TETRAS/ X(4),Y(4),Z(4), TETR(745),IPERM(3)
      DIMENSION C(28),S(28),E(16)
      WRITE(UNIT=*, FMT=*) ' *** NORLIN : ERR REL SUR DFINT'
C**** CHOIX DES PARAMETRES DE L'ELEMENT FINI. 
      DO 10 I=1,16
         CALL RANDOM(YFL)
         C(I)=4.*YFL-2.
   10 CONTINUE
C**** CALCUL DES COORDONNEES DE LA FONCTION D'INTERPOLATION.
      CALL FE2P16(C,S)
C**** VERIFICATIONS.
      DO 40 I=1,4
         J=MOD(I,4)+1
         K=MOD(J,4)+1
         L=MOD(K,4)+1 
         MIX=(Y(K)-Y(J))*(Z(L)-Z(J))-(Z(K)-Z(J))*(Y(L)-Y(J)) 
         MIY=(Z(K)-Z(J))*(X(L)-X(J))-(X(K)-X(J))*(Z(L)-Z(J)) 
         MIZ=(X(K)-X(J))*(Y(L)-Y(J))-(Y(K)-Y(J))*(X(L)-X(J)) 
         DJ=D1NOD(MIX,MIY,MIZ,J,C)
         DK=D1NOD(MIX,MIY,MIZ,K,C) 
         DL=D1NOD(MIX,MIY,MIZ,L,C) 
         DO 30 NU=1,4
   20       CALL RANDOM(YFL)
            B2=0.667*YFL
            CALL RANDOM(YFL)
            B3=0.667*YFL
            B1=1.-B2-B3
            IF (B1.LT.0.) GO TO 20
C****       INTERPOLATION LINEAIRE DE LA DERIVEE NORMALE DANS LA FACE. 
            DERLIN=B1*DJ+B2*DK+B3*DL
C****       DERIVEE DE LA FONCTION D'INTERPOLATION. 
            PX=B1*X(J)+B2*X(K)+B3*X(L) 
            PY=B1*Y(J)+B2*Y(K)+B3*Y(L) 
            PZ=B1*Z(J)+B2*Z(K)+B3*Z(L) 
            DER=DFINT(MIX,MIY,MIZ,PX,PY,PZ,S)
C****       COMPARAISON.
            IND=4*(I-1)+NU
            E(IND)=(DER-DERLIN)/DERLIN
   30    CONTINUE
   40 CONTINUE
      WRITE(UNIT=*, FMT=50)  (E(I), I=1,16)
   50 FORMAT(16E10.1) 
      RETURN
      END
C------------------------------------------------------------
C
      SUBROUTINE MOSAIQ
C
      DIMENSION PX(5),PY(5),PZ(5), GX(5),GY(5),GZ(5), S(28),C(28), E(20) 
C**** CHOIX DU PLAN DE RACCORDEMENT.
      DO 10 M=2,4
         CALL RANDOM (YFL) 
         PX(M)=4.*YFL - 2. 
         CALL RANDOM (YFL) 
         PY(M)=4.*YFL - 2. 
         CALL RANDOM (YFL) 
         PZ(M)=4.*YFL - 2.
   10 CONTINUE
C**** VECTEUR NORMAL AU PLAN DE RACCORDEMENT.
      RX=(PY(3)-PY(2))*(PZ(4)-PZ(2)) - (PZ(3)-PZ(2))*(PY(4)-PY(2))
      RY=(PZ(3)-PZ(2))*(PX(4)-PX(2)) - (PX(3)-PX(2))*(PZ(4)-PZ(2))
      RZ=(PX(3)-PX(2))*(PY(4)-PY(2)) - (PY(3)-PY(2))*(PX(4)-PX(2))
C**** CHOIX DE POINTS DANS LE PLAN DE RACCORDEMENT.
      DO 30 M=1,5
   20    CALL RANDOM (B2) 
         CALL RANDOM (B3) 
         B4=1.-B2-B3
         IF (B4.LE.0.) GO TO 20
         GX(M)=B2*PX(2) + B3*PX(3) + B4*PX(4)
         GY(M)=B2*PY(2) + B3*PY(3) + B4*PY(4)
         GZ(M)=B2*PZ(2) + B3*PZ(3) + B4*PZ(4)
   30 CONTINUE
C**** CHOIX DES PARAMETRES ATTACHES AU PLAN DE RACCORDEMENT.
      DO 40 M=5,16
         CALL RANDOM (YFL) 
         C(M)=4.*YFL - 2. 
   40 CONTINUE
C**** TEST DES DEUX TETRAEDRES M=1 ET M=5.
      DO 200 M=1,5,4
   50    CALL RANDOM (YFL) 
         PX(M)=4.*YFL - 2. 
         CALL RANDOM (YFL) 
         PY(M)=4.*YFL - 2. 
         CALL RANDOM (YFL) 
         PZ(M)=4.*YFL - 2.
         SC= RX*(PX(M) - (PX(2)+PX(3)+PX(4))/3. )
     1      +RY*(PY(M) - (PY(2)+PY(3)+PY(4))/3. )
     2      +RZ*(PZ(M) - (PZ(2)+PZ(3)+PZ(4))/3. )
         IF (M.EQ.1  .AND. SC.LT.0.1) GO TO 50
         IF (M.EQ.5  .AND. SC.GT.-0.1) GO TO 50
         WRITE(UNIT=*, FMT=*) ' *** MOSAIQUE : FINT, DFINT'
         CALL ELEM0 (PX(M),PY(M),PZ(M), PX(2),PY(2),PZ(2), 
     1               PX(3),PY(3),PZ(3), PX(4),PY(4),PZ(4)) 
         CALL ELEM1
         CALL ELEM2
         DO 70 N=1,4
            CALL RANDOM (YFL)
            C(N)=4.*YFL - 2.
   70    CONTINUE
         CALL FE2P16 (C, S)
         DO 80 N=1,5
            E(N   )= FINT (GX(N),GY(N),GZ(N), S)
            E(N+ 5)=DFINT (1., 0., 0., GX(N),GY(N),GZ(N), S) 
            E(N+10)=DFINT (0., 1., 0., GX(N),GY(N),GZ(N), S)
            E(N+15)=DFINT (0., 0., 1., GX(N),GY(N),GZ(N), S) 
   80    CONTINUE
         WRITE(UNIT=*, FMT=350)  (E(N), N=1,20)
  200 CONTINUE
      RETURN
  350 FORMAT (20F20.12)
      END
C---------------------------------------------------
C
      SUBROUTINE PERM (AX, AY, AZ)
C
      DIMENSION AX(4),AY(4),AZ(4), S(28), F0(4),FX(4),FY(4),FZ(4), 
     1          GX(15),GY(15),GZ(15), E(30), C(28) 
C**** CHOIX DES PARAMETRES.
      DO 10 M=1,4
         CALL RANDOM (YFL)
         F0(M)=2.*YFL - 1.
         CALL RANDOM (YFL)
         FX(M)=2.*YFL - 1.
         CALL RANDOM (YFL)
         FY(M)=2.*YFL - 1.
         CALL RANDOM (YFL)
         FZ(M)=2.*YFL - 1.
   10 CONTINUE
C**** CHOIX DE POINTS DANS LES FACES.
C**** FACE 1.
      DO 30 M=1,3
   20    CALL RANDOM (B2)
         CALL RANDOM (B3)
         B4=1.-B2-B3
         IF (B4.LE.0.) GO TO 20
         GX(M)=B2*AX(2) + B3*AX(3) + B4*AX(4) 
         GY(M)=B2*AY(2) + B3*AY(3) + B4*AY(4) 
         GZ(M)=B2*AZ(2) + B3*AZ(3) + B4*AZ(4)
   30 CONTINUE
C**** FACE 2.
      DO 50 M=4,6
   40    CALL RANDOM (B3)
         CALL RANDOM (B4)
         B1=1.-B3-B4
         IF (B1.LE.0.) GO TO 40
         GX(M)=B3*AX(3) + B4*AX(4) + B1*AX(1) 
         GY(M)=B3*AY(3) + B4*AY(4) + B1*AY(1) 
         GZ(M)=B3*AZ(3) + B4*AZ(4) + B1*AZ(1)
   50 CONTINUE
C**** FACE 3.
      DO 70 M=7,9
   60    CALL RANDOM (B4)
         CALL RANDOM (B1)
         B2=1.-B4-B1
         IF (B2.LE.0.) GO TO 60
         GX(M)=B4*AX(4) + B1*AX(1) + B2*AX(2) 
         GY(M)=B4*AY(4) + B1*AY(1) + B2*AY(2) 
         GZ(M)=B4*AZ(4) + B1*AZ(1) + B2*AZ(2)
   70 CONTINUE
C**** FACE 4.
      DO 90 M=10,12
   80    CALL RANDOM (B1)
         CALL RANDOM (B2)
         B3=1.-B1-B2
         IF (B3.LE.0.) GO TO 80
         GX(M)=B1*AX(1) + B2*AX(2) + B3*AX(3) 
         GY(M)=B1*AY(1) + B2*AY(2) + B3*AY(3) 
         GZ(M)=B1*AZ(1) + B2*AZ(2) + B3*AZ(3)
   90 CONTINUE
C**** CHOIX DE POINTS A L'INTERIEUR.
      DO 110 M=13,15
  100    CALL RANDOM (B1) 
         CALL RANDOM (B2) 
         CALL RANDOM (B3) 
         B4=1.-B1-B2-B3
         IF (B4.LE.0.) GO TO 100
         GX(M)=B1*AX(1) + B2*AX(2) + B3*AX(3) + B4*AX(4) 
         GY(M)=B1*AY(1) + B2*AY(2) + B3*AY(3) + B4*AY(4) 
         GZ(M)=B1*AZ(1) + B2*AZ(2) + B3*AZ(3) + B4*AZ(4)
  110 CONTINUE
C**** TEST DES 24 PERMUTATIONS DES SOMMETS DU TETRAEDRE.
      DO 200 I=1,4
         DO 190 J=1,4
            IF (J.EQ.I) GO TO 190
            DO 180 K=1,4
               IF (K.EQ.I .OR. K.EQ.J) GO TO 180
               DO 170 L=1,4
                  IF (L.EQ.I .OR. L.EQ.J .OR. L.EQ.K) GO TO 170
                  WRITE(UNIT=*, FMT=*) ' *** PERM : ',I,J,K,L,
     1                                 ' FINT, DFINT'
                  CALL ELEM0   (AX(I),AY(I),AZ(I), AX(J),AY(J),AZ(J),
     1                          AX(K),AY(K),AZ(K), AX(L),AY(L),AZ(L)) 
                  CALL ELEM1
                  CALL ELEM2 
                  C(1)=F0(I) 
                  C(2)=FX(I)
                  C(3)=FY(I) 
                  C(4)=FZ(I) 
                  C(5)=F0(J) 
                  C(6)=FX(J) 
                  C(7)=FY(J) 
                  C(8)=FZ(J) 
                  C(9)=F0(K) 
                  C(10)=FX(K) 
                  C(11)=FY(K) 
                  C(12)=FZ(K) 
                  C(13)=F0(L) 
                  C(14)=FX(L) 
                  C(15)=FY(L) 
                  C(16)=FZ(L)
                  CALL FE2P16 (C, S)
                  DO 150 M=1,15
                     IND=2*M-1
                     E(IND  )= FINT(GX(M),GY(M),GZ(M), S) 
                     E(IND+1)=DFINT(1., -3., 2., GX(M),GY(M),GZ(M), S) 
  150             CONTINUE
                  WRITE(UNIT=*, FMT=350)   (E(M), M=1,30)
  170          CONTINUE
  180       CONTINUE
  190    CONTINUE
  200 CONTINUE
      RETURN
  350 FORMAT (30F20.12) 
      END
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****
C****    VARIANTE II.                        MARCEL DELEZE
C****                                        LE CHENE 21
C****    FEVRIER 1978                        1633 MARSENS
C****    OCTOBRE 2008                        SWITZERLAND 
C****
C*************************************************************
C*************************************************************
C------------------------------------------------------------- 
C
      FUNCTION PT1 (J, WT,XT,YT,ZT)
C 
C-------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO—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=l.—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
      SUBROUTINE REFCT (WT,XT,YT,ZT, R0,R1,R2,R3)
C
C----------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO—REFERENCE. 
C**** VALEURS DES FONCTIONS RATIONNELLES R0, Rl, R2 ET R3
C**** 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 ARETE DU TETRAEDRE.
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT 
      R0=2.*X*Y*Y*Z*Z/((X+Y)*(X+Z))
      R =2.*X*X*Y*Y*Z*Z*W/((W+X)*(W+Y)*(W+Z)*(1.-Y)*(1.-Z)) 
      R1=R*(-1. - Y*Z*(1.+X)/((1.-Y)*(1.-Z)))
      R2=R*Z*(2.*X+Y)/(1.-Z)
      R3=R*Y*(2.*X+Z)/(1.-Y)
      RETURN
      END
C-------------------------------------------------------------------
C
      FUNCTION WT11 (MUX,MUY,MUZ, WT,XT,YT,ZT)
C
C-------------------------------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO-REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE WT(1,1,MU)
C**** AU POINT (XT,YT,ZT) DU TETRAEDRE DE REFERENCE.  WT=1.-XT-YT-ZT. 
C**** LA DIRECTION NORMALE DEPLACEE (MUX, MUY, MUZ) DEPEND DU TETRAEDRE 
C**** GENERIQUE.	NOUS SUPPOSONS QUE	MUX+MUY+MUZ = 1.
      REAL MUX,MUY,MUZ
      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 30
C**** DANS LES FACES 2, 3 ET 4, ON A 
      WT11=0.0
      RETURN
   30 CONTINUE
      IF (W.GT.EPS) GO TO 40
C**** A L'INTERIEUR DE LA FACE 1, ON A 
      WT11=2.*X*Y*Y*Z*Z/((X+Y)*(X+Z) ) 
      RETURN
   40 CONTINUE
C**** A L'INTERIEUR DU TETRAEDRE, ON A
      CALL REFCT(W,X,Y,Z,R0,R1,R2,R3) 
      WT11=R0 + MUX*R1 + MUY*R2 + MUZ*R3 
      RETURN
      END
c-----------------------------------------------------------------------------
C
      SUBROUTINE UT28 (NFCT, XT,YT,ZT, UT) 
C
C---------------------------–-----------------–----–-------------------------- 
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO—REFERENCE.
C**** VALEUR DE UT(1), ..., UT(NFCT) AU POINT (XT,YT,ZT) DU TETRAEDRE DE 
C**** REFERENCE.
C**** NFCT PEUT PRENDRE LES VALEURS 16 OU 28.
C**** LES DIRECTIONS NORMALES DEPLACEES MUX(I), MUY(I), MUZ(I), I=1(1)4, 
C**** DEPENDENT DU TETRAEDRE GENERIQUE.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3) 
      DIMENSION UT(28)
      NFCS=NFCT
      XL=XT
      YL=YT
      ZL=ZT
      WL=1.-XL-YL-ZL
      DO 90 I=1,4
C****	 BLOC 1.	FONCTIONS NUMEROS 1 A 16. 
         IND=4*I-3
         UT(IND  )=PT1(0,WL,XL,YL,ZL) 
         UT(IND+1)=PT1(1,WL,XL,YL,ZL) 
         UT(IND+2)=PT1(1,WL,YL,ZL,XL) 
         UT(IND+3)=PT1(1,WL,ZL,XL,YL) 
         IF (NFCS.LE.16) GO TO 80
C****	 BLOC 2.	FONCTIONS NUMEROS 17 A 28. 
         IND=3*I+14
         UT(IND  )=WT11(MUX(I),MUY(I),MUZ(I),WL,XL,YL,ZL)
         UT(IND+1)=WT11(MUY(I),MUZ(I),MUX(I),WL,YL,ZL,XL) 
         UT(IND+2)=WT11(MUZ(I),MUX(I),MUY(I),WL,ZL,XL,YL)
   80    CH=WL 
         WL=XL
         XL=YL 
         YL=ZL 
         ZL=CH
   90 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 PSEUDO—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
      SUBROUTINE DREFCT (TX,TY,TZ, WT,XT,YT,ZT, DR0,DR1,DR2,DR3)
C
C----------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO—REFERENCE.
C**** DERIVEES DES FONCTIONS RATIONNELLES R0, R1, R2 ET R3
C**** DANS LA DIRECTION (TX,TY,TZ) 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 ARETE DU TETRAEDRE.
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT 
      T1=TX
      T2=TY 
      T3=TZ 
      T0=-T1-T2-T3
C**** DERIVEE DE R0.
      P =2.*X*Y*Y*Z*Z
      DP =2.*Y*Z*(T1*Y*Z+2.*X*(T2*Z+T3*Y))
      Q =(X+Y)*(X+Z)
      DR0=(DP -  P*( (T1+T2)/(X+Y)+(T1+T3)/(X+Z) ))/Q
C**** DERIVEE DE R.
      P  =2.*X*X*Y*Y*Z*Z*W
      DP =2.*X*Y*Z*( T0*X*Y*Z+2.*W*(T1*Y*Z+T2*X*Z+T3*X*Y) )
      Q  =(W+X)*(W+Y)*(W+Z)*(1.-Y)*(1.-Z)
      R  =P/Q
      DR =(DP  -   P*( (T0+T1)/(W+X)+(T0+T2)/(W+Y)+(T0+T3)/(W+Z)
     1	               -T2/(1.-Y)-T3/(1.-Z)  ))/Q
C**** DERIVEE DE R1.
      P	 = R*Y*Z*(1.+X)
      DP =DR*Y*Z*(1.+X) + R*( T1*Y*Z+(1.+X)*(T2*Z+T3*Y) )
      Q	 =(1.-Y)*(1.-Z)
      DR1=-DR  -  (DP - P*( -T2/(1.-Y)-T3/(1.-Z)  ))/Q
C**** DERIVEE DE R2.
      P  =R*Z*(2.*X+Y)
      DP =DR*Z*(2.*X+Y) + R*( (T1*2.+T2)*Z+T3*(2.*X+Y) )
      Q  =1.-Z
      DR2=(DP - P*(-T3/Q))/Q
C**** DERIVEE DE R3.
      P  =R*Y*(2.*X+Z)
      DP =DR*Y*(2.*X+Z) + R*( Y*(2.*T1+T3)+T2*(2.*X+Z) )
      Q	 =1.-Y
      DR3=(DP - P*(-T2/Q))/Q
      RETURN
      END
C------------------------------------------------------------------
C
      FUNCTION DWT11 (TX,TY,TZ, MUX,MUY,MUZ, WT,XT,YT,ZT)
C
C------------------------------------------------------------------	 
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO-REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE WT(1,1,MU) DANS LA DIRECTION 
C**** (TX,TY,TZ) AU POINT (XT,YT,ZT) DU TETRAEDRE DE REFERENCE.
C**** WT=1.-XT-YT-ZT.
C**** LA DIRECTION NORMALE DEPLACEE (MUX, MUY, MUZ) DEPEND DU TETRAEDRE 
C**** GENERIQUE. NOUS SUPPOSONS QUE	MUX+MUY+MUZ = 1.
      REAL MUX,MUY,MUZ
      DATA  EPS /0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (Y.GT.EPS .AND. W+Y.GT.EPS .AND.
     1    Z.GT.EPS .AND. W+Z.GT.EPS) GO TO 30
C**** DANS LES FACES 3 ET 4, ON A
      DWT11=0.0
      RETURN
   30 CONTINUE
      T1=TX 
      T2=TY 
      T3=TZ
      IF (X.GT.EPS .AND. W+X.GT.EPS) GO TO 40
C**** DANS LA FACE 2, ON A
      DWT11=2.*T1*Y*Z
      RETURN
   40 CONTINUE
C**** A L'INTERIEUR DU TETRAEDRE ET A L'INTERIEUR DE LA FACE 1, ON A 
      CALL DREFCT(T1,T2,T3,W,X,Y,Z,DR0,DR1,DR2,DR3)
      DWT11=DR0 + MUX*DR1 + MUY*DR2 + MUZ*DR3
      RETURN
      END
C----------------------------------------------------------------------
C
      SUBROUTINE DUT28 (NFCT, TX,TY,TZ, XT,YT,ZT, DUT)
C
C-----------------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO-REFERENCE.
C**** DERIVEE DE UT(1), ..., UT(NFCT) DANS LA DIRECTION (TX,TY,TZ) 
C**** AU POINT (XT,YT,ZT) DU TETRAEDRE DE REFERENCE.
C**** NFCT PEUT PRENDRE LES VALEURS 16 OU 28.
C**** LES DIRECTIONS NORMALES DEPLACEES MUX(I), MUY(I), MUZ(I),  I=1(1)4,
C**** DEPENDENT DU TETRAEDRE GENERIQUE.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	      NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	      D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3) 
      DIMENSION DUT(28)
      NFCS=NFCT 
      XL=XT
      YL=YT
      ZL=ZT
      WL=1.-XL-YL-ZL
      T1=TX
      T2=TY
      T3=TZ
      T0=-T1-T2-T3
      DO 90 I=1,4
C****    BLOC 1.	FONCTIONS NUMEROS 1 A 16.
         IND=4*I-3
         DUT(IND  )=DPT1(0,T1,T2,T3,WL,XL,YL,ZL) 
         DUT(IND+1)=DPT1(1,T1,T2,T3,WL,XL,YL,ZL) 
         DUT(IND+2)=DPT1(1,T2,T3,T1,WL,YL,ZL,XL) 
         DUT(IND+3)=DPT1(1,T3,T1,T2,WL,ZL,XL,YL)
         IF (NFCS.LE.16) GO TO 80
C****    BLOC 2.	FONCTIONS NUMEROS 17 A 28.
         IND=3*I+14
         DUT(IND  )=DWT11(T1,T2,T3,MUX(I),MUY(I),MUZ(I),WL,XL,YL,ZL)
         DUT(IND+1)=DWT11(T2,T3,T1,MUY(I),MUZ(I),MUX(I),WL,YL,ZL,XL)
         DUT(IND+2)=DWT11(T3,T1,T2,MUZ(I),MUX(I),MUY(I),WL,ZL,XL,YL)
   80    CH=WL 
         WL=XL 
         XL=YL 
         YL=ZL 
         ZL=CH 
         CH=T0 
         T0=T1 
         T1=T2 
         T2=T3 
         T3=CH
   90 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 DU
C****                    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****	    NX(I,J), NY(I,J), NZ(I,J), J=1,2,3,  I=1(1)4, VECTEURS 
C****               PARALLELES A LA I-EME FACE ET NORMAUX AUX
C****	            ARETES DU TETRAEDRE.
C****       MUX(I), MUY(I), MUZ(I), I=1(1)4, DIRECTIONS NORMALES 
C****               DEPLACEES. CHAQUE VECTEUR A ETE MULTIPLIE PAR 
C****	            UN SCALAIRE TEL QUE LE PRODUIT SATISFASSE
C****               MUX(I) + MUY(I) + MUZ(I) = 1,	I=1(1)4.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	     D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION GEN(3),REF(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)
      CALL DECOMP(3,3,LL,GEN,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 60 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1 
         L=MOD(K,4)+1
C****	 VECTEUR NORMAL A LA I-EME FACE, PLACE DANS GEN. 
         GEN(1)=(Y(K)-Y(J))*(Z(L)-Z(J))-(Z(K)-Z(J))*(Y(L)-Y(J)) 
         GEN(2)=(Z(K)-Z(J))*(X(L)-X(J))-(X(K)-X(J))*(Z(L)-Z(J)) 
         GEN(3)=(X(K)-X(J))*(Y(L)-Y(J))-(Y(K)-Y(J))*(X(L)-X(J))
C****    VECTEURS PARALLELES A LA I-EME FACE ET NORMAUX AUX ARETES. 
         NX(I,1)=GEN(2)*(Z(L)-Z(K))-GEN(3)*(Y(L)-Y(K)) 
         NY(I,1)=GEN(3)*(X(L)-X(K))-GEN(1)*(Z(L)-Z(K)) 
         NZ(I,1)=GEN(1)*(Y(L)-Y(K))-GEN(2)*(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)=GEN(2)*(Z(J)-Z(L))-GEN(3)*(Y(J)-Y(L)) 
         NY(I,2)=GEN(3)*(X(J)-X(L))-GEN(1)*(Z(J)-Z(L))
         NZ(I,2)=GEN(1)*(Y(J)-Y(L))-GEN(2)*(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)=GEN(2)*(Z(K)-Z(J))-GEN(3)*(Y(K)-Y(J)) 
         NY(I,3)=GEN(3)*(X(K)-X(J))-GEN(1)*(Z(K)-Z(J)) 
         NZ(I,3)=GEN(1)*(Y(K)-Y(J))-GEN(2)*(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
C****	 DIRECTION NORMALE DEPLACEE RELATIVE A LA I-EME FACE.
         CALL SOLVE (3,3,LL,GEN,REF,IPERM)
         GO TO (10,20,30,40), I
   10    MUX(1)=REF(1) 
         MUY(1)=REF(2) 
         MUZ(1)=REF(3) 
         GO TO 50
   20    MUX(2)=REF(2) 
         MUY(2)=REF(3) 
         MUZ(2)=-REF(1)-REF(2)-REF(3)
         GO TO 50
   30	 MUX(3)=REF(3) 
         MUY(3)=-REF(1)-REF(2)-REF(3)
         MUZ(3)=REF(1) 
         GO TO 50
   40	 MUX(4)=-REF(1)-REF(2)-REF(3)
         MUY(4)=REF(1) 
         MUZ(4)=REF(2)
   50	 S=MUX(I)+MUY(I)+MUZ(I)
         MUX(I)=MUX(I)/S 
         MUY(I)=MUY(I)/S 
         MUZ(I)=MUZ(I)/S
   60 CONTINUE
      RETURN
      END
C---------------------------------------------------------------
C
      SUBROUTINE U28 (NFCT, PX,PY,PZ, U)
C 
C---------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** VALEUR DE U(1), ..., U(NFCT) AU POINT (PX,PY,PZ) DU TETRAEDRE. 
C**** NFCT PEUT PRENDRE LES VALEURS 16 OU 28.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	     NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	     D(4,3,3),D5(12),Q1(12,16),T(28,16), LL(3,3),IPERM(3) 
      DIMENSION U(28),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 U28 :',
     2                ' THE POINT IS OUTSIDE THE ELEMENT.'
      NERR=NERR+1
   20 CONTINUE
C**** VALEURS DES FONCTIONS DE BASE.
      CALL UT28(NFCT,XT,YT,ZT,U)
      RETURN
      END
C------------------------------------------------------------------- 
C
      SUBROUTINE DU28 (NFCT, VX,VY,VZ, PX,PY,PZ, DU)
C 
C-------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** DERIVEE DE U(1), .., U(NFCT) DANS LA DIRECTION (VX, VY, VZ) 
C**** AU POINT (PX,PY,PZ) DU TETRAEDRE GENERIQUE.
C**** NFCT PEUT PRENDRE LES VALEURS 16 OU 28.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	     NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	     D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION DU(28),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 DU28 :',
     2                ' THE POINT IS OUTSIDE THE ELEMENT.'
      NERR=NERR+1 
   20 CONTINUE
C**** DERIVEES DES FONCTIONS DE BASE.
      CALL DUT28(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, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	     NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	     D(4,3,3),D5(12),Q1(12,16), T(28,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 OU 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))
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------
C
      SUBROUTINE ELEM2
C 
C---------------------------------------------------------------
C**** MATRICE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** INITIALISATION DU BLOC Q1 DE Q. 
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ 
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION DU(28)
      DO 40 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1 
         L=MOD(K,4)+1 
         BX=0.5*(X(K)+X(L))
         BY=0.5*(Y(K)+Y(L))
         BZ=0.5*(Z(K)+Z(L))
         CALL DU28(16,NX(I,1),NY(I,1),NZ(I,1),BX,BY,BZ,DU)
         IND=3*I-2
         DO 10 M=1,16
            Q1(IND,M)=DU(M) 
   10    CONTINUE
         BX=0.5*(X(J)+X(L)) 
         BY=0.5*(Y(J)+Y(L)) 
         BZ=0.5*(Z(J)+Z(L))
         CALL DU28(16,NX(I,2),NY(I,2),NZ(I,2),BX,BY,BZ,DU)
         IND=IND+1
         DO 20 M=1,16
            Q1(IND,M)=DU(M) 
   20    CONTINUE
         BX=0.5*(X(K)+X(J)) 
         BY=0.5*(Y(K)+Y(J)) 
         BZ=0.5*(Z(K)+Z(J))
         CALL DU28(16,NX(I,3),NY(I,3),NZ(I,3),BX,BY,BZ,DU)
         IND=IND+1
         DO 30 M=1,16
            Q1(IND,M)=DU(M) 
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
C-------------------------------------------------------------------	
C
      SUBROUTINE FE2P28 (C, S)
C 
C-------------------------------------------------------------------
C**** ELEMENT GENERIQUE ELARGI.
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(28) VALEURS DES PARAMETRES DE L'ELEMENT.
C**** SORTIES : S(1), ..., S(28) 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, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION C(28),S(28)
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
      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(28)
      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**** SORTIES  :  S(1), ..., S(28) 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(28) SONT DES PLACES DE TRAVAIL. 
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION C(28),S(28)
C**** CALCUL DES PARAMETRES C(17), ..., C(28).
      DO 10 I=1,4
         J=MOD(I,4)+1
         K=MOD(J,4)+1
         L=MOD(K,4)+1
         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))
   10 CONTINUE
C**** CALCUL DES COORDONNEES S(1), ..., S(28).
      CALL FE2P28(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(28), 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(28),S(28)
      CALL U28(28,PX,PY,PZ,U)
      V=0.0
      DO 10 I=1,28
         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(28), 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(28),S(28)
      CALL DU28(28,VX,VY,VZ,PX,PY,PZ,DU)
      DV=0.0
      DO 10 I=1,28
         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, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION C(28), S(28)
      DATA C/28*0.0/
      DO 30 I=1,16
C****    CALCUL DES COEFFICIENTS T(1,I), ... , T(28,I) DU I-EME ELEMENT
C****    DE LA BASE D'HERMITE PAR RAPPORT A LA BASE D'INTERPOLATION DE
C****    L'ELEMENT GENERIQUE ELARGI.
         C(I)=1.0
         CALL FE2P16(C,S)
         DO 20 J=1,28
            T(J,I)=S(J)
   20    CONTINUE
         C(I)=0.0 
   30 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**** SORTIES : V(1), ..., V(16) VALEURS DES 16 FONCTIONS DE LA BASE
C****	               D'HERMITE AU POINT (PX, PY, PZ).
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION V(16), U(28)
      CALL U28(28,PX,PY,PZ,U)
      DO 20 I=1,16
         SUM=0.0
         DO 10 J=1,28
            SUM=SUM + T(J,I)*U(J)
   10    CONTINUE 
         V(I)=SUM 
   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**** 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**** SORTIES  :  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, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION DV(16), DU(28)
      CALL DU28(28,VX,VY,VZ,PX,PY,PZ,DU)
      DO 20 I=1,16
         SUM=0.0
         DO 10 J=1,28
            SUM=SUM + T(J,I)*DU(J)
   10    CONTINUE 
         DV(I)=SUM
   20 CONTINUE 
      RETURN 
      END
c-------------------------------------------------------------    
C
      SUBROUTINE DECOMP (NX,NN,A,SCALES,IPS,D) 
C
C     THIS ROUTINE DOES AN LU DECOMPOSITION ON A SQUARE MATRIX.
C     IT USES ROW SCALING AND ROW PIVOTING TO MAKE THE PROCESS
C     AS ACCURATE AS POSSIBLE.
C     THE PARAMETERS OF THE CALLING SEQUENCE ARE
C     N       THE ORDER OF THE MATRIX TO BE DECOMPOSED
C     NX      THE NUMBER OF LINES OF THE MATRIX A
C	      IN THE CALLING PROGRAM
C     A       THE MATRIX TO BE DECOMPOSED
C	      THE DECOMPOSED MATRIX. THE MAIN DIAGONAL AND EVERYTHING
C	      ABOVE IT MAKE UP THE UPPER TRIANGULAR MATRIX. THE LOWER
C	      TRIANGULAR MATRIX IS MADE UP OF EVERYTHING BELOW THE MAIN
C	      DIAGONAL. THE MAIN DIAGONAL FOR THE LOWER TRIANGULAR
C	      MATRIX IS ASSUMED TO BE ALL ONES.
C     D       THE DETERMINANT OF THE COEFFICIENT MATRIX
C     KERR    ERROR PRINT CONTROL AND RETURN ERROR FLAG
C             IF KERR SET NEGATIVE BY USER NO MESSAGE PRINTED
C             IF SET TO ZERO, ERROR MESSAGE WILL BE PRINTED.
C             ON RETURN, FOLLOWING CODE -
C             =1-MATRIX WITH ZERO ROW
C             =2-SINGULAR MATRIX
C             =3- SOLUTION OBTAINED
C     IPS    VECTOR CONTAINING PERMUTATED ROW NUMBERS OF A 
C
      DIMENSION A(NX,NX), SCALES(NX), IPS(NX)
      N=NN
      KERR=0
C
C     FIND THE SOLUTION IN THE CASE N=1
C
      IF(N.GT.1) GO TO 20 
      D=A(1,1)
      IPS(1)=1
      IF(D.EQ.0.) GO TO 10
      KERR=3 
      RETURN
   10 IF(KERR.EQ.0) CALL SING(1)
      KERR=1 
      RETURN
C
C	INITIALIZE
C
   20 ODD=1.
C
C	FIND THE SCALING FACTOR FOR EACH ROW AND STORE IN SCALES
C
      DO 40 I=1,N 
      IPS(I)=I 
      RMAX=0.
      DO 30 J=1,N
         RMAX=AMAX1(RMAX,ABS(A(I,J))) 
   30 CONTINUE
      IF(RMAX.EQ.0.) GO TO 110 
      SCALES(I)=1./RMAX
   40 CONTINUE
C
C     MAIN LOOP. GAUSSIAN ELIMINATION IS CARRIED OUT. LOOP IS GONE
C     THROUGH ONCE FOR EACH ELEMENT ON THE MAIN DIAGONAL . K IS THE
C     DIAGONAL ELEMENT NUMBER. LAST DIAGONAL ELEMENT IS AN EXCEPTION. 
C
      NM1=N-1
      DO 90 K=1,NM1 
         BIG=0.
         DO 50 I=K,N 
            IP=IPS(I) 
            SIZE=ABS(A(IP,K))*SCALES(IP)
            IF(BIG.GT.SIZE) GO TO 50
            BIG=SIZE
            IPIV=I
   50    CONTINUE
         IF(BIG.EQ.0.) GO TO 120
         IF(IPIV.EQ.K) GO TO 60
         J=IPS(K)
         IPS(K)=IPS(IPIV)
         IPS(IPIV)=J
         ODD=-ODD
   60    CONTINUE 
         KP=IPS(K) 
         PIVOT=A(KP,K)
         KP1=K+1
         DO 80 I=KP1,N
            IP=IPS(I)
            EM=A(IP,K)/PIVOT
            A(IP,K)=EM
            IF(EM.EQ.0.) GO TO 80
            DO 70 J=KP1,N
               A(IP,J)=A(IP,J)-EM*A(KP,J) 
   70       CONTINUE
   80    CONTINUE
   90 CONTINUE
      KP=IPS(N)
      IF(A(KP,N).EQ.0.) GO TO 120
C
C	CALCULATE THE DETERMINANT
C
      D=1.
      DO 100 I=1,N
         IP=IPS(I)
         D=D*A(IP,I)
  100 CONTINUE 
      D=D*ODD 
      KERR=3 
      RETURN
C
C	ERROR EXIST
C
  110 IF(KERR.EQ.0) CALL SING(1) 
      D=0.
      KERR=1
      RETURN
  120 IF(KERR.EQ.0) CALL SING(2) 
      D=0.
      KERR=2 
      RETURN 
      END
C
      SUBROUTINE SING(IWHY)
C
      GO TO (10,30),IWHY
   10 WRITE(UNIT=*, FMT=*) ' *** MATRIX WITH ZERO ROW IN DECOMPOSE'
      RETURN
   30 WRITE(UNIT=*, FMT=*) ' *** SINGULAR MATRIX IN DECOMPOSE.'
      WRITE(UNIT=*, FMT=*) '     ZERO DIVIDE IN SOLVE.'
      RETURN
      END
C
      SUBROUTINE SOLVE (NX,NN,A,B,X,IPS) 
C
C     THIS ROUTINE DOES FORWARD ELIMINATION AND BACK SUBSTITUTION
C     TO SOLVE THE EQUATION LU*X=B
C     THE PARAMETERS OF THE CALLING SEQUENCE ARE
C     N       THE ORDER OF LU
C     NX      THE NUMBER OF LINES OF THE MATRIX LU
C	      IN THE CALLING PROGRAM
C     A       THE DECOMPOSED MATRIX, L AND U ARE PACKED
C	      INTO ONE MATRIX BY DECOMP
C     B       THE RIGHT HAND SIDE OF THE EQUATION
C     IPS     VECTOR CONTAINING PERMUTATED ROW NUMBERS OF A
C	      THIS VECTOR IS BY DECOMP
C
      DIMENSION A(NX,NX), B(NX), X(NX), IPS(NX)
      N=NN
C
C     FIND THE SOLUTION IN THE CASE N=1
C
      IF(N.GT.1) GO TO 10
      X(1)=B(1)/A(1,1)
      RETURN
C
C     FORWARD ELIMINATION USING L AND STORING THE RESULT IN X
C
   10 NP1=N+1
      IP=IPS(1)
      X(1)=B(IP) 
      DO 30 I=2,N 
         IP=IPS(I)
         IM1=I-1
         SUM=0.
         DO 20 J=1,IM1
            SUM=SUM+A(IP,J)*X(J) 
   20    CONTINUE
         X(I)=B(IP)-SUM
   30 CONTINUE
C
C	BACK SUBSTITUTION USING U AND X
C
      IP=IPS(N) 
      X(N)=X(N)/A(IP,N)
      DO 50 IBACK=2,N
         I=NP1-IBACK
         IP=IPS(I)
         IP1=I+1 
         SUM=0.
         DO 40 J=IP1,N
               SUM=SUM+A(IP,J)*X(J)
   40    CONTINUE 
         X(I)=(X(I)-SUM)/A(IP,I)
   50 CONTINUE
      RETURN
      END
C------------------------------------------
      SUBROUTINE RANDOM(ALEAT)
C
C**** GENERATEUR DE NOMBRES ALEATOIRES UNIFORMEMENT DISTRIBUES
C**** ENTRE 0.0 ET 1.0
      DATA IA,IX,IC/6133,4409,2741/
      DATA M,RM/81799,81799./
      IX=MOD(IA*IX+IC,M)
      X=IX
      ALEAT=X/RM
      RETURN
      END
