s10forc3.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/.
!||====================================================================
!||    s10forc3               ../engine/source/elements/solid/solide10/s10forc3.f
!||--- called by ------------------------------------------------------
!||    forint                 ../engine/source/elements/forint.F
!||--- calls      -----------------------------------------------------
!||    boltst                 ../engine/source/elements/solid/solide/boltst.F
!||    mmain                  ../engine/source/materials/mat_share/mmain.F90
!||    nsvis_sm12             ../engine/source/elements/solid/solide10/nsvis_sm12.F
!||    s10_icp                ../engine/source/elements/solid/solide4_sfem/s10_icp.f
!||    s10bilan               ../engine/source/elements/solid/solide10/s10bilan.F
!||    s10coor3               ../engine/source/elements/solid/solide10/s10coor3.F
!||    s10cumu3               ../engine/source/elements/solid/solide10/s10cumu3.F
!||    s10cumu3p              ../engine/source/elements/solid/solide10/s10cumu3p.F
!||    s10defo3               ../engine/source/elements/solid/solide10/s10defo3.F
!||    s10defot3              ../engine/source/elements/solid/solide10/s10defot3.F
!||    s10deri3               ../engine/source/elements/solid/solide10/s10deri3.F
!||    s10derit3              ../engine/source/elements/solid/solide10/s10derit3.F
!||    s10divde12             ../engine/source/elements/solid/solide10/s10divde12.F
!||    s10dvm12               ../engine/source/elements/solid/solide10/s10dvm12.F
!||    s10fint3               ../engine/source/elements/solid/solide10/s10fint3.F
!||    s10for_distor          ../engine/source/elements/solid/solide10/s10for_distor.f
!||    s10get_x0              ../engine/source/elements/solid/solide10/s10get_x0.F90
!||    s10get_x3              ../engine/source/elements/solid/solide10/s10get_x3.F
!||    s10len3                ../engine/source/elements/solid/solide10/s10len3.F
!||    s10malla3              ../engine/source/elements/solid/solide10/s10malla3.F
!||    s10mallb3              ../engine/source/elements/solid/solide10/s10mallb3.F
!||    s10mallgeo3            ../engine/source/elements/solid/solide10/s10mallgeo3.f
!||    s10nx3                 ../engine/source/elements/solid/solide10/s10nx3.f
!||    s10nxt4                ../engine/source/elements/solid/solide10/s10nxt4.F
!||    s10pijto3              ../engine/source/elements/solid/solide10/s10pijto3.F
!||    s10rcoor12             ../engine/source/elements/solid/solide10/s10rcoor12.F
!||    s10sav12               ../engine/source/elements/solid/solide10/s10sav12.F
!||    s10sav3                ../engine/source/elements/solid/solide10/s10sav3.F
!||    s10sigp3               ../engine/source/elements/solid/solide4_sfem/s10sigp3.F
!||    s10therm               ../engine/source/elements/solid/solide10/s10therm.F
!||    s10upd11t12            ../engine/source/elements/solid/solide10/s10upd11t12.F
!||    s4therm_itet1          ../engine/source/elements/solid/solide4/s4therm-itet1.f
!||    s8e_sigp               ../engine/source/elements/solid/solide8e/s8e_sig.F
!||    sdistor_ini            ../engine/source/elements/solid/solide/sdistror_ini.f90
!||    sgcoor10               ../engine/source/elements/solid/solide10/sgcoor10.F
!||    sgeodel3               ../engine/source/elements/solid/solide/sgeodel3.F
!||    smallb3                ../engine/source/elements/solid/solide/smallb3.F
!||    smallgeo3              ../engine/source/elements/solid/solide/smallgeo3.F
!||    sordef12               ../engine/source/elements/solid/solidez/sordef12.F
!||    sordeft12              ../engine/source/elements/solid/solidez/sordeft12.F
!||    sorthdir3              ../engine/source/elements/solid/solide/sorthdir3.F
!||    sreploc3               ../engine/source/elements/solid/solide/sreploc3.F
!||    srho3                  ../engine/source/elements/solid/solide/srho3.F
!||    srota3                 ../engine/source/elements/solid/solide/srota3.F
!||    sroto12_sig            ../engine/source/elements/solid/solidez/sroto12_sig.F
!||    sstra3                 ../engine/source/elements/solid/solide/sstra3.F
!||    startime               ../engine/source/system/timer_mod.F90
!||    stoptime               ../engine/source/system/timer_mod.F90
!||    sxfillopt              ../engine/source/elements/solid/solide10/sxfillopt.F
!||--- uses       -----------------------------------------------------
!||    ale_connectivity_mod   ../common_source/modules/ale/ale_connectivity_mod.F
!||    dt_mod                 ../engine/source/modules/dt_mod.f
!||    elbufdef_mod           ../common_source/modules/mat_elem/elbufdef_mod.f90
!||    glob_therm_mod         ../common_source/modules/mat_elem/glob_therm_mod.F90
!||    mat_elem_mod           ../common_source/modules/mat_elem/mat_elem_mod.F90
!||    mmain_mod              ../engine/source/materials/mat_share/mmain.F90
!||    nlocal_reg_mod         ../common_source/modules/nlocal_reg_mod.F
!||    output_mod             ../common_source/modules/output/output_mod.f90
!||    s10get_x0_mod          ../engine/source/elements/solid/solide10/s10get_x0.F90
!||    sdistor_ini_mod        ../engine/source/elements/solid/solide/sdistror_ini.f90
!||    sensor_mod             ../common_source/modules/sensor_mod.F90
!||    table_mod              ../engine/share/modules/table_mod.F
!||    timer_mod              ../engine/source/system/timer_mod.F90
!||====================================================================
      SUBROUTINE s10forc3(TIMERS, OUTPUT, 
     1   ELBUF_TAB,       NG,         PM,         GEO,
     2   IXS,             X,          A,          V,
     3   MS,              W,          FLUX,       FLU1,
     4   VEUL,            FV,         ALE_CONNECT,IPARG,
     5   TF,              NPF,        BUFMAT,     PARTSAV,
     6   NLOC_DMG,        DT2T,       NELTST,     ITYPTST,
     7   STIFN,           FSKY,       IADS,       OFFSET,
     8   EANI,            IPARTS,     IXS10,      IADS10,
     9   NEL,             FX,         FY,         FZ,
     A   AR,              VR,         DR,         IPM,
     B   ISTRAIN,         ISOLNOD,    ITASK,      TEMP,
     C   FTHE,            FTHESKY,    IEXPAN,     STIFR,
     D   D,               GRESAV,     GRTH,       IGRTH,
     E   TABLE,           MSSA,       DMELS,      IGEO,
     F   XDP,             VOLN,       CONDN,      CONDNSKY,
     G   SFEM_NODVAR,     ITAGDN,     SENSORS,    IOUTPRT,
     H   MAT_ELEM,        H3D_STRAIN, DT,         IDEL7NOK,
     I   SNPC,            STF,        SBUFMAT,    SVIS,
     J   NSVOIS,          IDTMINS,    IDEL7NG,    MAXFUNC,
     V   USERL_AVAIL,     GLOB_THERM, impl_s,     idyna  ,
     W   S_SFEM_NODVAR)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE timer_mod
      USE output_mod, only : output_
      USE mmain_mod
      USE table_mod
      USE mat_elem_mod            
      USE nlocal_reg_mod
      USE ale_connectivity_mod
      USE sensor_mod
      USE dt_mod
      USE elbufdef_mod
      USE sdistor_ini_mod, ONLY : sdistor_ini
      use glob_therm_mod
      USE s10get_x0_mod, ONLY : s10get_x0
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      "com01_c.inc"
#include      "com04_c.inc"
#include      "com08_c.inc"
#include      "scr03_c.inc"
#include      "vect01_c.inc"
#include      "parit_c.inc"
#include      "param_c.inc"
#include      "timeri_c.inc"
#include      "scr18_c.inc"
#include      "scr05_c.inc"
C-----------------------------------------------
C   l o c a l   P a r a m e t e r s
C-----------------------------------------------
      INTEGER NPE
      PARAMETER (NPE=10)
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      type(timer_), INTENT(INOUT) :: timers
      TYPE(output_), INTENT(INOUT) :: OUTPUT
      INTEGER, INTENT(IN) :: S_SFEM_NODVAR
      INTEGER, INTENT(INOUT) :: IDEL7NOK
      INTEGER, INTENT(IN) :: SNPC
      INTEGER, INTENT(IN) :: STF
      INTEGER, INTENT(IN) :: SBUFMAT
      INTEGER, INTENT(IN) :: NSVOIS
      INTEGER, INTENT(IN) :: IDTMINS
      INTEGER ,INTENT(IN) :: IDEL7NG
      INTEGER ,INTENT(IN) :: MAXFUNC
      INTEGER, INTENT(IN) :: IMPL_S
      INTEGER, INTENT(IN) :: IDYNA
      INTEGER, INTENT(IN) :: USERL_AVAIL
      INTEGER IXS(NIXS,*),IPARG(NPARG,NGROUP),NPF(*),IADS(8,*),
     .        IPARTS(*),IXS10(6,*),IADS10(6,*),IPM(*),ITASK,GRTH(*),
     .        igrth(*),igeo(npropgi,*),ioutprt,h3d_strain
      INTEGER NELTST,ITYPTST,OFFSET,NEL,NG,ISTRAIN,ISOLNOD,IEXPAN,ITAGDN(*)      
      DOUBLE PRECISION  XDP(3,*)
      
      my_real DT2T
      my_real pm(npropm,*), geo(npropg,*), x(*), a(*), v(3,*), ms(*), w(*),
     .        flux(6,*),flu1(*), veul(*), fv(*), tf(*),
     .        bufmat(*),partsav(*),stifn(*), fsky(*),eani(*),
     .        ar(*),vr(*) ,dr(*) ,stifr(*),d(*), mssa(*)  ,dmels(*)
      my_real fx(mvsiz,10),fy(mvsiz,10),fz(mvsiz,10),
     .        temp(*), fthe(*), fthesky(*),gresav(*),voln(mvsiz),condn(*),
     .        condnsky(*),sfem_nodvar(s_sfem_nodvar)
      my_real, DIMENSION(MVSIZ,6), INTENT(INOUT) :: svis
      TYPE(ttable) TABLE(*)
      TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
      TYPE (NLOCAL_STR_)  , TARGET :: NLOC_DMG 
      TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
      TYPE (MAT_ELEM_) ,INTENT(INOUT) :: MAT_ELEM
      TYPE (SENSORS_)  ,INTENT(INOUT) :: SENSORS
      TYPE(dt_)        ,INTENT(INOUT) :: DT
      type (glob_therm_) ,intent(inout)   :: glob_therm
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,IP,LCO,NF1,NF2,IFLAG,IOFFS,IPTR,IPTS,IPTT,ILAY,IMAT
      INTEGER IBID,IBIDON(1),NNEGA,INDEX(MVSIZ),ITET,iel,ISM12_11
