i10lagm.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/.
!||====================================================================
!||    i10lagm       ../engine/source/interfaces/int16/i10lagm.F
!||--- called by ------------------------------------------------------
!||    i16main       ../engine/source/interfaces/int16/i16main.F
!||--- calls      -----------------------------------------------------
!||    i10lll        ../engine/source/interfaces/int16/i10lagm.F
!||--- uses       -----------------------------------------------------
!||    element_mod   ../common_source/modules/elements/element_mod.F90
!||====================================================================
      SUBROUTINE i10lagm(X    ,V      ,LLL     ,JLL   ,SLL   ,
     2                  XLL   ,CANDN  ,CANDE   ,I_STOK,IXS   ,
     3                  IXS10 ,IADLL  ,EMINX   ,NSV   ,NELEM ,
     4                  N_MUL_MX,ITASK  ,A     ,ITIED ,
     5                  NINT  ,NKMAX  ,COMNTAG )
      use element_mod , only : nixs
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.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      "task_c.inc"
#include      "com04_c.inc"
#include      "com08_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER I_STOK,N_MUL_MX,ITASK,ITIED,NINT,NKMAX ,
     .  LLL(*),JLL(*),SLL(*),CANDN(*),CANDE(*),COMNTAG(*),
     .  IXS(NIXS,*),IXS10(6,*),IADLL(*),NSV(*)  ,NELEM(*) 
