initwg_tri.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/.
!||====================================================================
!||    initwg_tri         ../starter/source/spmd/domain_decomposition/initwg_tri.F
!||--- called by ------------------------------------------------------
!||    initwg             ../starter/source/spmd/domain_decomposition/initwg.F
!||--- calls      -----------------------------------------------------
!||    bidon              ../starter/source/system/machine.F
!||    interlagran        ../starter/source/spmd/domain_decomposition/grid2mat.f
!||--- uses       -----------------------------------------------------
!||    ddweights_mod      ../starter/share/modules1/ddweights_mod.F
!||    mid_pid_mod        ../starter/share/modules1/mid_pid_mod.F
!||====================================================================
      SUBROUTINE initwg_tri(WD,PM,GEO,IXTG,IGEO,NUMELTG,IPM ,SIZE_IRUP,
     .            NUMMAT,NUMGEO,POIN_PART_TRI,MID_PID_TRI,IPARTG,
     .            OFF,BUFMAT,MID_OLD,PID_OLD,MLN_OLD,RECHERCHE,TELT_PRO,
     .            TABMP_L,MAT_PARAM)
C-----------------------------------------------
C            M o d u l e s
C-----------------------------------------------
      USE ddweights_mod
      USE mid_pid_mod
      USE matparam_def_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "param_c.inc"
#include      "com01_c.inc"
#include      "tablen_c.inc"
#include      "ddspmd_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  OFF,NUMELTG,TABMP_L,
     .  NUMMAT,NUMGEO, IXTG(NIXTG,*),IGEO(NPROPGI,*),
     .  IPM(NPROPMI,*)
      INTEGER, INTENT(IN) :: SIZE_IRUP
 
C     REAL OU REAL*8
      my_real
     .    pm(npropm,*), geo(npropg,*),bufmat(*)
 
      REAL WD(*)
      INTEGER MID_OLD,PID_OLD,MLN_OLD,RECHERCHE
      my_real TELT_PRO
 
      INTEGER, DIMENSION(*), INTENT(IN) :: IPARTG
      INTEGER, DIMENSION(2,*), INTENT(IN) :: POIN_PART_TRI
      TYPE(MID_PID_TYPE), DIMENSION(*), INTENT(INOUT) :: MID_PID_TRI
      TYPE(matparam_struct_) ,DIMENSION(NUMMAT), INTENT(IN) :: MAT_PARAM
C-----------------------------------------------
      INTEGER  NPN, MID, PID, JHBE, IGT, MLN,
     .    istrain, ithk, ihbe, ipla, issn, mtn, i, j, k,l,
     .    nfunc,mpt,npts,nptt,nptr,nptot,iflag,jsrot,
     .    i_mid,i_pid,i_mid_old,i_pid_old,puid,muid,
     .    elm_typ,elm_typ_old,ilaw,ilaw_old,test_mat,
     .    i_pro,isol2,muid_old,puid_old,
     .    test,nfunc1,nfunc2,nfail,irup2,ii,irup_tab(size_irup),
     .    indi,iad,indi2,mult
      INTEGER :: INDI3,COST_CHECK,POIN_PART,POIN_MID,POIN_PID
      INTEGER :: FLAG_NICE_NEWTON,FLAG_GURSON,FLAG_NON_LOCAL
      INTEGER :: SPECIAL_OPTION,SPE_I_1,SPE_I_2,SPE_I_3
     
      real
     .   wtype(9),fwihbe,fac8,
     .   tabmat(3),tabx(3),timmat,npt,telt,poids,w,
     .   batozmult,tmat,trup,tabrup(3),trup_local,tmatadd,
     .   wd_local,mult_spe,visc_prony
      
      my_real
     .        cc, a,b,a1,a2
      my_real
     .       invtref
       DATA wtype /1.6 ,1. ,1. ,.9 ,1.1 ,1.4 ,0.65 ,.9 ,2.0/     
C-----------------------------------------------
      CALL bidon()
      telt = 0
      nfunc = 0
!     DD_OPTIMIZATION = 0 --> default case, DD optimized for Broadwell processor
!     DD_OPTIMIZATION = 1 --> DD optimized for Skylake processor 
!     DD_OPTIMIZATION = 2 --> DD optimized for Sandy Bridge processor   
!     DD_OPTIMIZATION = 3 --> default case for ARM processor, DD optimized for ThunderX2 processor (ARM)   
      IF(dd_optimization==1) THEN
!       Skylake processor 
#include "weights_p4linux964_spmd_avx512.inc"
      ELSEIF(dd_optimization==2) THEN
!       Sandy Bridge processor   
#include "weights_p4linux964_spmd_sse3.inc"
      ELSEIF(dd_optimization==3) THEN
!       ThunderX2 processor (ARMV8.0)      
#include "weights_p4linuxa964_spmd.inc"
      ELSE