C-----
      INTEGER MXT(MVSIZ),NGL(MVSIZ),NGEO(MVSIZ)
      my_real
     . VD2(MVSIZ) , DVOL(MVSIZ),DELTAX(MVSIZ),
     . vis(mvsiz) , qvis(mvsiz), cxx(mvsiz) ,deltax2(mvsiz),
     . s1(mvsiz)  , s2(mvsiz)  , s3(mvsiz)  ,
     . s4(mvsiz)  , s5(mvsiz)  , s6(mvsiz)  ,
     . dxx(mvsiz) , dyy(mvsiz) , dzz(mvsiz) ,
     . d4(mvsiz)  , d5(mvsiz)  , d6(mvsiz)  , 
     . rx(mvsiz) , ry(mvsiz) , rz(mvsiz) ,
     . sx(mvsiz) , sy(mvsiz) , sz(mvsiz) ,
     . vdx(mvsiz), vdy(mvsiz), vdz(mvsiz),ssp_eq(mvsiz),aire(mvsiz),
     . conde(mvsiz),condeg(mvsiz), volg(mvsiz), jacgm(mvsiz)
C-----
C Arrays for material laws
      my_real
     .   sti(mvsiz),
     .   wxx(mvsiz) , wyy(mvsiz) , wzz(mvsiz),
     .   wxxg(mvsiz) , wyyg(mvsiz) , wzzg(mvsiz)
