i15tot1.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    i15tot1        ../engine/source/interfaces/int15/i15tot1.F
!||--- called by ------------------------------------------------------
!||    i15cmp         ../engine/source/interfaces/int15/i15cmp.F
!||--- uses       -----------------------------------------------------
!||    groupdef_mod   ../common_source/modules/groupdef_mod.F
!||====================================================================
      SUBROUTINE i15tot1(NOINT ,NDEB  ,NSC  ,X     ,KSURF ,
     2                  IGRSURF,BUFSF ,KSC  ,KSI   ,NOLD  ,
     3                  XP1   ,XP2    ,XP3  ,XP4   ,GX    ,
     4                  XTK   ,YTK    ,ZTK  ,NTX   ,NTY   ,
     5                  NTZ   ,PENET ,DEPTH ,XI    ,YI    ,
     6                  ZI    ,NXI   ,NYI   ,NZI   ,ANSMX ,
     7                  HOLD  ,IACTIV,ITAB  )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE groupdef_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
#include      "units_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER KSURF ,KSI(4,*) ,
     .        NOINT ,NDEB, NSC, IACTIV(4,*), KSC(*), ITAB(*)
C     REAL
      my_real 
     .  X(3,*)    ,BUFSF(*)  ,
     .  NOLD(3,*) ,XTK(4,*)  ,YTK(4,*) ,ZTK(4,*) ,
     .  ntx(4,*)  ,nty(4,*)  ,ntz(4,*) ,
     .  penet(4,*),depth(4,*),
     .  xi(4,*)   ,yi(4,*)  ,zi(4,*)  ,
     .  nxi(4,*)  ,nyi(4,*) ,nzi(4,*) ,
     .  xp1(3,mvsiz), xp2(3,mvsiz), xp3(3,mvsiz), xp4(3,mvsiz),
     .  gx(3,mvsiz), ansmx, hold(3,*)
      TYPE (SURF_)   , DIMENSION(NSURF)   :: IGRSURF
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER ADRBUF, I, IL, NLS, IDG,
     .        IN1, IN2, IN3, IN4,
     .        inside1,inside2
      my_real
     .   a, b, c, an, bn, cn, rot(9), dgr, expn,
     .   xg, yg, zg,
     .   ntn,
     .   n1, n2, n3, n, n11, n12, n13, n21, n22, n23, nr1, nr2,
     .   xkn1, ykn1, zkn1, sgnxkn, sgnykn, sgnzkn, xkn, ykn, zkn, eh,
     .   lambda1, lambda2, alp, bet,
     .   xh , yh , zh , mu, xh1, yh1, zh1, mu1, xh2, yh2, zh2, mu2,
     .   dx1, dy1, dz1, dx2, dy2, dz2, dx3, dy3, dz3,
     .   side1, side2, oldnx, oldny, oldnz, out, ps, nr,
     .   lambda, em, ep
      my_real
     .   x1(mvsiz),  y1(mvsiz),  z1(mvsiz),
     .   x2(mvsiz),  y2(mvsiz),  z2(mvsiz),
     .   x3(mvsiz),  y3(mvsiz),  z3(mvsiz),
     .   x12(mvsiz), y12(mvsiz), z12(mvsiz),
     .   x13(mvsiz), y13(mvsiz), z13(mvsiz),
     .   x23(mvsiz), y23(mvsiz), z23(mvsiz),
     .   xm(mvsiz) , ym(mvsiz) , zm(mvsiz) , 
     .   xtk2(mvsiz), ytk2(mvsiz) , ztk2(mvsiz) , 
     .   d(mvsiz)
C-----------------------------------------------
      adrbuf=igrsurf(ksurf)%IAD_BUFR
      a =bufsf(adrbuf+1)
      b =bufsf(adrbuf+2)
      c =bufsf(adrbuf+3)
      dgr =bufsf(adrbuf+36)
      idg =nint(dgr)
      dgr =idg
      dgr =max(two,dgr)
      expn=dgr/(dgr-1.)
      an=a**dgr
      bn=b**dgr
      cn=c**dgr
      an=one/max(em20,an)
      bn=one/max(em20,bn)
      cn=one/max(em20,cn)
      DO i=1,9
       rot(i)=bufsf(adrbuf+7+i-1)
      END DO
C-------------------------------
C     Passage au repere local :
C-------------------------------
      DO 75 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
      in1=ksi(1,il)
      in2=ksi(2,il)
      in3=ksi(3,il)
      in4=ksi(4,il)