!       DEFAULT CASE
#if ARCH_CPU
!       ThunderX2 processor (ARMV8.0)      
#include "weights_p4linuxa964_spmd.inc"
#elif 1
!       Broadwell processor
#include "weights_p4linux964_spmd.inc"
#endif
      ENDIF
      invtref = one/tpsref
      indi3 = 2
      DO i = 1, numeltg
        mid= ixtg(1,i)
        pid= ixtg(5,i)
        mln = nint(pm(19,abs(mid)))
        wd_local = wd(i+off)
        ! -----------------
        IF(recherche==1) THEN
         mid = mid_old
         pid = pid_old
         mln = mln_old
         wd_local = zero
        ENDIF
        ! -----------------
        npn = nint(geo(6,pid))
        ihbe = nint(geo(171,pid))
        ithk = nint(geo(35,pid))
        ipla = nint(geo(39,pid))
        npt = max(abs(npn),1)
        flag_non_local = 0  
        special_option = 0
        spe_i_1 = 1
        spe_i_2 = 1
        nfail = mat_param(mid)%NFAIL
        irup_tab = 0
        IF(nfail/=0) THEN          ! up to 6 failure models per material
         DO j=1,nfail
          irup_tab(j) = mat_param(mid)%FAIL(j)%IRUPT
         ENDDO
        ENDIF
        timmat = 0.
        trup = 0.
        tmatadd = 0.
        visc_prony = 0.
        mult = 0
        IF((mln<28).OR.(mln==32)) THEN
         irup2 = 0
        ELSE
         irup2 = 3
        ENDIF
 
        IF (mln==2.OR.mln==3) THEN
          cc = pm(43,mid)
          IF (cc/=0) THEN
           indi = 2
          ELSE
           indi = 1
          ENDIF
          IF (mat_param(abs(mid))%IVISC > 0) THEN 
            visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
          ENDIF
        ELSEIF ((mln==25).AND.(abs(npn)>0)) THEN
          iflag = nint( pm(40,mid))
          IF (iflag/=0) THEN
           indi = 2
          ELSE
           indi = 1
          ENDIF
           IF (mat_param(abs(mid))%IVISC > 0) THEN 
             visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
           ENDIF
C loi 36+86 en fonction des fct
        ELSEIF (mln==36.OR.(mln==86).AND.(abs(npn)>0)) THEN
          nfunc = nint(pm(40,mid))
          IF (nfunc<=2) THEN
           indi = 1
          ELSEIF (nfunc>2.AND.nfunc<=7) THEN
           indi = 2
          ELSEIF (nfunc>7) THEN
           indi = 3
          ENDIF
           IF (mat_param(abs(mid))%IVISC > 0) THEN 
             visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
           ENDIF
        ELSEIF((mln==42).OR.(mln==62).OR.(mln==69)) THEN  ! check prony option
           nfunc = 0   ! NPRONY
           IF (mln==42) nfunc = mat_param(abs(mid))%IPARAM(2)
           IF (mln==62) THEN
             iad=ipm(7,abs(mid))-1
             nfunc = nint(bufmat(iad+3))
           END IF
 
           IF(nfunc==0) THEN
            indi = 1
             IF (mat_param(abs(mid))%IVISC > 0) THEN 
               visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
             ENDIF
           ELSEIF(nfunc==1) THEN
            indi = 2
           ELSEIF(nfunc==2) THEN
            indi = 3
           ELSEIF(nfunc>2) THEN
            indi = 3
            mult = nfunc - 2
            indi2 = 2            
           ENDIF
        ELSEIF((mln==82)) THEN
           iad=ipm(7,abs(mid))-1
           nfunc=nint(bufmat(iad+1))
           IF(nfunc<=1) THEN
            indi = 1
            IF (mat_param(abs(mid))%IVISC > 0) THEN 
              visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
            ENDIF
           ELSEIF(nfunc==2) THEN
            indi = 2
           ELSEIF(nfunc==3) THEN
            indi = 3
           ELSEIF(nfunc>3) THEN
            indi = 3
            mult = nfunc - 3
            indi2 = 2            
           ENDIF
        ELSEIF(mln==104) THEN
           iad=ipm(7,abs(mid))-1
           flag_nice_newton=nint(bufmat(iad+11))
           IF(flag_nice_newton==2) THEN !   Newtow algo
             indi = 2
           ELSE     !   nice algo
             indi = 1
           ENDIF
           flag_gurson=nint(bufmat(iad+30))
           IF(flag_gurson/=0) THEN
             special_option=1
             spe_i_1 = 1
             spe_i_2 = 1
           ENDIF
           IF(flag_gurson==1) THEN
             spe_i_2 = 1
           ELSEIF(flag_gurson==2) THEN
             spe_i_2 = 2
           ELSEIF(flag_gurson==3) THEN
             spe_i_2 = 3
           ENDIF           
           flag_non_local = mat_param(abs(mid))%NLOC
        ELSE
         indi = 1
         IF (mat_param(abs(mid))%IVISC > 0) THEN 
           visc_prony = visc_prony_cost * mat_param(abs(mid))%VISC%IPARAM(1)
         ENDIF
        ENDIF 
 
        mult_spe = 0.
        spe_i_3 = 1
        IF(flag_non_local/=0) THEN
            spe_i_3 = 2
            mult_spe = npt
        ENDIF
 
        cost_check = 0