C Arrays for material laws (sph only)
      my_real
     .   muvoid(mvsiz)
C-----
C     void MMAIN
      my_real
     . sigy(mvsiz),et(mvsiz),gama(mvsiz,6),
     . r1_free(mvsiz),r3_free(mvsiz),r4_free(mvsiz)
C For solid elements
      INTEGER NC(MVSIZ,10),ICP,MX,IPLAW1
 
      DOUBLE PRECISION
     .  XX(MVSIZ,10), YY(MVSIZ,10), ZZ(MVSIZ,10),
     .  WXX0(MVSIZ) , WYY0(MVSIZ) , WZZ0(MVSIZ),
     .  XX0(MVSIZ,10), YY0(MVSIZ,10), ZZ0(MVSIZ,10),VOLDP(MVSIZ,5)
 
      my_real
     .  tx(mvsiz),ty(mvsiz),tz(mvsiz),off(mvsiz),volp(mvsiz,5),
     .  rhoo(mvsiz),offs(mvsiz),them(mvsiz,10),tempel(mvsiz),
     .  vx(mvsiz,10),vy(mvsiz,10),vz(mvsiz,10),
     .  px(mvsiz,10,5),py(mvsiz,10,5),pz(mvsiz,10,5),
     .  nx(mvsiz,10,5),vdxx(mvsiz,10),vdyy(mvsiz,10),vdzz(mvsiz,10),
     .  dxy(mvsiz),dyx(mvsiz),dyz(mvsiz),dzy(mvsiz),
     .  dzx(mvsiz),dxz(mvsiz),
     .  stig(mvsiz),  wip(5,5), alph(5,5), beta(5,5),bid(mvsiz),
     .  die(mvsiz), mbid(1),offg0(mvsiz),amu(mvsiz),sum,rho0_1,cns2
      my_real nxt4(mvsiz,4,4)
      my_real, DIMENSION(MVSIZ) ::
     .   e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z
C
      my_real
     .    vx0(mvsiz,10),vy0(mvsiz,10),vz0(mvsiz,10),
     .    mfxx(mvsiz,5),mfxy(mvsiz,5),mfyx(mvsiz,5),
     .    mfyy(mvsiz,5),mfyz(mvsiz,5),mfzy(mvsiz,5),
     .    mfzz(mvsiz,5),mfzx(mvsiz,5),mfxz(mvsiz,5),divde(mvsiz),
     .    nu(mvsiz),fac(mvsiz),facp(mvsiz),e0(mvsiz),c1,dvm(mvsiz),
     .    visp(mvsiz),facdb,rbid(mvsiz),sigp(nel,6),
     .    fld(mvsiz),sti_c(mvsiz),ll(mvsiz),offg(mvsiz),fqmax
c
      my_real varnl(nel),deltax4(mvsiz)
c     Flag Bolt Preloading
      INTEGER IBOLTP,NBPRELD,II(6),ISCTL,ISTAB(MVSIZ)
      INTEGER SZ_IX
      my_real, 
     .  DIMENSION(:), POINTER :: BPRELD
      my_real, dimension(mvsiz) :: fheat
C-----
      TYPE(g_bufel_) ,POINTER :: GBUF
      TYPE(l_bufel_) ,POINTER :: LBUF     
C--------------------------------------------------------
      DATA WIP / 1.  ,0.  ,0.  ,0.  ,0.  ,
     2           0.  ,0.  ,0.  ,0.  ,0.  ,
     3           0.  ,0.  ,0.  ,0.  ,0.  ,
     4           0.25,0.25,0.25,0.25,0.  ,
     5           0.45,0.45,0.45,0.45,-0.8/
C-----------------------------------------------
C   S o u r c e  L i n e s
C=======================================================================
      gbuf  => elbuf_tab(ng)%GBUF
      iboltp = iparg(72,ng)
      nbpreld = gbuf%G_BPRELD
      bpreld =>gbuf%BPRELD(1:nbpreld*nel)
      ism12_11 = elbuf_tab(ng)%BUFLY(1)%L_SIGL
C-----
      sz_ix=numelq+numels+nsvois
      nf1=nft+1
      nf2=nf1-numels8
      ibid = 0
      ibidon  = 0
      ioffs=0
       ipts = 1
       iptt = 1
       ilay = 1
      IF(isrot == 1) THEN
        iisrot=1
        nf2=1
      END IF
      icp = iparg(10,ng)
      DO i=lft,llt
        offs(i)=ep20
      END DO