C     REAL
      my_real
     .  x(3,*),v(3,*),xll(*),
     .  eminx(6,*),a(3,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, K, IE, IS, IC, III(MVSIZ, 7), LLT, NFT, LE, FIRST, LAST, 
     .        I10
      my_real
     .   XX(MVSIZ,7),YY(MVSIZ,7),ZZ(MVSIZ,7),
     . dist
C-----------------------------------------------
C
C
C       | M | Lt| | a    | M ao
C       |---+---| |    = |
C       | L | 0 | | lambda   | bo
C
C        [M] a + [L]t la = [M] ao
C        [L] a = bo
C
C         a = -[M]-1[L]t la + ao
C        [L][M]-1[L]t la  = [L] ao - bo
C
C on pose:
C        [H] = [L][M]-1[L]t
C        b  = [L] ao - bo
C
C        [H] lambda = b
C
C        a = ao - [M]-1[L]t la
C-----------------------------------------------
C
C        lambda              : LAMBDA(NC)
C        ao              : A(NUMNOD)
C        L               : XLL(NK,NC)
C        M               : MAS(NUMNOD)
C        [L][M]-1[L]t la : HLA(NC)
C        [L] ao - b      : B(NC)
C        [M]-1[L]t la    : LTLA(NUMNOD)
C
C        NC : number of contacts
C        NK: Number of node for contact (8+1.16+1.8+8.16+16)
C
C        IC : contact number (1,NC)
C        IK : local node number for a contact (1,NK)
C        I  : global node number (1,NUMNOD)
C
C        IADLL(NC)        : IAD = IADLL(IC)
C        LLL(NC*(7,21))  : I  = LLL(IAD+1,2...IADNEXT-1)
C-----------------------------------------------
C  evaluation of b:
C
C         Vs = Somme(Ni Vi)
C         Vs_ + dt As = Somme(Ni Vi_) + Somme(dt Ni Ai)
C         Somme(dt Ni Ai) - dt As =  Vs_ -Somme(Ni Vi_)
C         [L] = dt {N1,N2,..,N15,-1}
C         bo = [L] a = -[L]/dt v_
C         b  = [L] ao - bo
C         b  = [L] ao + [L]/dt v_ = [L] (v_ + ao dt)/dt
C-----------------------------------------------
C                 b  = [L] vo+/dt   +   vout
C-----------------------------------------------
C-----------------------------------------------------------------------
C     loop over contact candidates
C-----------------------------------------------------------------------
      first = 1 + i_stok * itask / nthread
      last = i_stok*(itask+1) / nthread
      llt = 0
      nft=llt+1
      DO ic=first,last
       le = cande(ic)
       ie = nelem(le)
       i10 = ie - numels8
C-----------------------------------------------------------------------
C      test si shell 16  
C-----------------------------------------------------------------------
       IF(i10.ge .1.AND.i10.le .numels10)THEN
        is = nsv(candn(ic))
        dist = -1.e30
        dist = max(eminx(1,le)-x(1,is)-dt2*(v(1,is)+dt12*a(1,is)),
     .             x(1,is)+dt2*(v(1,is)+dt12*a(1,is))-eminx(4,le),dist)
        dist = max(eminx(2,le)-x(2,is)-dt2*(v(2,is)+dt12*a(2,is)),
     .             x(2,is)+dt2*(v(2,is)+dt12*a(2,is))-eminx(5,le),dist)
        dist = max(eminx(3,le)-x(3,is)-dt2*(v(3,is)+dt12*a(3,is)),
     .             x(3,is)+dt2*(v(3,is)+dt12*a(3,is))-eminx(6,le),dist)
c        IF (DIST<0.) CANDN(I) = -CANDN(I)
C-----------------------------------------------------------------------
C       test if inside the box
C-----------------------------------------------------------------------
        IF(dist.le .zero)THEN
c
c      print *, "in la boite",XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
c
          llt = llt+1
          iii(llt,7)=is
          xx(llt,7)=x(1,is)
          yy(llt,7)=x(2,is)
          zz(llt,7)=x(3,is)
          iii(llt,1)=ixs(2,ie)
          iii(llt,2)=ixs(4,ie)
          iii(llt,3)=ixs(6,ie)
          iii(llt,4)=ixs(7,ie)
          DO k=1,6
            iii(llt,k+8)=ixs10(k,i10)
          ENDDO
          DO k=1,7
            i = iii(llt,k)
            xx(llt,k)=x(1,i)
            yy(llt,k)=x(2,i)
            zz(llt,k)=x(3,i)
          ENDDO
c
C-----------------------------------------------------------------------
C     calculation of [L] by packet of mvsiz
C-----------------------------------------------------------------------
          IF(llt==mvsiz-1)THEN
            CALL i10lll(
     1       llt ,lll       ,jll       ,sll       ,xll     ,v       ,
     2       xx  ,yy        ,zz        ,iii       ,iadll   ,
     3       n_mul_mx ,a    ,x         ,itied     ,nint    ,nkmax   ,
     4       comntag )
            nft=llt+1
            llt = 0
          ENDIF
        ELSE
c debug
          k=0
        ENDIF
       ENDIF
      ENDDO
C-----------------------------------------------------------------------
C     calculation of [L] for last packet
C-----------------------------------------------------------------------
      IF(llt/=0) CALL i10lll(
     1       llt ,lll       ,jll       ,sll       ,xll     ,v       ,
     2       xx  ,yy        ,zz        ,iii       ,iadll   ,
     3       n_mul_mx ,a    ,x         ,itied     ,nint    ,nkmax   ,
     4       comntag )
C
C-----------------------------------------------
      RETURN
      END
!||====================================================================
!||    i10lll        ../engine/source/interfaces/int16/i10lagm.F
!||--- called by ------------------------------------------------------
!||    i10lagm       ../engine/source/interfaces/int16/i10lagm.F
!||--- calls      -----------------------------------------------------
!||    ancmsg        ../engine/source/output/message/message.F
!||    arret         ../engine/source/system/arret.F
!||    i10rst        ../engine/source/interfaces/int16/i10lagm.F
!||--- uses       -----------------------------------------------------
!||    message_mod   ../engine/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE i10lll(LLT  ,LLL       ,JLL  ,SLL  ,XLL  ,V     ,
     2                  XX   ,YY        ,ZZ   ,III  ,IADLL ,
     3                  N_MUL_MX,A      ,X    ,ITIED,NINT ,NKMAX ,
     4                  COMNTAG )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_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      "com08_c.inc"
      COMMON /lagglob/n_mult
      INTEGER N_MULT
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT,N_MUL_MX,ITIED,NINT ,NKMAX   
      INTEGER LLL(*),JLL(*),SLL(*),COMNTAG(*),
     .        iii(mvsiz,7),iadll(*)
C     REAL
      my_real
     .  xll(*),v(3,*),a(3,*)
      my_real
     .   xx(mvsiz,7),yy(mvsiz,7),zz(mvsiz,7),x(3,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, IK, NK, IAD, NN
      my_real
     .   vx,vy,vz,vn,aa
      my_real
     .   r(mvsiz),s(mvsiz),t(mvsiz),
     .   nx(mvsiz), ny(mvsiz), nz(mvsiz),
     .   ni(mvsiz,7) 
C-----------------------------------------------
C      calculation de r,s,t
C-----------------------------------------------
c
c      print *, "XX(1,1),XX(1,9)",XX(1,1),XX(1,9)
c
      CALL i10rst(llt   ,r     ,s     ,t     ,ni    ,
     2            nx    ,ny    ,nz    ,xx    ,yy    ,zz    )
C-----------------------------------------------
C      calculation of [L]
C-----------------------------------------------
      IF(itied==0)THEN
       DO i=1,llt
C-----------------------------------------------
C       Test if contact
C-----------------------------------------------
        IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
     .     r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
C
          nk = 7
          vx = zero
          vy = zero
          vz = zero
          DO ik=1,nk
            vx = vx - (v(1,iii(i,ik))+dt12*a(1,iii(i,ik)))*ni(i,ik)
            vy = vy - (v(2,iii(i,ik))+dt12*a(2,iii(i,ik)))*ni(i,ik)
            vz = vz - (v(3,iii(i,ik))+dt12*a(3,iii(i,ik)))*ni(i,ik)
          ENDDO
c
c
          vn = nx(i)*vx + ny(i)*vy + nz(i)*vz
C-----------------------------------------------
C         Test if incoming velocity in S
C-----------------------------------------------
          IF(s(i)*vn<=zero)THEN
c
c      print *, "velocity entrante",vn
c      print *, "s = ",S(I)
c
           aa = one/sqrt(nx(i)*nx(i)+ny(i)*ny(i)+nz(i)*nz(i))
           nx(i) = nx(i)*aa
           ny(i) = ny(i)*aa
           nz(i) = nz(i)*aa
#include "lockon.inc"
            n_mult=n_mult+1
            IF(n_mult>n_mul_mx)THEN
#include "lockoff.inc"
              CALL ancmsg(msgid=84,anmode=aninfo)
              CALL arret(2)
            ENDIF
            iadll(n_mult+1)=iadll(n_mult) + 21
            IF(iadll(n_mult+1)-1>nkmax)THEN
#include "lockoff.inc"
              CALL ancmsg(msgid=84,anmode=aninfo)
              CALL arret(2)
            ENDIF
            iad = iadll(n_mult) - 1
            DO ik=1,7
                lll(iad+ik)    = iii(i,ik)
                jll(iad+ik)    = 1
                sll(iad+ik)    = 0
                xll(iad+ik)    = nx(i)*ni(i,ik)
                lll(iad+ik+7)  = iii(i,ik)
                jll(iad+ik+7)  = 2
                sll(iad+ik+7)  = 0
                xll(iad+ik+7)  = ny(i)*ni(i,ik)
                lll(iad+ik+14) = iii(i,ik)
                jll(iad+ik+14) = 3
                sll(iad+ik+14) = 0
                xll(iad+ik+14) = nz(i)*ni(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
            sll(iad+7) = nint
            sll(iad+14) = nint
            sll(iad+21) = nint
#include "lockoff.inc"
          ENDIF
        ENDIF
       ENDDO
      ELSEIF(itied==1)THEN
C-----------------------------------------------
C      ITIED = 1
C-----------------------------------------------
       DO i=1,llt
C-----------------------------------------------
C       Test if contact
C-----------------------------------------------
        IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
     .     r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
C
          nk = 7
          vx = zero
          vy = zero
          vz = zero
          DO ik=1,nk
            vx = vx - (v(1,iii(i,ik))+dt12*a(1,iii(i,ik)))*ni(i,ik)
            vy = vy - (v(2,iii(i,ik))+dt12*a(2,iii(i,ik)))*ni(i,ik)
            vz = vz - (v(3,iii(i,ik))+dt12*a(3,iii(i,ik)))*ni(i,ik)
          ENDDO
c
c      print *, "vx,vy,vz s-m",vx,vy,vz
c      print *, "nx,ny,nz ", NX(I),NY(I),NZ(I)
c
          vn = nx(i)*vx + ny(i)*vy + nz(i)*vz
C-----------------------------------------------
C         Test if incoming velocity in S
C-----------------------------------------------
          IF(s(i)*vn<=zero)THEN
c          
c      print *, "velocity entrante",vn
c      print *, "s = ",S(I)
c
#include "lockon.inc"
            IF(n_mult+3>n_mul_mx)THEN
#include "lockoff.inc"
              CALL ancmsg(msgid=84,anmode=aninfo)
              CALL arret(2)
            ENDIF
            IF(iadll(n_mult+1)-1+7*3>nkmax)THEN
#include "lockoff.inc"
              CALL ancmsg(msgid=84,anmode=aninfo)
              CALL arret(2)
            ENDIF
C
            n_mult=n_mult+1
            iadll(n_mult+1)=iadll(n_mult) + 7
            iad = iadll(n_mult) - 1
            DO ik=1,7
                lll(iad+ik) = iii(i,ik)
                jll(iad+ik) = 1
                sll(iad+ik) = 0
                xll(iad+ik) = ni(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
            sll(iad+7) = nint
C
            n_mult=n_mult+1
            iadll(n_mult+1)=iadll(n_mult) + 7
            iad = iadll(n_mult) - 1
            DO ik=1,7
                lll(iad+ik) = iii(i,ik)
                jll(iad+ik) = 2
                sll(iad+ik) = 0
                xll(iad+ik) = ni(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
            sll(iad+7) = nint
C
            n_mult=n_mult+1
            iadll(n_mult+1)=iadll(n_mult) + 7
            iad = iadll(n_mult) - 1
            DO ik=1,7
                lll(iad+ik) = iii(i,ik)
                jll(iad+ik) = 3
                sll(iad+ik) = 0
                xll(iad+ik) = ni(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
            sll(iad+7) = nint
#include "lockoff.inc"
          ENDIF
        ENDIF
       ENDDO
      ELSE
C-----------------------------------------------
C      ITIED = 2
C-----------------------------------------------
       DO i=1,llt
C-----------------------------------------------
C       Test if contact
C-----------------------------------------------
        IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
     .     r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
C
           nk = 7
C-----------------------------------------------
c      print *, "s = ",S(I)
c
#include "lockon.inc"
            IF(n_mult+3>n_mul_mx)THEN
#include "lockoff.inc"
              CALL ancmsg(msgid=84,anmode=aninfo)
              CALL arret(2)
            ENDIF
            IF(iadll(n_mult+1)-1+7*3>nkmax)THEN
#include "lockoff.inc"
              CALL ancmsg(msgid=84,anmode=aninfo)
              CALL arret(2)
            ENDIF
            n_mult=n_mult+1
            iadll(n_mult+1)=iadll(n_mult) + 7
            iad = iadll(n_mult) - 1
            DO ik=1,7
                lll(iad+ik) = iii(i,ik)
                jll(iad+ik) = 1
                sll(iad+ik) = 0
                xll(iad+ik) = ni(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
            sll(iad+7) = nint
C
            n_mult=n_mult+1
            iadll(n_mult+1)=iadll(n_mult) + 7
            iad = iadll(n_mult) - 1
            DO ik=1,7
                lll(iad+ik) = iii(i,ik)
                jll(iad+ik) = 2
                sll(iad+ik) = 0
                xll(iad+ik) = ni(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
            sll(iad+7) = nint
C
            n_mult=n_mult+1
            iadll(n_mult+1)=iadll(n_mult) + 7
            iad = iadll(n_mult) - 1
            DO ik=1,7
                lll(iad+ik) = iii(i,ik)
                jll(iad+ik) = 3
                sll(iad+ik) = 0
                xll(iad+ik) = ni(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
            sll(iad+7) = nint
C
#include "lockoff.inc"
        ENDIF
       ENDDO
      ENDIF
c
c      print *, "r,s,t",r(1),s(1),t(1)
C
      RETURN
      END
C
!||====================================================================
!||    i10rst    ../engine/source/interfaces/int16/i10lagm.F
!||--- called by ------------------------------------------------------
!||    i10lll    ../engine/source/interfaces/int16/i10lagm.F
!||--- calls      -----------------------------------------------------
!||    i10deri   ../engine/source/interfaces/int16/i10lagm.F
!||    i10ni     ../engine/source/interfaces/int16/i10lagm.F
!||    i10rstn   ../engine/source/interfaces/int16/i10lagm.F
!||====================================================================
      SUBROUTINE i10rst(LLT,R  ,S     ,T     ,NI    ,
     2            NX  ,NY  ,NZ  ,XX    ,YY    ,ZZ    )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT
C     REAL
      my_real
     .   XX(MVSIZ,7),YY(MVSIZ,7),ZZ(MVSIZ,7)
      my_real
     .   R(MVSIZ),S(MVSIZ),T(MVSIZ),NI(MVSIZ,7) ,
     .   NX(MVSIZ),NY(MVSIZ),NZ(MVSIZ)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, ITER, NITERMAX, JTER, NJTERMAX, CONV
      my_real
     .   drdx(mvsiz),drdy(mvsiz),drdz(mvsiz),
     .   dsdx(mvsiz),dsdy(mvsiz),dsdz(mvsiz),     
     .   dtdx(mvsiz),dtdy(mvsiz),dtdz(mvsiz),
     .   dxdr(mvsiz),dydr(mvsiz),dzdr(mvsiz),
     .   dxdt(mvsiz),dydt(mvsiz),dzdt(mvsiz),
     .   rr(mvsiz),ss(mvsiz),tt(mvsiz)     
C-----------------------------------------------
C
C      r=s=t=0
C
C  +---> iter 
C  |                  
C  |     Ni(r,s,t) =  
C  |     dNi/dr    =  
C  |     ...         _ 
C  |                \
C  |     dx/dr    = /_ (xi * dNi/dr)
C  |     ...        
C  | 
C  |            [dx/dr dy/dr dz/dr]            
C  |      [J] = |dx/ds dy/ds dz/ds|                
C  |            [dx/dt dy/dt dz/dt]               
C  | 
C  | +--> jter                                 
C  | |                  _ 
C  | |                 \
C  | |      x(r,s,t) = /_ (xi * Ni(r,s,t))
C  | |      ...        
C  | |
C  | |      |r|   |r|      -1  |xs-x(r,s,t)|
C  | |      {s} = {s} + [J]    {ys-y(r,s,t)}
C  | |      |t|   |t|          |zs-z(r,s,t)|
C  | |                   
C  | |      Ni(r,s,t) =  
C  +-+---
C-----------------------------------------------
       nitermax = 3
       njtermax = 3
       conv = 0
C
       DO i=1,llt
         rr(i) = zero
         ss(i) = zero
         tt(i) = zero
       ENDDO
C-----------------------------------------------
C      calculation de r,s,t   et   Ni(r,s,t)
C-----------------------------------------------
       DO iter=1,nitermax
c
c      print *, "iter",iter
c
C-----------------------------------------------
C          calculation de Ni(r,s,t); [J]; [J]-1
C-----------------------------------------------
           CALL i10deri(llt,rr ,ss    ,tt    ,ni    ,
     2            drdx  ,drdy  ,drdz  ,dsdx  ,dsdy  ,dsdz  ,
     3            dtdx  ,dtdy  ,dtdz  ,dxdr  ,dydr  ,dzdr  ,
     4            dxdt  ,dydt  ,dzdt  ,xx    ,yy    ,zz    )
C
           DO jter=1,njtermax
c
c      print *, "jter",jter
c
C-----------------------------------------------
C            calculation de r,s,t new
C-----------------------------------------------
             CALL i10rstn(llt,rr,ss   ,tt    ,ni    ,conv  ,
     2            drdx  ,drdy  ,drdz  ,dsdx  ,dsdy  ,dsdz  ,
     3            dtdx  ,dtdy  ,dtdz  ,xx    ,yy    ,zz    ,
     4            r     ,s     ,t     )
c
c      print *, "r,s,t",r(1),s(1),t(1)
c      print *, "rr,ss,tt",rr(1),ss(1),tt(1)
c
C-----------------------------------------------
C            calculation de Ni(-1<r<1 , -1<s<1 , -1<t<1)
C-----------------------------------------------
             CALL i10ni(llt,rr ,ss ,tt ,ni  )
C  parity problem if convergence depends on mvsiz !!!!!!!
C             IF(CONV/=0)RETURN
C
           ENDDO
      ENDDO
C
       DO i=1,llt
         nx(i) = dydt(i)*dzdr(i) - dzdt(i)*dydr(i)
         ny(i) = dzdt(i)*dxdr(i) - dxdt(i)*dzdr(i)
         nz(i) = dxdt(i)*dydr(i) - dydt(i)*dxdr(i)
       ENDDO
C
      RETURN
      END
!||====================================================================
!||    i10ni    ../engine/source/interfaces/int16/i10lagm.F
!||--- called by ------------------------------------------------------
!||    i10rst   ../engine/source/interfaces/int16/i10lagm.F
!||====================================================================
      SUBROUTINE i10ni(LLT,RR ,SS ,TT ,NI    )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT
      my_real
     .   RR(MVSIZ),SS(MVSIZ),TT(MVSIZ),NI(MVSIZ,7)  
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  u_m_r,u_p_r,u_m_s,u_p_s,u_m_t,u_p_t,
     .  ums_umt,ums_upt,ups_umt,ups_upt,
     .  umr_ums,umr_ups,upr_ums,upr_ups,
     .  umt_umr,umt_upr,upt_umr,upt_upr,
     .  a,r05,s05,t05
C-----------------------------------------------------------------------
C     calculation of Ni
C-----------------------------------------------------------------------
      DO i=1,llt
C
        r05 = half*rr(i)
        s05 = half*ss(i)
        t05 = half*tt(i)
C
        u_m_r = half - r05
        u_p_r = half + r05
C
        u_m_s = half - s05
        u_p_s = half + s05
C
        u_m_t = half - t05
        u_p_t = half + t05
C
        ums_umt = u_m_s * u_m_t
        ums_upt = u_m_s * u_p_t
        ups_umt = u_p_s * u_m_t
        ups_upt = u_p_s * u_p_t
C
        umr_ums = u_m_r * u_m_s
        umr_ups = u_m_r * u_p_s
        upr_ums = u_p_r * u_m_s
        upr_ups = u_p_r * u_p_s
C
        umt_umr = u_m_t * u_m_r
        umt_upr = u_m_t * u_p_r
        upt_umr = u_p_t * u_m_r
        upt_upr = u_p_t * u_p_r
C
        a = -rr(i)-tt(i)-one
        ni(i,1) = u_m_r * ums_umt * a
        ni(i,2) = u_m_r * ums_upt * a
        ni(i,3) = u_p_r * ums_upt * a
        ni(i,4) = u_p_r * ums_umt * a
        ni(i,5) = u_m_r * ups_umt * a
        ni(i,6) = u_m_r * ups_upt * a
C------------------------------------
        ni(i,7) = -1.
C------------------------------------
      ENDDO
C-----------------------------------------------
      RETURN
      END
!||====================================================================
!||    i10rstn   ../engine/source/interfaces/int16/i10lagm.F
!||--- called by ------------------------------------------------------
!||    i10rst    ../engine/source/interfaces/int16/i10lagm.F
!||====================================================================
      SUBROUTINE i10rstn(LLT,RR,SS    ,TT    ,NI    ,CONV  ,
     2            DRDX  ,DRDY  ,DRDZ  ,DSDX  ,DSDY  ,DSDZ  ,
     3            DTDX  ,DTDY  ,DTDZ  ,XX    ,YY    ,ZZ    ,
     4            R     ,S     ,T     )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
c#include      "implicit_f.inc"
      implicit none
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
#include      "constant.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT,CONV
      my_real
     .   r(mvsiz),s(mvsiz),t(mvsiz),ni(mvsiz,7) ,
     .   rr(mvsiz),ss(mvsiz),tt(mvsiz),
     .   xx(mvsiz,7) ,yy(mvsiz,7) ,zz(mvsiz,7) ,
     .   drdx(mvsiz),drdy(mvsiz),drdz(mvsiz),
     .   dsdx(mvsiz),dsdy(mvsiz),dsdz(mvsiz),     
     .   dtdx(mvsiz),dtdy(mvsiz),dtdz(mvsiz)     
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .   dx ,dy,dz,dr ,ds,dt,err
C
      err=zero
C-----------------------------------------------
      DO i=1,llt
C
        dx = xx(i,7)
     +     - ni(i, 1)*xx(i, 1) - ni(i, 2)*xx(i, 2) - ni(i, 3)*xx(i, 3)
     +     - ni(i, 4)*xx(i, 4) - ni(i, 5)*xx(i, 5) - ni(i, 6)*xx(i, 6)
        dy = yy(i,7)
     +     - ni(i, 1)*yy(i, 1) - ni(i, 2)*yy(i, 2) - ni(i, 3)*yy(i, 3)
     +     - ni(i, 4)*yy(i, 4) - ni(i, 5)*yy(i, 5) - ni(i, 6)*yy(i, 6)
        dz = zz(i,7)
     +     - ni(i, 1)*zz(i, 1) - ni(i, 2)*zz(i, 2) - ni(i, 3)*zz(i, 3)
     +     - ni(i, 4)*zz(i, 4) - ni(i, 5)*zz(i, 5) - ni(i, 6)*zz(i, 6)
C
        dr = drdx(i)*dx + drdy(i)*dy + drdz(i)*dz
        ds = dsdx(i)*dx + dsdy(i)*dy + dsdz(i)*dz 
        dt = dtdx(i)*dx + dtdy(i)*dy + dtdz(i)*dz
C
        rr(i) = rr(i) + dr
        ss(i) = ss(i) + ds
        tt(i) = tt(i) + dt
C
        r(i) = rr(i)
        s(i) = ss(i)
        t(i) = tt(i)
C
        IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
     .     r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
          err = max(err,abs(dr),abs(ds),abs(dt))
        ELSE
          rr(i) = max(min(rr(i),one),-one)
          ss(i) = max(min(ss(i),one),-one)
          tt(i) = max(min(tt(i),one),-one)
        ENDIF
c
C
      ENDDO
C
      IF(err<em4) conv = 1
C-----------------------------------------------
      RETURN
      END
!||====================================================================
!||    i10deri   ../engine/source/interfaces/int16/i10lagm.F
!||--- called by ------------------------------------------------------
!||    i10rst    ../engine/source/interfaces/int16/i10lagm.F
!||====================================================================
      SUBROUTINE i10deri(LLT,RR,SS ,TT    ,NI    ,
     2   DRDX  ,DRDY  ,DRDZ  ,DSDX  ,DSDY  ,DSDZ  ,
     3   DTDX  ,DTDY  ,DTDZ  ,DXDR  ,DYDR  ,DZDR  ,
     4   DXDT  ,DYDT  ,DZDT  ,XX    ,YY    ,ZZ    )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT
      my_real
     .   DXDR(MVSIZ), DYDR(MVSIZ), DZDR(MVSIZ),
     .   DXDT(MVSIZ), DYDT(MVSIZ), DZDT(MVSIZ),
     .   DRDX(MVSIZ), DSDX(MVSIZ), DTDX(MVSIZ),
     .   DRDY(MVSIZ), DSDY(MVSIZ), DTDY(MVSIZ),
     .   DRDZ(MVSIZ), DSDZ(MVSIZ), DTDZ(MVSIZ),
     .   XX(MVSIZ,7) ,YY(MVSIZ,7),ZZ(MVSIZ,7),
     .   NI(MVSIZ,7) ,RR(MVSIZ) ,SS(MVSIZ) ,TT(MVSIZ)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .   dxds(mvsiz), dyds(mvsiz), dzds(mvsiz),
     .   dnidr(10),dnids(10),dnidt(10),
     . d, det(mvsiz)
      my_real
     .  u_m_r,u_p_r,u_m_s,u_p_s,u_m_t,u_p_t,
     .  ums_umt,ums_upt,ups_umt,ups_upt,
     .  umr_ums,umr_ups,upr_ums,upr_ups,
     .  umt_umr,umt_upr,upt_umr,upt_upr,
     .  a,r05,s05,t05
C-----------------------------------------------
C/*
C
C
C*/
C-----------------------------------------------
C
C-----------------------------------------------
C                _ 
C               \
C    x(r,s,t) = /_ (xi * Ni(r,s,t))
C                _ 
C               \
C    y(r,s,t) = /_ (yi * Ni(r,s,t))
C                _ 
C               \
C    z(r,s,t) = /_ zi * Ni(r,s,t))
C        
C                _ 
C               \
C    dx/dr    = /_ (xi * dNi/dr)
C    ...        
C
C          [dx/dr dy/dr dz/dr]            
C    [J] = |dx/ds dy/ds dz/ds|                
C          [dx/dt dy/dt dz/dt]               
C                                     
C    |r|   |r|      -1  |xs-x|
C    {s} = {s} + [J]    {ys-y}
C    |t|   |t|          |zs-z|
C
C-----------------------------------------------------------------------
C     Ni; dNi/dr; dNi/ds; dNi/dt
C-----------------------------------------------------------------------
      DO i=1,llt
        r05 = half*rr(i)
        s05 = half*ss(i)
        t05 = half*tt(i)
C
        u_m_r = half - r05
        u_p_r = half + r05
C
        u_m_s = half - s05
        u_p_s = half + s05
C
        u_m_t = half - t05
        u_p_t = half + t05
C
        ums_umt = u_m_s * u_m_t
        ums_upt = u_m_s * u_p_t
        ups_umt = u_p_s * u_m_t
        ups_upt = u_p_s * u_p_t
C
        umr_ums = u_m_r * u_m_s
        umr_ups = u_m_r * u_p_s
        upr_ums = u_p_r * u_m_s
        upr_ups = u_p_r * u_p_s
C
        umt_umr = u_m_t * u_m_r
        umt_upr = u_m_t * u_p_r
        upt_umr = u_p_t * u_m_r
        upt_upr = u_p_t * u_p_r
C
        a = -rr(i)-tt(i)-one
        ni(i,1) = u_m_r * ums_umt * a
        ni(i,2) = u_m_r * ums_upt * a
        ni(i,3) = u_p_r * ums_upt * a
        ni(i,4) = u_p_r * ums_umt * a
        ni(i,5) = u_m_r * ups_umt * a
        ni(i,6) = u_m_r * ups_upt * a
C
        a = -t05 - rr(i)
        dnidr(1) = -ums_umt * a 
        dnidr(5) = -ups_umt * a
        dnidr(2) = -ums_upt * a
        dnidr(6) = -ups_upt * a
        dnidr(3) =  ums_upt * a
        dnidr(4) =  ums_umt * a
C
        dnids(1) = -umt_umr * a
        dnids(5) =  umt_umr * a
        dnids(2) = -upt_umr * a
        dnids(6) =  upt_umr * a
        dnids(3) = -upt_upr * a
        dnids(4) = -umt_upr * a
C
        dnidt(1) = -umr_ums * a
        dnidt(5) = -umr_ups * a
        dnidt(2) =  umr_ums * a
        dnidt(6) =  umr_ups * a
        dnidt(3) =  upr_ums * a
        dnidt(4) = -upr_ums * a
C------------------------------------
        ni(i,7) = -1.
C-----------------------------------------------------------------------
C     dx/dr; dx/ds; dx/dt
C-----------------------------------------------------------------------
        dxdr(i) = dnidr(1)*xx(i,1) + dnidr(2)*xx(i,2) + dnidr(3)*xx(i,3)
     +          + dnidr(4)*xx(i,4) + dnidr(5)*xx(i,5) + dnidr(6)*xx(i,6)
C  
        dxds(i) = dnids(1)*xx(i,1) + dnids(2)*xx(i,2) + dnids(3)*xx(i,3)
     +          + dnids(4)*xx(i,4) + dnids(5)*xx(i,5) + dnids(6)*xx(i,6)
C 
        dxdt(i) = dnidt(1)*xx(i,1) + dnidt(2)*xx(i,2) + dnidt(3)*xx(i,3)
     +          + dnidt(4)*xx(i,4) + dnidt(5)*xx(i,5) + dnidt(6)*xx(i,6)
C-----------------------------------------------------------------------
C     dy/dr; dy/ds; dy/dt
C-----------------------------------------------------------------------
        dydr(i) = dnidr(1)*yy(i,1) + dnidr(2)*yy(i,2) + dnidr(3)*yy(i,3)
     +          + dnidr(4)*yy(i,4) + dnidr(5)*yy(i,5) + dnidr(6)*yy(i,6)
C  
        dyds(i) = dnids(1)*yy(i,1) + dnids(2)*yy(i,2) + dnids(3)*yy(i,3)
     +          + dnids(4)*yy(i,4) + dnids(5)*yy(i,5) + dnids(6)*yy(i,6)
C 
        dydt(i) = dnidt(1)*yy(i,1) + dnidt(2)*yy(i,2) + dnidt(3)*yy(i,3)
     +          + dnidt(4)*yy(i,4) + dnidt(5)*yy(i,5) + dnidt(6)*yy(i,6)
C-----------------------------------------------------------------------
C     dz/dr; dz/ds; dz/dt
C-----------------------------------------------------------------------
        dzdr(i) = dnidr(1)*zz(i,1) + dnidr(2)*zz(i,2) + dnidr(3)*zz(i,3)
     +          + dnidr(4)*zz(i,4) + dnidr(5)*zz(i,5) + dnidr(6)*zz(i,6)
C  
        dzds(i) = dnids(1)*zz(i,1) + dnids(2)*zz(i,2) + dnids(3)*zz(i,3)
     +          + dnids(4)*zz(i,4) + dnids(5)*zz(i,5) + dnids(6)*zz(i,6)
C 
        dzdt(i) = dnidt(1)*zz(i,1) + dnidt(2)*zz(i,2) + dnidt(3)*zz(i,3)
     +          + dnidt(4)*zz(i,4) + dnidt(5)*zz(i,5) + dnidt(6)*zz(i,6)
C-----------------------------------------------------------------------
C          -1
C       [J]          inversion of the jacobian
C-----------------------------------------------------------------------
        drdx(i)=dyds(i)*dzdt(i)-dzds(i)*dydt(i)
        drdy(i)=dzds(i)*dxdt(i)-dxds(i)*dzdt(i)
        drdz(i)=dxds(i)*dydt(i)-dyds(i)*dxdt(i)
C
        dsdz(i)=dxdt(i)*dydr(i)-dydt(i)*dxdr(i)
        dsdy(i)=dzdt(i)*dxdr(i)-dxdt(i)*dzdr(i)
        dsdx(i)=dydt(i)*dzdr(i)-dzdt(i)*dydr(i)
C
        dtdx(i)=dydr(i)*dzds(i)-dzdr(i)*dyds(i)
        dtdy(i)=dzdr(i)*dxds(i)-dxdr(i)*dzds(i)
        dtdz(i)=dxdr(i)*dyds(i)-dydr(i)*dxds(i)
C
        det(i) = dxdr(i) * drdx(i)
     .         + dydr(i) * drdy(i)
     .         + dzdr(i) * drdz(i)
C
c
c
      ENDDO
C
      DO i=1,llt
C-----------------------------------------------------------------------
C          -1            
C       [J] Inversion of the Jacobian Suite
C-----------------------------------------------------------------------
        d = one/det(i)
        drdx(i)=d*drdx(i)
        dsdx(i)=d*dsdx(i)
        dtdx(i)=d*dtdx(i)
C
        drdy(i)=d*drdy(i)
        dsdy(i)=d*dsdy(i)
        dtdy(i)=d*dtdy(i)
C
        drdz(i)=d*drdz(i)
        dsdz(i)=d*dsdz(i)
        dtdz(i)=d*dtdz(i)
C
c
c      print *, "DRDX(I),DRDY(I),DRDZ(I)",DRDX(I),DRDY(I),DRDZ(I)
c      print *, "DSDX(I),DSDY(I),DSDZ(I)",DSDX(I),DSDY(I),DSDZ(I)
c      print *, "DTDX(I),DTDY(I),DTDZ(I)",DTDX(I),DTDY(I),DTDZ(I)
c
      ENDDO
C-----------------------------------------------------------------------
      RETURN
      END