!****************************************
!       ---------------------------
!       SHELL3N ELEMENT
!       ---------------------------
        ! check if the (mid,pid) cost must be initialized from a previous run
        IF(recherche==0.AND.test_poids/=0) THEN
          poin_part = ipartg(i)
          poin_mid = poin_part_tri(1,poin_part)
          poin_pid = poin_part_tri(2,poin_part)
        ! if POIN_MID==0 and POIN_PID == 0, the element cost in the .ddw file is 0 --> must be initialized 
        ! from the .inc file 
          IF(poin_mid/=0.AND.poin_pid/=0) THEN
           IF(mid_pid_tri(poin_mid)%COST1D(poin_pid)/=zero) THEN
                cost_check = 1
                telt = mid_pid_tri(poin_mid)%COST1D(poin_pid)
           ENDIF
          ENDIF
        ENDIF
        ! the (mid,pid) cost must be initialized from .inc file
        IF(cost_check==0) THEN
         IF( ddweights(1,2,mid)/=zero)THEN
!       Compute time according to integration points:
            a1 = ddweights(1,2,mid) * tpsref 
            a2 = ddweights(2,2,mid) * tpsref 
 
            IF (a2 /=zero)THEN
!       Compute line function A1 ="Time for 1 int Point" - A2="time for 5 int points"
              a = (a2-a1)/4
              b = a1-a
              timmat = a*npt + b
            ELSE
              timmat = a1*npt
            ENDIF
!           --------------
!           Failure
            IF(nfail/=0) THEN
             DO j=1,nfail
              a1 = rupture_shell(irup_tab(j),irup2+1)
              a2 = rupture_shell(irup_tab(j),irup2+3)
              IF (a2 /=zero)THEN
               a = (a2-a1)/4
               b = a1-a
               trup = trup + a*npt + b
              ELSE
               trup = trup + a1*npt
              ENDIF
             ENDDO
            ENDIF
!           --------------
         ELSE
          IF(ithk==2)THEN
            ithk = 0
          ELSEIF(mln==32)THEN
            ithk = 1
          ENDIF
          istrain = nint(geo(11,pid))
          IF(mln==19.OR.mln>=25)istrain = 1
          issn = nint(geo(3,pid))
!       ow test elem delete
          IF (wd_local==0.) THEN
           IF(abs(npn)>0) THEN
              tabx(1) = 1.
              tabx(2) = 3.
              tabx(3) = 5.
              ! *******--------*******
              DO j=1,3
                IF(mult/=0) tmatadd = mult * 
     .              (tritnl(min(mln,maxlaw),j,indi) - tritnl(min(mln,maxlaw),j,indi2) )
                  IF(special_option/=0) tmatadd = tmatadd + shtnl_option(spe_i_1,spe_i_2) 
                tabmat(j) = tritnl(min(mln,maxlaw),j,indi) + tmatadd
              ENDDO
              ! *******--------*******
              npt = abs(npn)
              CALL interlagran(tabmat,tabx,3,npt,timmat)
!              ----------------
!               Failure 
                IF(nfail/=0) THEN
                 DO j=1,nfail
                  DO ii=1,3
                   tabrup(ii) = rupture_shell(irup_tab(j),irup2+ii)
                  ENDDO
                  CALL interlagran(tabrup,tabx,3,npt,trup_local)
                  trup = trup + trup_local
                 ENDDO
                ENDIF          !   <--- fin nfail/=0
!               ---------------- 
 
            ELSE
!       0pt d integration doit etre traite a part
              ! *******--------*******
              IF(mult/=0) tmatadd = mult * 
     .              (tritnl(min(mln,maxlaw),0,indi) - tritnl(min(mln,maxlaw),0,indi2) )
              timmat =  tritnl(min(mln,maxlaw),0,indi)  + tmatadd
              ! *******--------*******
!             ----------------
!             Failure incompatible N=0
              IF(nfail/=0) THEN
               trup = 0.
              ENDIF          !   <--- fin NFAIL/=0
!             ----------------  
            ENDIF
          ENDIF
         ENDIF
         IF(mln/=0)THEN
              telt = tritelt(1) 
         ENDIF
        ENDIF   !   fin TRITNL_OLD(MID,PID)
!****************************************
 
        ! --------------------
        IF(recherche==0) THEN
         IF((wd_local==0.).AND.(mln/=0))THEN
          poids = (telt + timmat + trup + mult_spe*nlocal_option(spe_i_3) + visc_prony) * invtref
          wd(i+off) = poids
 
          poin_part = ipartg(i)
          poin_mid = poin_part_tri(1,poin_part)
          poin_pid = poin_part_tri(2,poin_part)
          IF(poin_mid/=0.AND.poin_pid/=0) 
     .    mid_pid_tri(poin_mid)%COST1D(poin_pid) = telt + timmat + trup +
     .           mult_spe*nlocal_option(spe_i_3)
         ELSE
          wd(i+off) = 0.0001
         END IF
        ELSE
         telt_pro = telt + timmat + trup + mult_spe*nlocal_option(spe_i_3)
        ENDIF
        ! --------------------
      ENDDO
      RETURN
      END