C      
      DO ip=1,3
        DO j=1,5
          alph(j,ip)=zero
          beta(j,ip)=zero
        END DO
      END DO
c
      alph(1,4)=zep5854102
      alph(2,4)=zep5854102
      alph(3,4)=zep5854102
      alph(4,4)=zep5854102
      alph(5,4)=zero
      alph(1,5)=half
      alph(2,5)=half
      alph(3,5)=half
      alph(4,5)=half
      alph(5,5)=fourth
      beta(1,4)=zep1381966
      beta(2,4)=zep1381966
      beta(3,4)=zep1381966
      beta(4,4)=zep1381966
      beta(5,4)=zero
      beta(1,5)=one_over_6
      beta(2,5)=one_over_6
      beta(3,5)=one_over_6
      beta(4,5)=one_over_6
      beta(5,5)=fourth
C   
      tempel(:) = zero
      fheat(:)  = zero
      IF (jthe < 0) them(lft:llt,1:10) = zero
C      
      IF (icp==1) THEN
        mx = ixs(1,nf1)
        nu(lft:llt)=min(half,pm(21,mx))
        facp(lft:llt)=one
      ELSEIF (icp==2) THEN
        mx = ixs(1,nf1)
        nu(lft:llt)=min(half,pm(21,mx))
        c1 =pm(32,mx)
        e0(lft:llt) =three*(one-two*nu(lft:llt))*c1
        sigp=zero   
        IF (gbuf%G_PLA>0) THEN
          CALL s8e_sigp(elbuf_tab(ng),sigp,         nel)
        END IF 
        CALL s10sigp3(sigp,e0 ,gbuf%PLA,facp ,gbuf%G_PLA,nel )
      END IF
C-----------
      CALL s10coor3(
     1   x,           ixs(1,nf1),  ixs10(1,nf2),v,
     2   w,           xx,          yy,          zz,
     3   vx,          vy,          vz,          vdxx,
     4   vdyy,        vdzz,        vdx,         vdy,
     5   vdz,         vd2,         vis,         gbuf%OFF,
     6   off,         gbuf%SMSTR,  nc,          ngl,
     7   mxt,         ngeo,        fx,          fy,
     8   fz,          stig,        gbuf%SIG,    gbuf%EINT,
     9   gbuf%RHO,    gbuf%QVIS,   gbuf%PLA,    gbuf%EPSD,
     a   vr,          dr,          d,           wxxg,
     b   wyyg,        wzzg,        gbuf%G_PLA,  xdp,
     c   nel,         condeg,      gbuf%G_EPSD, jale,
     d   ismstr,      jeul,        jlag,        israt,
     e   isrot)
C     
      iplaw1 = 0
      cns2 = zero
      IF (ism12_11>0 .AND.idtmin(1)==3) THEN
         mx = ixs(1,nf1)
         rho0_1 =pm( 1,mx)
         IF (pm(21,mx)>0.49) iplaw1=1
         IF (iplaw1==1) THEN
           facdb = one- zep02
           facdb = min(facdb,two*pm(21,mx))
           facp(lft:llt)=facdb
           visp(lft:llt)=two
           cns2 = zep02
           IF (igeo(35,ngeo(1))>0) cns2=cns2-abs(geo(17,ngeo(1)))
         END IF 
      ELSEIF (ismstr==10.AND.mtn==1) THEN
         mx = ixs(1,nf1)
         rho0_1 =pm( 1,mx)
         IF (pm(21,mx)>0.49) THEN
           visp(lft:llt)=two
           cns2 = zep02
           IF (igeo(35,ngeo(1))>0) cns2=cns2-abs(geo(17,ngeo(1)))
         END IF 
      END IF 
      isctl = igeo(97,ngeo(1))
      IF (isrot == 1) isctl = 0 ! uncompatible with itetra4=1
C     
      CALL s10nx3(
     1   nx,      nel,     npt)
      IF(jthe < 0 .AND. isolnod == 4) CALL s10nxt4(nxt4,nel)
C-----------
C GATHER NODAL VARIABLES FOR TOTAL STRAIN CASE.
C-----------
      IF (ismstr >= 10.AND.ismstr <= 12) THEN
       CALL sgcoor10(
     1   xx,        yy,        zz,        x,
     2   xdp,       xx0,       yy0,       zz0,
     3   vx0,       vy0,       vz0,       gbuf%SMSTR,
     4   nc,        d,         gbuf%OFF,  offg0,
     5   nel,       mtn,       ismstr)
C----------------------
       IF (ismstr ==  11) THEN
        CALL s10derit3(
     1   volp,       deltax,     deltax2,    xx0,
     2   yy0,        zz0,        px,         py,
     3   pz,         nx,         rx,         ry,
     4   rz,         sx,         sy,         sz,
     5   tx,         ty,         tz,         wip(1,npt),
     6   alph(1,npt),beta(1,npt),voln,       volg,
     7   voldp,      nel,        gbuf%OFF,   npt)
        CALL s10len3(
     1   volp,         ngl,          deltax,       deltax2,
     2   px,           py,           pz,           volg,
     3   gbuf%VOL,     rx,           ry,           rz,
     4   sx,           sy,           sz,           tx,
     5   ty,           tz,           nc,           nel,
     6   mxt,          pm,           gbuf%ISMS,    gbuf%DT_PITER,
     7   npt,          iint,         isrot,        iformdt)
       ELSE