C     Passage au repere local.
      xg=x(1,in1)-bufsf(adrbuf+4)
      yg=x(2,in1)-bufsf(adrbuf+5)
      zg=x(3,in1)-bufsf(adrbuf+6)
      xp1(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
      xp1(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
      xp1(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
      xg=x(1,in2)-bufsf(adrbuf+4)
      yg=x(2,in2)-bufsf(adrbuf+5)
      zg=x(3,in2)-bufsf(adrbuf+6)
      xp2(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
      xp2(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
      xp2(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
      xg=x(1,in3)-bufsf(adrbuf+4)
      yg=x(2,in3)-bufsf(adrbuf+5)
      zg=x(3,in3)-bufsf(adrbuf+6)
      xp3(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
      xp3(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
      xp3(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
      xg=x(1,in4)-bufsf(adrbuf+4)
      yg=x(2,in4)-bufsf(adrbuf+5)
      zg=x(3,in4)-bufsf(adrbuf+6)
      xp4(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
      xp4(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
      xp4(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
      gx(1,i)=fourth*(xp1(1,i)+xp2(1,i)+xp3(1,i)+xp4(1,i))
      gx(2,i)=fourth*(xp1(2,i)+xp2(2,i)+xp3(2,i)+xp4(2,i))
      gx(3,i)=fourth*(xp1(3,i)+xp2(3,i)+xp3(3,i)+xp4(3,i))     
 75   CONTINUE
C-------------------------------
C     1ER TRIANGLE.
C-------------------------------
      DO 100 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      x1(i)=gx(1,i)
      y1(i)=gx(2,i)
      z1(i)=gx(3,i)
      x2(i)=xp1(1,i)
      y2(i)=xp1(2,i)
      z2(i)=xp1(3,i)
      x3(i)=xp2(1,i)
      y3(i)=xp2(2,i)
      z3(i)=xp2(3,i)
      x12(i)=x2(i)-x1(i)
      y12(i)=y2(i)-y1(i)
      z12(i)=z2(i)-z1(i)
      x13(i)=x3(i)-x1(i)
      y13(i)=y3(i)-y1(i)
      z13(i)=z3(i)-z1(i)
      n1=y12(i)*z13(i)-z12(i)*y13(i)
      n2=z12(i)*x13(i)-x12(i)*z13(i)
      n3=x12(i)*y13(i)-y12(i)*x13(i)
      ntn=one/max(em20,sqrt(n1*n1+n2*n2+n3*n3))
      ntx(1,i)=ntn*n1
      nty(1,i)=ntn*n2
      ntz(1,i)=ntn*n3
      d(i)    =-ntx(1,i)*x1(i)-nty(1,i)*y1(i)-ntz(1,i)*z1(i)
C-----
      x23(i)=x3(i)-x2(i)
      y23(i)=y3(i)-y2(i)
      z23(i)=z3(i)-z2(i)
 100  CONTINUE
C-------------------------------
C     POINTS K,H SUR E,L / LA DISTANCE D(K,H) EST LOCALEMENT MAXIMUM
C-------------------------------
      DO 125 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
      eh  =(abs(ntx(1,i)/(dgr*an))**expn)*an
     .    +(abs(nty(1,i)/(dgr*bn))**expn)*bn
     .    +(abs(ntz(1,i)/(dgr*cn))**expn)*cn
C     X EST DU SIGNE DE LAMBDA*NTX, IDEM EN Y ET Z.
C     LAMBDA1=EXP(LOG(1./EH)/EXPN)
      lambda1=(max(em20,eh))**(-one/expn)
      xh1=abs(lambda1*ntx(1,i)/(dgr*an))**(one/(dgr-one))
      IF (ntx(1,i)<zero) xh1=-xh1
      yh1=abs(lambda1*nty(1,i)/(dgr*bn))**(one/(dgr-one))
      IF (nty(1,i)<zero) yh1=-yh1
      zh1=abs(lambda1*ntz(1,i)/(dgr*cn))**(one/(dgr-one))
      IF (ntz(1,i)<zero) zh1=-zh1
      mu1   =-ntx(1,i)*xh1-nty(1,i)*yh1-ntz(1,i)*zh1-d(i)
C     LAMBDA2=-LAMBDA1
      xh2 =-xh1
      yh2 =-yh1
      zh2 =-zh1
C     MU2   =-NTX(1,I)*XH2-NTY(1,I)*YH2-NTZ(1,I)*ZH2-D(I)
      mu2=-mu1-two*d(i)
      xtk(1,i) =xh1+mu1*ntx(1,i)
      ytk(1,i) =yh1+mu1*nty(1,i)
      ztk(1,i) =zh1+mu1*ntz(1,i)
      xtk2(i)  =xh2+mu2*ntx(1,i)
      ytk2(i)  =yh2+mu2*nty(1,i)
      ztk2(i)  =zh2+mu2*ntz(1,i)
 125  CONTINUE
C-------------------------------
C     RAMENER LE PT PK SUR LE TRIANGLE.
C-------------------------------
      DO 150 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      dx1=xtk(1,i)-x1(i)
      dy1=ytk(1,i)-y1(i)
      dz1=ztk(1,i)-z1(i)
      dx2=xtk(1,i)-x2(i)
      dy2=ytk(1,i)-y2(i)
      dz2=ztk(1,i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(1,i)
     .     +(dz1*dx2-dz2*dx1)*nty(1,i)
     .     +(dx1*dy2-dy1*dx2)*ntz(1,i)
      IF (out<0.) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk(1,i)=alp*x1(i)+bet*x2(i)
        ytk(1,i)=alp*y1(i)+bet*y2(i)
        ztk(1,i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk(1,i)-x3(i)
      dy3=ytk(1,i)-y3(i)
      dz3=ztk(1,i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(1,i)
     .     +(dz2*dx3-dz3*dx2)*nty(1,i)
     .     +(dx2*dy3-dy2*dx3)*ntz(1,i)
      IF (out<zero) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk(1,i)=alp*x2(i)+bet*x3(i)
        ytk(1,i)=alp*y2(i)+bet*y3(i)
        ztk(1,i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(1,i)
     .     +(dz3*dx1-dz1*dx3)*nty(1,i)
     .     +(dx3*dy1-dy3*dx1)*ntz(1,i)
      IF (out<zero) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk(1,i)=alp*x3(i)+bet*x1(i)
        ytk(1,i)=alp*y3(i)+bet*y1(i)
        ztk(1,i)=alp*z3(i)+bet*z1(i)
      ENDIF 
C-----
      dx1=xtk2(i)-x1(i)
      dy1=ytk2(i)-y1(i)
      dz1=ztk2(i)-z1(i)
      dx2=xtk2(i)-x2(i)
      dy2=ytk2(i)-y2(i)
      dz2=ztk2(i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(1,i)
     .      +(dz1*dx2-dz2*dx1)*nty(1,i)
     .      +(dx1*dy2-dy1*dx2)*ntz(1,i)
      IF (out<zero) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk2(i)=alp*x1(i)+bet*x2(i)
        ytk2(i)=alp*y1(i)+bet*y2(i)
        ztk2(i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk2(i)-x3(i)
      dy3=ytk2(i)-y3(i)
      dz3=ztk2(i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(1,i)
     .      +(dz2*dx3-dz3*dx2)*nty(1,i)
     .      +(dx2*dy3-dy2*dx3)*ntz(1,i)
      IF (out<zero) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk2(i)=alp*x2(i)+bet*x3(i)
        ytk2(i)=alp*y2(i)+bet*y3(i)
        ztk2(i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(1,i)
     .      +(dz3*dx1-dz1*dx3)*nty(1,i)
     .      +(dx3*dy1-dy3*dx1)*ntz(1,i)
      IF (out<zero) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk2(i)=alp*x3(i)+bet*x1(i)
        ytk2(i)=alp*y3(i)+bet*y1(i)
        ztk2(i)=alp*z3(i)+bet*z1(i)
      ENDIF 
C-------------------------------
 150  CONTINUE
C-------------------------------
C     PROJECTION DE PK SUR L ET PENETRATION.
C-------------------------------
      DO 175 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
      IF (iactiv(1,il)==-1) GOTO 175
C-----
      xh =xtk(1,i)
      yh =ytk(1,i)
      zh =ztk(1,i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one
      ykn1 =yh**(idg-1)
      sgnykn=-one
      IF (ykn1*yh>=zero) sgnykn=one
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one
      n11 =sgnxkn*xkn1*an
      n12 =sgnykn*ykn1*bn
      n13 =sgnzkn*zkn1*cn
      nr1 =n11*n11+n12*n12+n13*n13
      nr1 =sqrt(nr1)
      em =n11*xtk(1,i)+n12*ytk(1,i)+n13*ztk(1,i)
      IF (em<=one) THEN
        lambda1=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr1)
        inside1=1
      ELSE
        inside1=0
      ENDIF
C-----
      xh =xtk2(i)
      yh =ytk2(i)
      zh =ztk2(i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one
      ykn1 =yh**(idg-1)
      sgnykn=-one
      IF (ykn1*yh>=zero) sgnykn=one
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one
      n21 =sgnxkn*xkn1*an
      n22 =sgnykn*ykn1*bn
      n23 =sgnzkn*zkn1*cn
      nr2 =n21*n21+n22*n22+n23*n23
      nr2 =sqrt(nr2)
      em =n21*xtk2(i)+n22*ytk2(i)+n23*ztk2(i)
      IF (em<=one) THEN
        lambda2=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr2)
        inside2=1
      ELSE
        inside2=0
      ENDIF
C-----
      IF (inside1==0.AND.inside2==0) THEN
        iactiv(1,il)=0
      ELSE
C-----
      IF (iactiv(1,il)==0) THEN
        IF (inside1/=0.AND.inside2/=0) THEN
         IF (abs(lambda1)>=abs(lambda2)) THEN
          xm(i)=xtk(1,i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(1,i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(1,i)-lambda1*n13/max(em20,nr1)
         ELSE
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(1,i)=xtk2(i)
          ytk(1,i)=ytk2(i)
          ztk(1,i)=ztk2(i)
         ENDIF
        ELSEIF(inside1/=0) THEN
          xm(i)=xtk(1,i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(1,i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(1,i)-lambda1*n13/max(em20,nr1)
        ELSEIF(inside2/=0) THEN
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(1,i)=xtk2(i)
          ytk(1,i)=ytk2(i)
          ztk(1,i)=ztk2(i)
        ENDIF
C-----
      ELSE
        xh=hold(1,4*(il-1)+1)
        yh=hold(2,4*(il-1)+1)
        zh=hold(3,4*(il-1)+1)
        n1=nold(1,4*(il-1)+1)
        n2=nold(2,4*(il-1)+1)
        n3=nold(3,4*(il-1)+1)
        lambda1=(xh-xtk(1,i))*n1
     .         +(yh-ytk(1,i))*n2
     .         +(zh-ztk(1,i))*n3
        lambda2=(xh-xtk2(i))*n1
     .         +(yh-ytk2(i))*n2
     .         +(zh-ztk2(i))*n3
        IF (inside1/=0.AND.inside2/=0) THEN
          IF (abs(lambda1)>=abs(lambda2)) THEN
           xm(i)=xtk(1,i)+lambda1*n1
           ym(i)=ytk(1,i)+lambda1*n2
           zm(i)=ztk(1,i)+lambda1*n3
          ELSE
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(1,i)=xtk2(i)
           ytk(1,i)=ytk2(i)
           ztk(1,i)=ztk2(i)
          ENDIF
        ELSEIF(inside1/=0) THEN
           xm(i)=xtk(1,i)+lambda1*n1
           ym(i)=ytk(1,i)+lambda1*n2
           zm(i)=ztk(1,i)+lambda1*n3
        ELSEIF(inside2/=0) THEN
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(1,i)=xtk2(i)
           ytk(1,i)=ytk2(i)
           ztk(1,i)=ztk2(i)
        ENDIF
C       one more iteration.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        em=n1*xm(i)+n2*ym(i)+n3*zm(i)
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))        
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        nr =n1*n1+n2*n2+n3*n3
        nr =one/max(em20,sqrt(nr))
        n1 =n1*nr
        n2 =n2*nr
        n3 =n3*nr
        lambda1=(xm(i)-xtk(1,i))*n1
     .         +(ym(i)-ytk(1,i))*n2
     .         +(zm(i)-ztk(1,i))*n3
        xm(i)=xtk(1,i)+lambda1*n1
        ym(i)=ytk(1,i)+lambda1*n2
        zm(i)=ztk(1,i)+lambda1*n3
      ENDIF
      iactiv(1,il)=iactiv(1,il)+1
      ENDIF
 175  CONTINUE
C-----
#include "vectorize.inc"
      DO 200 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(1,il)<=0) GOTO 200
C-----
C      projection radiale de M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one        
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       em=n1*xm(i)+n2*ym(i)+n3*zm(i)
C------
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))        
C      normale a l'ellipsoide en M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one        
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       nr =n1*n1+n2*n2+n3*n3
       nr =one/max(em20,sqrt(nr))
       nxi(1,i)=n1*nr
       nyi(1,i)=n2*nr
       nzi(1,i)=n3*nr
       xi(1,i)=xm(i)
       yi(1,i)=ym(i)
       zi(1,i)=zm(i)
       depth(1,i)=xm(i)*nxi(1,i)+ym(i)*nyi(1,i)+zm(i)*nzi(1,i)
C------
       penet(1,i)=(xtk(1,i)-xm(i))*nxi(1,i)
     .           +(ytk(1,i)-ym(i))*nyi(1,i)
     .           +(ztk(1,i)-zm(i))*nzi(1,i)
       penet(1,i)=-penet(1,i)
       penet(1,i)=max(zero,penet(1,i))
       IF (depth(1,i)-penet(1,i)<em10*depth(1,i)) THEN
         penet(1,i)=zero
         iactiv(1,il)=-2
       ENDIF
       ansmx=max(ansmx,penet(1,i))
 200  CONTINUE
C-------------------------------
C     MESSAGES DESACTIVATION
C-------------------------------
      DO 210 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(1,il)==-2) THEN
         in1=itab(ksi(1,il))
         in2=itab(ksi(2,il))
         in3=itab(ksi(3,il))
         in4=itab(ksi(4,il))
#include "lockon.inc"
         WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(istdo,'(A,A,A,/,2I8,A,4I8,A)')' 3NODES FACET',
     .         ' FROM 4NODES ELEMENT/SEGMENT,',                 
     .         ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
     .         in1,in2,'ISOCENTER(',in1,in2,in3,in4,')'
         WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(iout,'(A,A,A,/,2I8,A,4I8,A)') ' 3NODES FACET',
     .         ' FROM 4NODES ELEMENT/SEGMENT,',                 
     .         ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
     .         in1,in2,'ISOCENTER(',in1,in2,in3,in4,')'
#include "lockoff.inc"
         iactiv(1,il)=-1
      ENDIF
 210  CONTINUE
C-------------------------------
C     2EME TRIANGLE.
C-------------------------------
      DO 225 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      x1(i)=gx(1,i)
      y1(i)=gx(2,i)
      z1(i)=gx(3,i)
      x2(i)=xp2(1,i)
      y2(i)=xp2(2,i)
      z2(i)=xp2(3,i)
      x3(i)=xp3(1,i)
      y3(i)=xp3(2,i)
      z3(i)=xp3(3,i)
      x12(i)=x2(i)-x1(i)
      y12(i)=y2(i)-y1(i)
      z12(i)=z2(i)-z1(i)
      x13(i)=x3(i)-x1(i)
      y13(i)=y3(i)-y1(i)
      z13(i)=z3(i)-z1(i)
      n1=y12(i)*z13(i)-z12(i)*y13(i)
      n2=z12(i)*x13(i)-x12(i)*z13(i)
      n3=x12(i)*y13(i)-y12(i)*x13(i)
      ntn=one/max(em20,sqrt(n1*n1+n2*n2+n3*n3))
      ntx(2,i)=ntn*n1
      nty(2,i)=ntn*n2
      ntz(2,i)=ntn*n3
      d(i)    =-ntx(2,i)*x1(i)-nty(2,i)*y1(i)-ntz(2,i)*z1(i)
C-----
      x23(i)=x3(i)-x2(i)
      y23(i)=y3(i)-y2(i)
      z23(i)=z3(i)-z2(i)
 225  CONTINUE
C-------------------------------
C     POINTS K,H SUR E,L / LA DISTANCE D(K,H) EST LOCALEMENT MAXIMUM
C-------------------------------
      DO 250 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      eh  =(abs(ntx(2,i)/(dgr*an))**expn)*an
     .    +(abs(nty(2,i)/(dgr*bn))**expn)*bn
     .    +(abs(ntz(2,i)/(dgr*cn))**expn)*cn
C     X EST DU SIGNE DE LAMBDA*NTX, IDEM EN Y ET Z.
C     LAMBDA1=EXP(LOG(1./EH)/EXPN)
      lambda1=(max(em20,eh))**(-one/expn)
      xh1   =abs(lambda1*ntx(2,i)/(dgr*an))**(one/(dgr-one))
      IF (ntx(2,i)<zero) xh1=-xh1
      yh1   =abs(lambda1*nty(2,i)/(dgr*bn))**(one/(dgr-one))
      IF (nty(2,i)<zero) yh1=-yh1
      zh1   =abs(lambda1*ntz(2,i)/(dgr*cn))**(one/(dgr-one))
      IF (ntz(2,i)<zero) zh1=-zh1
      mu1   =-ntx(2,i)*xh1-nty(2,i)*yh1-ntz(2,i)*zh1-d(i)
C     LAMBDA2=-LAMBDA1
      xh2   =-xh1
      yh2   =-yh1
      zh2   =-zh1
C      MU2   =-NTX(2,I)*XH2-NTY(2,I)*YH2-NTZ(2,I)*ZH2-D(I)
      mu2=-mu1-two*d(i)
      xtk(2,i) =xh1+mu1*ntx(2,i)
      ytk(2,i) =yh1+mu1*nty(2,i)
      ztk(2,i) =zh1+mu1*ntz(2,i)
      xtk2(i)  =xh2+mu2*ntx(2,i)
      ytk2(i)  =yh2+mu2*nty(2,i)
      ztk2(i)  =zh2+mu2*ntz(2,i)
C-------------------------------
 250  CONTINUE
C-------------------------------
C     RAMENER LE PT PK SUR LE TRIANGLE.
C-------------------------------
      DO 275 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      dx1=xtk(2,i)-x1(i)
      dy1=ytk(2,i)-y1(i)
      dz1=ztk(2,i)-z1(i)
      dx2=xtk(2,i)-x2(i)
      dy2=ytk(2,i)-y2(i)
      dz2=ztk(2,i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(2,i)
     .     +(dz1*dx2-dz2*dx1)*nty(2,i)
     .     +(dx1*dy2-dy1*dx2)*ntz(2,i)
      IF (out<zero) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(2,i)=alp*x1(i)+bet*x2(i)
        ytk(2,i)=alp*y1(i)+bet*y2(i)
        ztk(2,i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk(2,i)-x3(i)
      dy3=ytk(2,i)-y3(i)
      dz3=ztk(2,i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(2,i)
     .     +(dz2*dx3-dz3*dx2)*nty(2,i)
     .     +(dx2*dy3-dy2*dx3)*ntz(2,i)
      IF (out<zero) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk(2,i)=alp*x2(i)+bet*x3(i)
        ytk(2,i)=alp*y2(i)+bet*y3(i)
        ztk(2,i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(2,i)
     .     +(dz3*dx1-dz1*dx3)*nty(2,i)
     .     +(dx3*dy1-dy3*dx1)*ntz(2,i)
      IF (out<zero) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(2,i)=alp*x3(i)+bet*x1(i)
        ytk(2,i)=alp*y3(i)+bet*y1(i)
        ztk(2,i)=alp*z3(i)+bet*z1(i)
      ENDIF 
C-----
      dx1=xtk2(i)-x1(i)
      dy1=ytk2(i)-y1(i)
      dz1=ztk2(i)-z1(i)
      dx2=xtk2(i)-x2(i)
      dy2=ytk2(i)-y2(i)
      dz2=ztk2(i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(2,i)
     .      +(dz1*dx2-dz2*dx1)*nty(2,i)
     .      +(dx1*dy2-dy1*dx2)*ntz(2,i)
      IF (out<zero) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk2(i)=alp*x1(i)+bet*x2(i)
        ytk2(i)=alp*y1(i)+bet*y2(i)
        ztk2(i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk2(i)-x3(i)
      dy3=ytk2(i)-y3(i)
      dz3=ztk2(i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(2,i)
     .      +(dz2*dx3-dz3*dx2)*nty(2,i)
     .      +(dx2*dy3-dy2*dx3)*ntz(2,i)
      IF (out<zero) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk2(i)=alp*x2(i)+bet*x3(i)
        ytk2(i)=alp*y2(i)+bet*y3(i)
        ztk2(i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(2,i)
     .      +(dz3*dx1-dz1*dx3)*nty(2,i)
     .      +(dx3*dy1-dy3*dx1)*ntz(2,i)
      IF (out<zero) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk2(i)=alp*x3(i)+bet*x1(i)
        ytk2(i)=alp*y3(i)+bet*y1(i)
        ztk2(i)=alp*z3(i)+bet*z1(i)
      ENDIF 
C-------------------------------
 275  CONTINUE
C-------------------------------
C     PROJECTION DE PK SUR L ET PENETRATION.
C-------------------------------
      DO 300 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
      IF (iactiv(2,il)==-1) GOTO 300
C-----
      xh =xtk(2,i)
      yh =ytk(2,i)
      zh =ztk(2,i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one
      ykn1 =yh**(idg-1)
      sgnykn=-1.
      IF (ykn1*yh>=0.) sgnykn=1.
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one
      n11 =sgnxkn*xkn1*an
      n12 =sgnykn*ykn1*bn
      n13 =sgnzkn*zkn1*cn
      nr1 =n11*n11+n12*n12+n13*n13
      nr1 =sqrt(nr1)
      em =n11*xtk(2,i)+n12*ytk(2,i)+n13*ztk(2,i)
      IF (em<=one) THEN
        lambda1=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr1)
        inside1=1
      ELSE
        inside1=0
      ENDIF
C-----
      xh =xtk2(i)
      yh =ytk2(i)
      zh =ztk2(i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one
      ykn1 =yh**(idg-1)
      sgnykn=-one
      IF (ykn1*yh>=zero) sgnykn=one     
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one
      n21 =sgnxkn*xkn1*an
      n22 =sgnykn*ykn1*bn
      n23 =sgnzkn*zkn1*cn
      nr2  =n21*n21+n22*n22+n23*n23
      nr2 =sqrt(nr2)
      em =n21*xtk2(i)+n22*ytk2(i)+n23*ztk2(i)
      IF (em<=one) THEN
        lambda2=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr2)
        inside2=1
      ELSE
        inside2=0
      ENDIF
C-----
      IF (inside1==0.AND.inside2==0) THEN
        iactiv(2,il)=0
      ELSE
C-----
      IF (iactiv(2,il)==0) THEN
        IF (inside1/=0.AND.inside2/=0) THEN
         IF (abs(lambda1)>=abs(lambda2)) THEN
          xm(i)=xtk(2,i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(2,i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(2,i)-lambda1*n13/max(em20,nr1)
         ELSE
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(2,i)=xtk2(i)
          ytk(2,i)=ytk2(i)
          ztk(2,i)=ztk2(i)
         ENDIF
        ELSEIF(inside1/=0) THEN
          xm(i)=xtk(2,i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(2,i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(2,i)-lambda1*n13/max(em20,nr1)
        ELSEIF(inside2/=0) THEN
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(2,i)=xtk2(i)
          ytk(2,i)=ytk2(i)
          ztk(2,i)=ztk2(i)
        ENDIF
C-----
      ELSE
        xh=hold(1,4*(il-1)+2)
        yh=hold(2,4*(il-1)+2)
        zh=hold(3,4*(il-1)+2)
        n1=nold(1,4*(il-1)+2)
        n2=nold(2,4*(il-1)+2)
        n3=nold(3,4*(il-1)+2)
        lambda1=(xh-xtk(2,i))*n1
     .         +(yh-ytk(2,i))*n2
     .         +(zh-ztk(2,i))*n3
        lambda2=(xh-xtk2(i))*n1
     .         +(yh-ytk2(i))*n2
     .         +(zh-ztk2(i))*n3
        IF (inside1/=0.AND.inside2/=0) THEN
         IF (abs(lambda1)>=abs(lambda2)) THEN
           xm(i)=xtk(2,i)+lambda1*n1
           ym(i)=ytk(2,i)+lambda1*n2
           zm(i)=ztk(2,i)+lambda1*n3
          ELSE
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(2,i)=xtk2(i)
           ytk(2,i)=ytk2(i)
           ztk(2,i)=ztk2(i)
          ENDIF
        ELSEIF(inside1/=0) THEN
           xm(i)=xtk(2,i)+lambda1*n1
           ym(i)=ytk(2,i)+lambda1*n2
           zm(i)=ztk(2,i)+lambda1*n3
        ELSEIF(inside2/=0) THEN
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(2,i)=xtk2(i)
           ytk(2,i)=ytk2(i)
           ztk(2,i)=ztk2(i)
        ENDIF
C       one more iteration.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        em=n1*xm(i)+n2*ym(i)+n3*zm(i)
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))        
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one        
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        nr =n1*n1+n2*n2+n3*n3
        nr =one/max(em20,sqrt(nr))
        n1 =n1*nr
        n2 =n2*nr
        n3 =n3*nr
        lambda1=(xm(i)-xtk(2,i))*n1
     .         +(ym(i)-ytk(2,i))*n2
     .         +(zm(i)-ztk(2,i))*n3
        xm(i)=xtk(2,i)+lambda1*n1
        ym(i)=ytk(2,i)+lambda1*n2
        zm(i)=ztk(2,i)+lambda1*n3
      ENDIF 
      iactiv(2,il)=iactiv(2,il)+1
      ENDIF
 300  CONTINUE
C-----
#include "vectorize.inc"
      DO 325 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(2,il)<=0) GOTO 325
C-----
C      projection radiale de M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one        
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       em=n1*xm(i)+n2*ym(i)+n3*zm(i)
C------
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))        
C------
C      normale a l'ellipsoide en M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
       zkn1 =zm(i)**(idg-1)
       sgnzkn=-one
       IF (zkn1*zm(i)>=zero) sgnzkn=one        
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       nr =n1*n1+n2*n2+n3*n3
       nr =one/max(em20,sqrt(nr))
       nxi(2,i)=n1*nr
       nyi(2,i)=n2*nr
       nzi(2,i)=n3*nr
       xi(2,i)=xm(i)
       yi(2,i)=ym(i)
       zi(2,i)=zm(i)
       depth(2,i)=xm(i)*nxi(2,i)+ym(i)*nyi(2,i)+zm(i)*nzi(2,i)
C------
       penet(2,i)=(xtk(2,i)-xm(i))*nxi(2,i)
     .          +(ytk(2,i)-ym(i))*nyi(2,i)
     .          +(ztk(2,i)-zm(i))*nzi(2,i)
       penet(2,i)=-penet(2,i)
       penet(2,i)=max(zero,penet(2,i))
       IF (depth(2,i)-penet(2,i)<em10*depth(2,i)) THEN
         penet(2,i)=zero
         iactiv(2,il)=-2
       ENDIF
       ansmx=max(ansmx,penet(2,i))
 325  CONTINUE
C-------------------------------
C     MESSAGES DESACTIVATION
C-------------------------------
      DO 335 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(2,il)==-2) THEN
         in1=itab(ksi(1,il))
         in2=itab(ksi(2,il))
         in3=itab(ksi(3,il))
         in4=itab(ksi(4,il))
#include "lockon.inc"
         WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(istdo,'(A,A,A,/,2I8,A,4I8,A)')' 3NODES FACET',
     .         ' FROM 4NODES ELEMENT/SEGMENT,',                 
     .         ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
     .         in2,in3,'ISOCENTER(',in1,in2,in3,in4,')'
         WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(iout,'(A,A,A,/,2I8,A,4I8,A)') ' 3NODES FACET',
     .         ' FROM 4NODES ELEMENT/SEGMENT,',                 
     .         ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
     .         in2,in3,'ISOCENTER(',in1,in2,in3,in4,')'
#include "lockoff.inc"
         iactiv(2,il)=-1
      ENDIF
 335  CONTINUE
C-------------------------------
C     3EME TRIANGLE.
C-------------------------------
      DO 350 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      x1(i)=gx(1,i)
      y1(i)=gx(2,i)
      z1(i)=gx(3,i)
      x2(i)=xp3(1,i)
      y2(i)=xp3(2,i)
      z2(i)=xp3(3,i)
      x3(i)=xp4(1,i)
      y3(i)=xp4(2,i)
      z3(i)=xp4(3,i)
      x12(i)=x2(i)-x1(i)
      y12(i)=y2(i)-y1(i)
      z12(i)=z2(i)-z1(i)
      x13(i)=x3(i)-x1(i)
      y13(i)=y3(i)-y1(i)
      z13(i)=z3(i)-z1(i)
      n1=y12(i)*z13(i)-z12(i)*y13(i)
      n2=z12(i)*x13(i)-x12(i)*z13(i)
      n3=x12(i)*y13(i)-y12(i)*x13(i)
      ntn=one/max(em20,sqrt(n1*n1+n2*n2+n3*n3))
      ntx(3,i)=ntn*n1
      nty(3,i)=ntn*n2
      ntz(3,i)=ntn*n3
      d(i)    =-ntx(3,i)*x1(i)-nty(3,i)*y1(i)-ntz(3,i)*z1(i)
C-----
      x23(i)=x3(i)-x2(i)
      y23(i)=y3(i)-y2(i)
      z23(i)=z3(i)-z2(i)
 350  CONTINUE
C-------------------------------
C     POINTS K,H SUR E,L / LA DISTANCE D(K,H) EST LOCALEMENT MAXIMUM
C-------------------------------
      DO 375 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      eh  =(abs(ntx(3,i)/(dgr*an))**expn)*an
     .    +(abs(nty(3,i)/(dgr*bn))**expn)*bn
     .    +(abs(ntz(3,i)/(dgr*cn))**expn)*cn
C     X EST DU SIGNE DE LAMBDA*NTX, IDEM EN Y ET Z.
C     LAMBDA1=EXP(LOG(1./EH)/EXPN)
      lambda1=(max(em20,eh))**(-one/expn)
      xh1   =abs(lambda1*ntx(3,i)/(dgr*an))**(one/(dgr-one))
      IF (ntx(3,i)<zero) xh1=-xh1
      yh1   =abs(lambda1*nty(3,i)/(dgr*bn))**(one/(dgr-one))
      IF (nty(3,i)<zero) yh1=-yh1
      zh1   =abs(lambda1*ntz(3,i)/(dgr*cn))**(one/(dgr-one))
      IF (ntz(3,i)<zero) zh1=-zh1
      mu1   =-ntx(3,i)*xh1-nty(3,i)*yh1-ntz(3,i)*zh1-d(i)
C     LAMBDA2=-LAMBDA1
      xh2   =-xh1
      yh2   =-yh1
      zh2   =-zh1
C      MU2   =-NTX(3,I)*XH2-NTY(3,I)*YH2-NTZ(3,I)*ZH2-D(I)
      mu2=-mu1-two*d(i)
      xtk(3,i) =xh1+mu1*ntx(3,i)
      ytk(3,i) =yh1+mu1*nty(3,i)
      ztk(3,i) =zh1+mu1*ntz(3,i)
      xtk2(i)=xh2+mu2*ntx(3,i)
      ytk2(i)=yh2+mu2*nty(3,i)
      ztk2(i)=zh2+mu2*ntz(3,i)
C-------------------------------
 375  CONTINUE
C-------------------------------
C     RAMENER LE PT PK SUR LE TRIANGLE.
C-------------------------------
      DO 400 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      dx1=xtk(3,i)-x1(i)
      dy1=ytk(3,i)-y1(i)
      dz1=ztk(3,i)-z1(i)
      dx2=xtk(3,i)-x2(i)
      dy2=ytk(3,i)-y2(i)
      dz2=ztk(3,i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(3,i)
     .     +(dz1*dx2-dz2*dx1)*nty(3,i)
     .     +(dx1*dy2-dy1*dx2)*ntz(3,i)
      IF (out<zero) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(3,i)=alp*x1(i)+bet*x2(i)
        ytk(3,i)=alp*y1(i)+bet*y2(i)
        ztk(3,i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk(3,i)-x3(i)
      dy3=ytk(3,i)-y3(i)
      dz3=ztk(3,i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(3,i)
     .     +(dz2*dx3-dz3*dx2)*nty(3,i)
     .     +(dx2*dy3-dy2*dx3)*ntz(3,i)
      IF (out<zero) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(3,i)=alp*x2(i)+bet*x3(i)
        ytk(3,i)=alp*y2(i)+bet*y3(i)
        ztk(3,i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(3,i)
     .     +(dz3*dx1-dz1*dx3)*nty(3,i)
     .     +(dx3*dy1-dy3*dx1)*ntz(3,i)
      IF (out<0.) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(3,i)=alp*x3(i)+bet*x1(i)
        ytk(3,i)=alp*y3(i)+bet*y1(i)
        ztk(3,i)=alp*z3(i)+bet*z1(i)
      ENDIF 
C-----
      dx1=xtk2(i)-x1(i)
      dy1=ytk2(i)-y1(i)
      dz1=ztk2(i)-z1(i)
      dx2=xtk2(i)-x2(i)
      dy2=ytk2(i)-y2(i)
      dz2=ztk2(i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(3,i)
     .      +(dz1*dx2-dz2*dx1)*nty(3,i)
     .      +(dx1*dy2-dy1*dx2)*ntz(3,i)
      IF (out<zero) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk2(i)=alp*x1(i)+bet*x2(i)
        ytk2(i)=alp*y1(i)+bet*y2(i)
        ztk2(i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk2(i)-x3(i)
      dy3=ytk2(i)-y3(i)
      dz3=ztk2(i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(3,i)
     .      +(dz2*dx3-dz3*dx2)*nty(3,i)
     .      +(dx2*dy3-dy2*dx3)*ntz(3,i)
      IF (out<zero) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk2(i)=alp*x2(i)+bet*x3(i)
        ytk2(i)=alp*y2(i)+bet*y3(i)
        ztk2(i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(3,i)
     .      +(dz3*dx1-dz1*dx3)*nty(3,i)
     .      +(dx3*dy1-dy3*dx1)*ntz(3,i)
      IF (out<zero) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk2(i)=alp*x3(i)+bet*x1(i)
        ytk2(i)=alp*y3(i)+bet*y1(i)
        ztk2(i)=alp*z3(i)+bet*z1(i)
      ENDIF
C-------------------------------
 400  CONTINUE
C-------------------------------
C     PROJECTION DE PK SUR L ET PENETRATION.
C-------------------------------
      DO 425 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
      IF (iactiv(3,il)==-1) GOTO 425
C-----
      xh =xtk(3,i)
      yh =ytk(3,i)
      zh =ztk(3,i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one      
      ykn1 =yh**(idg-1)
      sgnykn=-one
      IF (ykn1*yh>=zero) sgnykn=one     
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one      
      n11 =sgnxkn*xkn1*an
      n12 =sgnykn*ykn1*bn
      n13 =sgnzkn*zkn1*cn
      nr1 =n11*n11+n12*n12+n13*n13
      nr1 =sqrt(nr1)
      em =n11*xtk(3,i)+n12*ytk(3,i)+n13*ztk(3,i)
      IF (em<=one) THEN
        lambda1=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr1)
        inside1=1
      ELSE
        inside1=0
      ENDIF
C-----
      xh =xtk2(i)
      yh =ytk2(i)
      zh =ztk2(i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one      
      ykn1 =yh**(idg-1)
      sgnykn=-one
      IF (ykn1*yh>=zero) sgnykn=one     
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one      
      n21 =sgnxkn*xkn1*an
      n22 =sgnykn*ykn1*bn
      n23 =sgnzkn*zkn1*cn
      nr2 =n21*n21+n22*n22+n23*n23
      nr2 =sqrt(nr2)
      em =n21*xtk2(i)+n22*ytk2(i)+n23*ztk2(i)
       IF (em<=one) THEN
        lambda2=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr2)    
        inside2=1
      ELSE
        inside2=0
      ENDIF
C-----
      IF (inside1==0.AND.inside2==0) THEN
        iactiv(3,il)=0
      ELSE
C-----
      IF (iactiv(3,il)==0) THEN
        IF (inside1/=0.AND.inside2/=0) THEN
         IF (abs(lambda1)>=abs(lambda2)) THEN
          xm(i)=xtk(3,i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(3,i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(3,i)-lambda1*n13/max(em20,nr1)
         ELSE
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(3,i)=xtk2(i)
          ytk(3,i)=ytk2(i)
          ztk(3,i)=ztk2(i)
         ENDIF
        ELSEIF(inside1/=0) THEN
          xm(i)=xtk(3,i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(3,i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(3,i)-lambda1*n13/max(em20,nr1)
        ELSEIF(inside2/=0) THEN
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(3,i)=xtk2(i)
          ytk(3,i)=ytk2(i)
          ztk(3,i)=ztk2(i)
        ENDIF
C-----
      ELSE
        xh=hold(1,4*(il-1)+3)
        yh=hold(2,4*(il-1)+3)
        zh=hold(3,4*(il-1)+3)
        n1=nold(1,4*(il-1)+3)
        n2=nold(2,4*(il-1)+3)
        n3=nold(3,4*(il-1)+3)
        lambda1=(xh-xtk(3,i))*n1
     .         +(yh-ytk(3,i))*n2
     .         +(zh-ztk(3,i))*n3
        lambda2=(xh-xtk2(i))*n1
     .         +(yh-ytk2(i))*n2
     .         +(zh-ztk2(i))*n3
        IF (inside1/=0.AND.inside2/=0) THEN
         IF (abs(lambda1)>=abs(lambda2)) THEN
           xm(i)=xtk(3,i)+lambda1*n1
           ym(i)=ytk(3,i)+lambda1*n2
           zm(i)=ztk(3,i)+lambda1*n3
          ELSE
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(3,i)=xtk2(i)
           ytk(3,i)=ytk2(i)
           ztk(3,i)=ztk2(i)
          ENDIF
        ELSEIF(inside1/=0) THEN
           xm(i)=xtk(3,i)+lambda1*n1
           ym(i)=ytk(3,i)+lambda1*n2
           zm(i)=ztk(3,i)+lambda1*n3
        ELSEIF(inside2/=0) THEN
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(3,i)=xtk2(i)
           ytk(3,i)=ytk2(i)
           ztk(3,i)=ztk2(i)
        ENDIF
C       one more iteration.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        em=n1*xm(i)+n2*ym(i)+n3*zm(i)
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))        
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one        
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        nr =n1*n1+n2*n2+n3*n3
        nr =one/max(em20,sqrt(nr))
        n1 =n1*nr
        n2 =n2*nr
        n3 =n3*nr
        lambda1=(xm(i)-xtk(3,i))*n1
     .         +(ym(i)-ytk(3,i))*n2
     .         +(zm(i)-ztk(3,i))*n3
        xm(i)=xtk(3,i)+lambda1*n1
        ym(i)=ytk(3,i)+lambda1*n2
        zm(i)=ztk(3,i)+lambda1*n3
      ENDIF 
      iactiv(3,il)=iactiv(3,il)+1
      ENDIF
 425  CONTINUE
C-----
#include "vectorize.inc"
      DO 450 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(3,il)<=0) GOTO 450
C-----
C      projection radiale de M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one        
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       em=n1*xm(i)+n2*ym(i)+n3*zm(i)
C------
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))        
C------
C      normale a l'ellipsoide en M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one        
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       nr =n1*n1+n2*n2+n3*n3
       nr =one/max(em20,sqrt(nr))
       nxi(3,i)=n1*nr
       nyi(3,i)=n2*nr
       nzi(3,i)=n3*nr
       xi(3,i)=xm(i)
       yi(3,i)=ym(i)
       zi(3,i)=zm(i)
       depth(3,i)=xm(i)*nxi(3,i)+ym(i)*nyi(3,i)+zm(i)*nzi(3,i)
C------
       penet(3,i)=(xtk(3,i)-xm(i))*nxi(3,i)
     .          +(ytk(3,i)-ym(i))*nyi(3,i)
     .          +(ztk(3,i)-zm(i))*nzi(3,i)
       penet(3,i)=-penet(3,i)
       penet(3,i)=max(zero,penet(3,i))
       IF (depth(3,i)-penet(3,i)<em10*depth(3,i)) THEN
         penet(3,i)=zero
         iactiv(3,il)=-2
       ENDIF
       ansmx=max(ansmx,penet(3,i))
 450  CONTINUE
C-------------------------------
C     MESSAGES DESACTIVATION
C-------------------------------
      DO 460 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(3,il)==-2) THEN
         in1=itab(ksi(1,il))
         in2=itab(ksi(2,il))
         in3=itab(ksi(3,il))
         in4=itab(ksi(4,il))
#include "lockon.inc"
         WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(istdo,'(A,A,A,/,2I8,A,4I8,A)')' 3NODES FACET',
     .         ' FROM 4NODES ELEMENT/SEGMENT,',                 
     .         ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
     .         in3,in4,'ISOCENTER(',in1,in2,in3,in4,')'
         WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(iout,'(A,A,A,/,2I8,A,4I8,A)') ' 3NODES FACET',
     .         ' FROM 4NODES ELEMENT/SEGMENT,',                 
     .         ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
     .         in3,in4,'ISOCENTER(',in1,in2,in3,in4,')'
#include "lockoff.inc"
         iactiv(3,il)=-1
      ENDIF
 460  CONTINUE
C-------------------------------
C     4EME TRIANGLE.
C-------------------------------
      DO 475 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      x1(i)=gx(1,i)
      y1(i)=gx(2,i)
      z1(i)=gx(3,i)
      x2(i)=xp4(1,i)
      y2(i)=xp4(2,i)
      z2(i)=xp4(3,i)
      x3(i)=xp1(1,i)
      y3(i)=xp1(2,i)
      z3(i)=xp1(3,i)
      x12(i)=x2(i)-x1(i)
      y12(i)=y2(i)-y1(i)
      z12(i)=z2(i)-z1(i)
      x13(i)=x3(i)-x1(i)
      y13(i)=y3(i)-y1(i)
      z13(i)=z3(i)-z1(i)
      n1=y12(i)*z13(i)-z12(i)*y13(i)
      n2=z12(i)*x13(i)-x12(i)*z13(i)
      n3=x12(i)*y13(i)-y12(i)*x13(i)
      ntn=one/max(em20,sqrt(n1*n1+n2*n2+n3*n3))
      ntx(4,i)=ntn*n1
      nty(4,i)=ntn*n2
      ntz(4,i)=ntn*n3
      d(i)    =-ntx(4,i)*x1(i)-nty(4,i)*y1(i)-ntz(4,i)*z1(i)
C-----
      x23(i)=x3(i)-x2(i)
      y23(i)=y3(i)-y2(i)
      z23(i)=z3(i)-z2(i)
 475  CONTINUE
C-------------------------------
C     POINTS K,H SUR E,L / LA DISTANCE D(K,H) EST LOCALEMENT MAXIMUM
C-------------------------------
      DO 500 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      eh  =(abs(ntx(4,i)/(dgr*an))**expn)*an
     .    +(abs(nty(4,i)/(dgr*bn))**expn)*bn
     .    +(abs(ntz(4,i)/(dgr*cn))**expn)*cn
C     X EST DU SIGNE DE LAMBDA*NTX, IDEM EN Y ET Z.
C     LAMBDA1=EXP(LOG(1./EH)/EXPN)
      lambda1=(max(em20,eh))**(-one/expn)
      xh1   =abs(lambda1*ntx(4,i)/(dgr*an))**(one/(dgr-one))
      IF (ntx(4,i)<zero) xh1=-xh1
      yh1   =abs(lambda1*nty(4,i)/(dgr*bn))**(one/(dgr-one))
      IF (nty(4,i)<zero) yh1=-yh1
      zh1   =abs(lambda1*ntz(4,i)/(dgr*cn))**(one/(dgr-one))
      IF (ntz(4,i)<zero) zh1=-zh1
      mu1   =-ntx(4,i)*xh1-nty(4,i)*yh1-ntz(4,i)*zh1-d(i)
C     LAMBDA2=-LAMBDA1
      xh2   =-xh1
      yh2   =-yh1
      zh2   =-zh1
C      MU2   =-NTX(4,I)*XH2-NTY(4,I)*YH2-NTZ(4,I)*ZH2-D(I)
      mu2=-mu1-two*d(i)
      xtk(4,i) =xh1+mu1*ntx(4,i)
      ytk(4,i) =yh1+mu1*nty(4,i)
      ztk(4,i) =zh1+mu1*ntz(4,i)
      xtk2(i)=xh2+mu2*ntx(4,i)
      ytk2(i)=yh2+mu2*nty(4,i)
      ztk2(i)=zh2+mu2*ntz(4,i)
C-------------------------------
 500  CONTINUE
C-------------------------------
C     RAMENER LE PT PK SUR LE TRIANGLE.
C-------------------------------
      DO 525 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      dx1=xtk(4,i)-x1(i)
      dy1=ytk(4,i)-y1(i)
      dz1=ztk(4,i)-z1(i)
      dx2=xtk(4,i)-x2(i)
      dy2=ytk(4,i)-y2(i)
      dz2=ztk(4,i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(4,i)
     .     +(dz1*dx2-dz2*dx1)*nty(4,i)
     .     +(dx1*dy2-dy1*dx2)*ntz(4,i)
      IF (out<zero) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(4,i)=alp*x1(i)+bet*x2(i)
        ytk(4,i)=alp*y1(i)+bet*y2(i)
        ztk(4,i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk(4,i)-x3(i)
      dy3=ytk(4,i)-y3(i)
      dz3=ztk(4,i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(4,i)
     .     +(dz2*dx3-dz3*dx2)*nty(4,i)
     .     +(dx2*dy3-dy2*dx3)*ntz(4,i)
      IF (out<0.) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(4,i)=alp*x2(i)+bet*x3(i)
        ytk(4,i)=alp*y2(i)+bet*y3(i)
        ztk(4,i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(4,i)
     .     +(dz3*dx1-dz1*dx3)*nty(4,i)
     .     +(dx3*dy1-dy3*dx1)*ntz(4,i)
      IF (out<0.) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(4,i)=alp*x3(i)+bet*x1(i)
        ytk(4,i)=alp*y3(i)+bet*y1(i)
        ztk(4,i)=alp*z3(i)+bet*z1(i)
      ENDIF 
C-----
      dx1=xtk2(i)-x1(i)
      dy1=ytk2(i)-y1(i)
      dz1=ztk2(i)-z1(i)
      dx2=xtk2(i)-x2(i)
      dy2=ytk2(i)-y2(i)
      dz2=ztk2(i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(4,i)
     .      +(dz1*dx2-dz2*dx1)*nty(4,i)
     .      +(dx1*dy2-dy1*dx2)*ntz(4,i)
      IF (out<zero) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk2(i)=alp*x1(i)+bet*x2(i)
        ytk2(i)=alp*y1(i)+bet*y2(i)
        ztk2(i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk2(i)-x3(i)
      dy3=ytk2(i)-y3(i)
      dz3=ztk2(i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(4,i)
     .      +(dz2*dx3-dz3*dx2)*nty(4,i)
     .      +(dx2*dy3-dy2*dx3)*ntz(4,i)
      IF (out<zero) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk2(i)=alp*x2(i)+bet*x3(i)
        ytk2(i)=alp*y2(i)+bet*y3(i)
        ztk2(i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(4,i)
     .      +(dz3*dx1-dz1*dx3)*nty(4,i)
     .      +(dx3*dy1-dy3*dx1)*ntz(4,i)
      IF (out<zero) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=one-bet
        xtk2(i)=alp*x3(i)+bet*x1(i)
        ytk2(i)=alp*y3(i)+bet*y1(i)
        ztk2(i)=alp*z3(i)+bet*z1(i)
      ENDIF 
C-------------------------------
 525  CONTINUE
C-------------------------------
C     PROJECTION DE PK SUR L ET PENETRATION.
C-------------------------------
      DO 550 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
      IF (iactiv(4,il)==-1) GOTO 550
C-----
      xh =xtk(4,i)
      yh =ytk(4,i)
      zh =ztk(4,i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one
      ykn1 =yh**(idg-1)
      sgnykn=-one
      IF (ykn1*yh>=zero) sgnykn=one
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one
      n11 =sgnxkn*xkn1*an
      n12 =sgnykn*ykn1*bn
      n13 =sgnzkn*zkn1*cn
      nr1 =n11*n11+n12*n12+n13*n13
      nr1 =sqrt(nr1)
      em =n11*xtk(4,i)+n12*ytk(4,i)+n13*ztk(4,i)
      IF (em<=one) THEN
        lambda1=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr1)
        inside1=1
      ELSE
        inside1=0
      ENDIF
C-----
      xh =xtk2(i)
      yh =ytk2(i)
      zh =ztk2(i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one
      ykn1 =yh**(idg-1)
      sgnykn=-one
      IF (ykn1*yh>=one) sgnykn=one
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one
      n21 =sgnxkn*xkn1*an
      n22 =sgnykn*ykn1*bn
      n23 =sgnzkn*zkn1*cn
      nr2 =n21*n21+n22*n22+n23*n23
      nr2 =sqrt(nr2)
      em =n21*xtk2(i)+n22*ytk2(i)+n23*ztk2(i)
      IF (em<=one) THEN
        lambda2=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr2)
        inside2=1
      ELSE
        inside2=0
      ENDIF
C-----
      IF (inside1==0.AND.inside2==0) THEN
        iactiv(4,il)=0
      ELSE
C-----
      IF (iactiv(4,il)==0) THEN
        IF (inside1/=0.AND.inside2/=0) THEN
         IF (abs(lambda1)>=abs(lambda2)) THEN
          xm(i)=xtk(4,i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(4,i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(4,i)-lambda1*n13/max(em20,nr1)
         ELSE
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(4,i)=xtk2(i)
          ytk(4,i)=ytk2(i)
          ztk(4,i)=ztk2(i)
         ENDIF
        ELSEIF(inside1/=0) THEN
          xm(i)=xtk(4,i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(4,i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(4,i)-lambda1*n13/max(em20,nr1)
        ELSEIF(inside2/=0) THEN
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(4,i)=xtk2(i)
          ytk(4,i)=ytk2(i)
          ztk(4,i)=ztk2(i)
        ENDIF
C-----
      ELSE
        xh=hold(1,4*(il-1)+4)
        yh=hold(2,4*(il-1)+4)
        zh=hold(3,4*(il-1)+4)
        n1=nold(1,4*(il-1)+4)
        n2=nold(2,4*(il-1)+4)
        n3=nold(3,4*(il-1)+4)
        lambda1=(xh-xtk(4,i))*n1
     .         +(yh-ytk(4,i))*n2
     .         +(zh-ztk(4,i))*n3
        lambda2=(xh-xtk2(i))*n1
     .         +(yh-ytk2(i))*n2
     .         +(zh-ztk2(i))*n3
        IF (inside1/=0.AND.inside2/=0) THEN
         IF (abs(lambda1)>=abs(lambda2)) THEN
           xm(i)=xtk(4,i)+lambda1*n1
           ym(i)=ytk(4,i)+lambda1*n2
           zm(i)=ztk(4,i)+lambda1*n3
          ELSE
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(4,i)=xtk2(i)
           ytk(4,i)=ytk2(i)
           ztk(4,i)=ztk2(i)
          ENDIF
        ELSEIF(inside1/=0) THEN
           xm(i)=xtk(4,i)+lambda1*n1
           ym(i)=ytk(4,i)+lambda1*n2
           zm(i)=ztk(4,i)+lambda1*n3
        ELSEIF(inside2/=0) THEN
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(4,i)=xtk2(i)
           ytk(4,i)=ytk2(i)
           ztk(4,i)=ztk2(i)
        ENDIF
C       one more iteration.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        em=n1*xm(i)+n2*ym(i)+n3*zm(i)
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one 
        IF (ykn1*ym(i)>=zero) sgnykn=one
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        nr =n1*n1+n2*n2+n3*n3
        nr =1./max(em20,sqrt(nr))
        n1 =n1*nr
        n2 =n2*nr
        n3 =n3*nr
        lambda1=(xm(i)-xtk(4,i))*n1
     .         +(ym(i)-ytk(4,i))*n2
     .         +(zm(i)-ztk(4,i))*n3
        xm(i)=xtk(4,i)+lambda1*n1
        ym(i)=ytk(4,i)+lambda1*n2
        zm(i)=ztk(4,i)+lambda1*n3
      ENDIF 
      iactiv(4,il)=iactiv(4,il)+1
      ENDIF
 550  CONTINUE
C-----
#include "vectorize.inc"
      DO 575 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(4,il)<=0) GOTO 575
C-----
C      projection radiale de M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       em=n1*xm(i)+n2*ym(i)+n3*zm(i)
C------
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))
C------
C      normale a l'ellipsoide en M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       nr =n1*n1+n2*n2+n3*n3
       nr =one/max(em20,sqrt(nr))
       nxi(4,i)=n1*nr
       nyi(4,i)=n2*nr
       nzi(4,i)=n3*nr
       xi(4,i)=xm(i)
       yi(4,i)=ym(i)
       zi(4,i)=zm(i)
       depth(4,i)=xm(i)*nxi(4,i)+ym(i)*nyi(4,i)+zm(i)*nzi(4,i)
C------
       penet(4,i)=(xtk(4,i)-xm(i))*nxi(4,i)
     .          +(ytk(4,i)-ym(i))*nyi(4,i)
     .          +(ztk(4,i)-zm(i))*nzi(4,i)
       penet(4,i)=-penet(4,i)
       penet(4,i)=max(zero,penet(4,i))
       IF (depth(4,i)-penet(4,i)<em10*depth(4,i)) THEN
         penet(4,i)=zero
         iactiv(4,il)=-2
       ENDIF
       ansmx=max(ansmx,penet(4,i))
 575  CONTINUE
C-------------------------------
C     MESSAGES DESACTIVATION
C-------------------------------
      DO 585 nls=ndeb+1,min(ndeb+mvsiz,nsc)
      il =ksc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(4,il)==-2) THEN
         in1=itab(ksi(1,il))
         in2=itab(ksi(2,il))
         in3=itab(ksi(3,il))
         in4=itab(ksi(4,il))
#include "lockon.inc"
         WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(istdo,'(A,A,A,/,2I8,A,4I8,A)')' 3NODES FACET',
     .         ' FROM 4NODES ELEMENT/SEGMENT,',                 
     .         ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
     .         in4,in1,'ISOCENTER(',in1,in2,in3,in4,')'
         WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(iout,'(a,a,a,/,2i8,a,4i8,a)') ' 3nodes facet',
     .         ' from 4nodes element/segment,',                 
     .         ' de-activated from interface; facet vertices:',
     .         IN4,IN1,'isocenter(',IN1,IN2,IN3,IN4,')'
#include "lockoff.inc"
         IACTIV(4,IL)=-1
      ENDIF
 585  CONTINUE
C------------------------------------------------------------
      RETURN
      END