C!!!!!!calcul local rep for ISMSTR 10 to 11 (offg>un) 
        IF (ismstr == 12.AND.ism12_11==0.AND.idtmin(1)==3) THEN
         CALL s10rcoor12(
     1   gbuf%OFF,x,       xdp,     nc,
     2   e1x,     e2x,     e3x,     e1y,
     3   e2y,     e3y,     e1z,     e2z,
     4   e3z,     nel)
        END IF !(ISMSTR == 10.AND.ISORTH/=0.AND.IDTMIN(1)==3)
        ibid = 1
        DO ip=1,npt
         lbuf => elbuf_tab(ng)%BUFLY(ibid)%LBUF(ip,ibid,ibid)
         CALL s10pijto3(px(1,1,ip),py(1,1,ip),pz(1,1,ip),lbuf%PIJ,llt)
        ENDDO
       END IF !(ISMSTR ==  11) THEN
       DO ip=1,npt
        CALL s10defot3(
     1   px(1,1,ip),py(1,1,ip),pz(1,1,ip),vx0,
     2   vy0,       vz0,       mfxx(1,ip),mfxy(1,ip),
     3   mfxz(1,ip),mfyx(1,ip),mfyy(1,ip),mfyz(1,ip),
     4   mfzx(1,ip),mfzy(1,ip),mfzz(1,ip),nel)
       END DO 
       IF (ismstr == 12.AND.ism12_11==0.AND.idtmin(1)==3) THEN
         DO ip=1,npt
          CALL sordeft12(lft,llt,mfxx(1,ip), mfxy(1,ip), mfxz(1,ip),
     .         mfyx(1,ip), mfyy(1,ip), mfyz(1,ip),
     .         mfzx(1,ip), mfzy(1,ip), mfzz(1,ip),
     .         e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z,gbuf%OFF)
         END DO 
       ENDIF
      ENDIF     ! ISMSTR == 10
C-----------
      IF (ismstr /= 11) THEN
      CALL s10deri3(
     1   off,        volp,       ngl,        deltax,
     2   deltax2,    xx,         yy,         zz,
     3   px,         py,         pz,         nx,
     4   rx,         ry,         rz,         sx,
     5   sy,         sz,         tx,         ty,
     6   tz,         wip(1,npt), alph(1,npt),beta(1,npt),
     7   voln,       volg,       voldp,      nc,
     8   gbuf%SMSTR, gbuf%OFF,   nel,        npt,
     9   ismstr,     jlag)
C
        CALL s10len3(
     1   volp,         ngl,          deltax,       deltax2,
     2   px,           py,           pz,           volg,
     3   gbuf%VOL,     rx,           ry,           rz,
     4   sx,           sy,           sz,           tx,
     5   ty,           tz,           nc,           nel,
     6   mxt,          pm,           gbuf%ISMS,    gbuf%DT_PITER,
     7   npt,          iint,         isrot,        iformdt)
 
        IF (iplaw1>0) THEN
           CALL s10dvm12(
     .               px, py, pz, vx, vy, vz, 
     .               dvm ,gbuf%OFF, npt ,nel)
        END IF
C--- for /CST avoid multi-printout----
       IF (ismstr == 2 .OR.ismstr ==12) THEN
        DO ip=1,npt
         iptr = ip
         lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(iptr,ipts,iptt)
         DO i=lft,llt
          IF (gbuf%OFF(i)==two) lbuf%OFF(i)=gbuf%OFF(i)
         ENDDO
        ENDDO
       END IF 
      END IF 
C --------------------------
C  --- UPDATE REF. CONFIGURATION 
C --------------------------
      IF (ismstr <= 3.OR.(ismstr==4.AND.jlag>0)) THEN
        CALL s10sav3(
     1   gbuf%OFF,  gbuf%SMSTR,nc,        xx,
     2   yy,        zz,        nel)
      END IF !(ISMSTR == 2) THEN
C-----------------------------
      IF (isorth == 0) THEN            
        DO i=lft,llt                                            
          gama(i,1) = one                               
          gama(i,2) = zero                                
          gama(i,3) = zero             
          gama(i,4) = zero                                
          gama(i,5) = one                                
          gama(i,6) = zero             
        ENDDO                          
      ELSE                             
        CALL sreploc3(
     1   rx,      ry,      rz,      sx,
     2   sy,      sz,      tx,      ty,
     3   tz,      e1x,     e2x,     e3x,
     4   e1y,     e2y,     e3y,     e1z,
     5   e2z,     e3z,     llt)
        CALL sorthdir3(
     1   rx,       ry,       rz,       sx,
     2   sy,       sz,       tx,       ty,
     3   tz,       e1x,      e2x,      e3x,
     4   e1y,      e2y,      e3y,      e1z,
     5   e2z,      e3z,      gbuf%GAMA,gama,
     6   nel,      irep)
      ENDIF
      IF(icp >0 .AND. ismstr/=10) THEN
         DO i=lft,llt
          IF(gbuf%OFF(i) == zero) cycle
          sum=sfem_nodvar(nc(i,1))+sfem_nodvar(nc(i,2))+sfem_nodvar(nc(i,3))+sfem_nodvar(nc(i,4))
          jacgm(i)=fourth*sum
         ENDDO
      ENDIF
C-----------------------------
C     POINTS D' INTEGRATION 
C-----------------------------
      DO ip=1,npt
        iptr = ip
        lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(iptr,ipts,iptt)
       IF (ioffs == 1)THEN
        DO i=lft,llt
         IF (offs(i)<=two) lbuf%OFF(i)=offs(i)
        ENDDO
       END IF
C
        CALL s10defo3(
     1   px(1,1,ip),py(1,1,ip),pz(1,1,ip),vx,
     2   vy,        vz,        dxx,       dxy,
     3   dxz,       dyx,       dyy,       dyz,
     4   dzx,       dzy,       dzz,       d4,
     5   d5,        d6,        wxx,       wyy,
     6   wzz,       volp(1,ip),voln,      lbuf%RHO,
     7   rhoo,      nel,       jhbe,      isrot)
 
         IF (ismstr == 12.AND.ism12_11==0.AND.idtmin(1)==3) THEN
          CALL sordef12(lft,llt,dxx, dyy, dzz,
     .         d4, d5, d6,
     .         e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z,offg0)
         ENDIF
         IF (icp>0) THEN
           IF (ismstr==10) THEN
             DO i=lft,llt
               IF(gbuf%OFF(i) == zero) cycle
               jacgm(i)=sfem_nodvar(nc(i,ip))
             ENDDO
           END IF
           CALL s10_icp(
     1   gbuf%OFF,   jacgm,      facp,       nu,
     2   mfxx(1,ip), mfxy(1,ip), mfxz(1,ip), mfyx(1,ip),
     3   mfyy(1,ip), mfyz(1,ip), mfzx(1,ip), mfzy(1,ip),
     4   mfzz(1,ip), lbuf%VOL,   voln,       lbuf%VOL0DP,
     5   voldp(1,ip),nel,        ismstr)
         ENDIF
C   
        divde(1:nel) = dt1*(dxx(1:nel)+ dyy(1:nel)+ dzz(1:nel))  
        IF (iplaw1>0) CALL s10divde12(dvm ,divde ,facp,gbuf%OFF,nel)
        CALL srho3(
     1   pm,         lbuf%VOL,   lbuf%RHO,   lbuf%EINT,
     2   divde,      flux(1,nf1),flu1(nf1),  voln,
     3   dvol,       ngl,        mxt,        off,
     4   0,          gbuf%TAG22, voldp(1,ip),lbuf%VOL0DP,
     5   amu,        gbuf%OFF,   nel,        mtn,
     6   jale,       ismstr,     jeul,       jlag)
c
       IF (ismstr == 12.AND.ism12_11==0.AND.idtmin(1)==3) THEN
          CALL sroto12_sig(lft,llt,lbuf%SIG,nel,
     .                     e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z,offg0)
!! temporary replaced by (the same) SROTO12_SIG() in order to do not affect
!! the other multidimensional buffer ARRAYS which are still not modified
!!          CALL SROTO12(LFT,LLT,LBUF%SIG,
!!     .                 E1X, E1Y, E1Z, E2X, E2Y, E2Z, E3X, E3Y, E3Z,OFFG0)
       ENDIF
        CALL srota3(
     1   lbuf%SIG,s1,      s2,      s3,
     2   s4,      s5,      s6,      wxx,
     3   wyy,     wzz,     nel,     mtn,
     4   ismstr)
C-----------------------------
C       SMALL STRAIN
C-----------------------------
        CALL s10malla3(
     1   gbuf%OFF,   off,        wxx,        wyy,
     2   wzz,        wxxg,       wyyg,       wzzg,
     3   wip(ip,npt),nel,        ismstr,     jlag)
           
        IF(jthe < 0 ) THEN                                        
          DO i=lft,llt                                            
            tempel(i)= zero                                       
          ENDDO                                                   
          IF(isolnod == 10) THEN      
          DO j = 1,10                                             
            DO i=lft,llt                                          
              tempel(i)= tempel(i) +  nx(i,j,ip)*temp(nc(i,j))       
            ENDDO                                                 
          ENDDO                                                   
          ELSEIF(isolnod == 4) THEN
            DO j = 1,4                                             
              DO i=lft,llt                                          
                tempel(i)= tempel(i) +  nxt4(i,j,ip)*temp(nc(i,j))
              ENDDO 
            ENDDO 
        ENDIF                                                     
        ENDIF                                                     
C------------------------------------------------------
C     CALCUL DES CONTRAINTES SUIVANT LOIS CONSTITUTIVES
C------------------------------------------------------
      IF ((itask==0).AND.(imon_mat==1)) CALL startime(timers,35)
C
      IF(iboltp /= 0) CALL boltst(
     .                ip,        bpreld,    lbuf%SIG,tt,        
     .                nel   ,npt   ,sensors%NSENSOR,sensors%SENSOR_TAB)
C
        CALL mmain(timers, output,
     1   elbuf_tab,   ng,          pm,          geo,
     2                ale_connect, ixs,         iparg,
     3   v,           tf,          npf,         bufmat,
     4   sti,         x,           dt2t,        neltst,
     5   ityptst,     offset,      nel,         w,
     6   off,         ngeo,        mxt,         ngl,
     7   voln,        vd2,         dvol,        deltax,
     8   vis,         qvis,        cxx,         s1,
     9   s2,          s3,          s4,          s5,
     a   s6,          dxx,         dyy,         dzz,
     b   d4,          d5,          d6,          wxx,
     c   wyy,         wzz,         rx,          ry,
     d   rz,          sx,          sy,          sz,
     e   vdx,         vdy,         vdz,         muvoid,
     f   ssp_eq,      aire,        sigy,        et,
     g   r1_free,     lbuf%PLA,    r3_free,     amu,
     h   mfxx(1,ip),  mfxy(1,ip),  mfxz(1,ip),  mfyx(1,ip),
     i   mfyy(1,ip),  mfyz(1,ip),  mfzx(1,ip),  mfzy(1,ip),
     j   mfzz(1,ip),  ipm,         gama,        bid,
     k   bid,         bid,         bid,         bid,
     l   bid,         bid,         istrain,     tempel,
     m   die,         iexpan,      ilay,        mssa,
     n   dmels,       iptr,        ipts,        iptt,
     o   table,       bid,         bid,         bid,
     p   bid,         iparg(1,ng), igeo,        conde,
     q   itask,       nloc_dmg,    varnl,       mat_elem,
     r   h3d_strain,  jplasol,     jsph,        mvsiz,
     *   snpc,        stf,         sbufmat,     glob_therm, 
     s   svis,        sz_ix,       iresp,
     t   n2d,         th_strain,   ngroup,      tt, 
     .   dt1,         ntable,      numelq,      nummat,
     .   numgeo,      numnod,      numels,
     .   idel7nok,    idtmin,      maxfunc,
     .   imon_mat,    userl_avail, impl_s,
     .   idyna,       dt,          fheat      ,sensors)
C
      IF ((itask==0).AND.(imon_mat==1)) CALL stoptime(timers,35)
C-----------  return to global system    
       IF (ismstr == 12.AND.ism12_11==0.AND.idtmin(1)==3) THEN
         CALL sroto12_sig(lft,llt,lbuf%SIG,nel,
     .                 e1x,e2x,e3x,e1y,e2y,e3y,e1z,e2z,e3z,offg0)
!! temporary replaced by (the same) SROTO12_SIG() in order to do not affect
!! the other multidimensional buffer ARRAYS which are still not modified
!!         CALL SROTO12(LFT,LLT,LBUF%SIG,
!!     .                 E1X,E2X,E3X,E1Y,E2Y,E3Y,E1Z,E2Z,E3Z,OFFG0)
         IF (istrain == 1) THEN 
          CALL sordef12(lft,llt,dxx, dxy, dxz,
     .         d4, d5, d6,
     .         e1x,e2x,e3x,e1y,e2y,e3y,e1z,e2z,e3z,offg0)
         ENDIF
       ENDIF
        IF (istrain == 1) CALL sstra3(
     1   dxx,      dyy,      dzz,      d4,
     2   d5,       d6,       lbuf%STRA,wxx,
     3   wyy,      wzz,      off,      nel,
     4   jcvt)
C--------------------------
        iflag=mod(ncycle,ncpri)
        IF(ioutprt>0)THEN
           CALL s10bilan(partsav,lbuf%EINT,lbuf%RHO,lbuf%RK,lbuf%VOL,
     .                   vx, vy, vz,nx(1,1,ip),voln,iparts,
     .                   gresav,grth,igrth,iexpan,lbuf%EINTTH,
     .                   gbuf%FILL,xx,yy,zz,itask,iparg(1,ng),gbuf%OFF,sensors,
     .                   nel, elbuf_tab(ng)%BUFLY(ilay)%L_WPLA, lbuf%WPLA)
        ENDIF
C-------------------------
        IF (cns2>zero)
     .     CALL nsvis_sm12(visp  ,cns2,cxx  ,voln ,dxx     ,
     .                     dyy     ,dzz    ,d4    ,d5  ,d6   ,
     .                     lbuf%VOL,rho0_1,sti   ,nel  ,svis ) 
C----------------------------
C       INTERNAL FORCES
C----------------------------
        CALL s10fint3(
     1   lbuf%SIG,   px(1,1,ip), py(1,1,ip), pz(1,1,ip),
     2   fx,         fy,         fz,         voln,
     3   qvis,       sti,        stig,       lbuf%EINT,
     4   lbuf%RHO,   lbuf%QVIS,  lbuf%PLA,   lbuf%EPSD,
     5   gbuf%EPSD,  gbuf%SIG,   gbuf%EINT,  gbuf%RHO,
     6   gbuf%QVIS,  gbuf%PLA,   wip(ip,npt),gbuf%G_PLA,
     7   nel,        conde,      condeg,     gbuf%G_EPSD,
     8   israt,      svis ,glob_therm%NODADT_THERM)
C
        DO i=lft,llt
         IF (lbuf%OFF(i) > one .AND. gbuf%OFF(i) == one) THEN
C          switched to small strain
           offs(i)=min(lbuf%OFF(i),offs(i))
           ioffs  =1
         END IF
        ENDDO
C
        IF (jthe < 0 .AND. isolnod == 10) THEN
          imat = mxt(1)
          IF (mat_elem%MAT_PARAM(imat)%HEAT_FLAG == 1) THEN
            CALL s10therm(
     1           pm,        imat,       nc,        voln,
     2           px(1,1,ip),py(1,1,ip),pz(1,1,ip),nx(1,1,ip),
     3           dt1,       temp,      tempel,    fheat,
     4           them,      gbuf%OFF,  lbuf%OFF,  nel,
     5           glob_therm%THEACCFACT)
          ELSE
            CALL s10therm(
     1           pm,        imat,       nc,        voln,
     2           px(1,1,ip),py(1,1,ip),pz(1,1,ip),nx(1,1,ip),
     3           dt1,       temp,      tempel,    die,
     4           them,      gbuf%OFF,  lbuf%OFF,  nel,
     5           glob_therm%THEACCFACT)
          END IF
        ENDIF
c----------
      ENDDO   ! end integration point loop
 
      IF (jthe < 0 .AND. isolnod == 4) THEN
          imat = mxt(1)
          IF (mat_elem%MAT_PARAM(imat)%HEAT_FLAG == 1) THEN
            CALL s4therm_itet1(pm   ,imat   ,nc   ,nel  ,
     .                         xx   ,yy     ,zz   ,dt1  ,fheat ,
     .                         temp ,them   ,gbuf%OFF   ,lbuf%OFF,
     .                         glob_therm%THEACCFACT)
          ELSE
            CALL s4therm_itet1(pm   ,imat   ,nc   ,nel  ,
     .                         xx   ,yy     ,zz   ,dt1  ,die ,
     .                         temp ,them   ,gbuf%OFF   ,lbuf%OFF,
     .                         glob_therm%THEACCFACT)
          END IF
      ENDIF
c-----------------------------
      IF (jlag+jale+jeul /= 0) THEN
c
C-----  small strain
C------ correction of GBUF%SMSTR will be done only from next cycle
        CALL s10mallb3(
     1   gbuf%SMSTR,gbuf%OFF,  wxxg,      wyyg,
     2   wzzg,      nel,       ismstr,    jlag)
        IF (ioffs == 1)THEN
          DO i=lft,llt
C           switch to small strain
            IF (offs(i)<=two) gbuf%OFF(i) = offs(i)
          END DO
C
        ipts = 1
        iptt = 1
        ilay = 1
          DO ip=1,npt
           iptr = ip
           lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(iptr,ipts,iptt)
            DO i=lft,llt
              IF (gbuf%OFF(i) > one) lbuf%OFF(i)=gbuf%OFF(i)
            END DO
          END DO
        END IF
c
        itet=1
        CALL s10mallgeo3(ngl,gbuf%OFF ,volg ,deltax, gbuf%VOL ,
     .                       rx , ry , rz , 
     .                       sx , sy , sz ,
     .                       tx , ty , tz ,deltax4,geo(1,ngeo(1)),
     .                       nel,npt,ismstr,isrot,dt)
        rbid(lft:llt)=zero
        CALL sgeodel3(ngl,gbuf%OFF,volg,deltax4,gbuf%VOL,geo(1,ngeo(1)),rbid,dt,nel,idel7nok)
        CALL smallb3(gbuf%OFF,off,nel,ismstr)
        CALL smallgeo3(ngl, gbuf%OFF ,volg ,deltax4, gbuf%VOL ,itet, nel, ismstr,dt)
C-----------!!!-seperate LBUF%SIGL=LBUF%SIG,w/ npt>0 ISM12_11
       IF (ismstr == 12.AND.idtmin(1)==3) THEN
        ioffs =0
         DO i=lft,llt
          IF(gbuf%OFF(i)/=offg0(i).AND.abs(gbuf%OFF(i)) > one ) ioffs=1
         ENDDO
         IF (ioffs == 1) THEN
           CALL s10sav12(
     1   gbuf%OFF,  offg0,     gbuf%SMSTR,nc,
     2   xx,        yy,        zz,        nel)
           IF (ism12_11>0 .AND. isorth == 0) THEN
             CALL s10upd11t12(
     1   elbuf_tab(ng),gbuf%OFF,     offg0,        nc,
     2   xx,           yy,           zz,           nel,
     3   npt)
           END IF 
          ipts = 1
          iptt = 1
          ilay = 1
          DO ip=1,npt
           iptr = ip
           lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(iptr,ipts,iptt)
            DO i=lft,llt
              IF (abs(gbuf%OFF(i)) > one) lbuf%OFF(i)=gbuf%OFF(i)
            END DO
          END DO
         END IF !(IOFFS == 1) THEN
       END IF 
C----------------------------
C     distortion control  
C----------------------------
      IF (isctl > 0) THEN
         offg(1:nel) = min(off(1:nel),abs(gbuf%OFF(1:nel))) ! no offg for tetra10 
         CALL sdistor_ini(                                      
     1                nel      ,sti_c    ,npropm     ,nummat  ,       
     2                ismstr   ,mxt      ,istab      ,pm      ,       
     3                gbuf%SIG ,gbuf%RHO ,cxx        ,offg    ,       
     4                gbuf%OFF ,ll       ,voln       ,fld     ,       
     5                cns2     ,fqmax    )
         IF (ismstr>=11)  ! XX,YY,ZZ now current conf.
     *     CALL s10get_x3(
     1          x,     xdp,      dr,    numnod,
     2         xx,      yy,      zz,        nc,
     3      isrot,   iresp,    nel )
         IF (ismstr<10)  ! remove limitation of Ismstr>=
     *     CALL s10get_x0(
     1        nel,  numnod,       x,     xdp,         
     2          d,     xx0,     yy0,     zz0,
     3         nc)
         CALL s10for_distor(
     .         stig,   fld  ,   sti_c,
     .          xx ,     yy ,     zz ,   
     .          vx ,     vy ,     vz ,     
     .          fx ,     fy ,     fz ,     
     .          xx0,     yy0,     zz0,
     .         cns2,   istab,    ll  ,
     .        fqmax,    nel ,gbuf%EINT_DISTOR,
     .          dt1)
      ENDIF
c-----------------------------
      IF(nfilsol/=0) CALL sxfillopt(
     1   npe,      gbuf%FILL,stig,     fx,
     2   fy,       fz,       nel)
C-----
        IF (iparit == 0) THEN
          CALL s10cumu3(
     1   gbuf%OFF,  a,         nc,        stifn,
     2   stig,      fx,        fy,        fz,
     3   deltax2,   them,      fthe,      ar,
     4   x,         stifr,     gbuf%SMSTR,condn,
     5   condeg,    itagdn,    nel,       ismstr,
     6   jthe,      isrot ,glob_therm%NODADT_THERM)
          ELSE
          CALL s10cumu3p(
     1   gbuf%OFF,  stig,      fsky,      fsky,
     2   iads,      fx,        fy,        fz,
     3   deltax2,   iads10,    nc,        them,
     4   fthesky,   ar,        x,         gbuf%SMSTR,
     5   condnsky,  condeg,    itagdn,    nel,
     6   nft,       ismstr,    jthe,      isrot,glob_therm%NODADT_THERM)
        ENDIF
c-----
      ENDIF
c-----------
      RETURN
      END