spgrtails.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/.
!||====================================================================
!||    spgrtails              ../starter/source/elements/sph/spgrtails.F
!||--- called by ------------------------------------------------------
!||    lectur                 ../starter/source/starter/lectur.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                 ../starter/source/output/message/message.F
!||    fretitl2               ../starter/source/starter/freform.F
!||    get_u_geo              ../starter/source/user_interface/uaccess.F
!||    zeroin                 ../starter/source/system/zeroin.F
!||--- uses       -----------------------------------------------------
!||    message_mod            ../starter/share/message_module/message_mod.F
!||    r2r_mod                ../starter/share/modules1/r2r_mod.F
!||====================================================================
      SUBROUTINE spgrtails(KXSP   , IPARG , PM  ,IPART ,
     2                     IPARTSP, EADD,   ND  , CEPSP, DD_IAD,
     3                     IDX,   IXSP, IPM  , IGEO,
     4                     SPBUF,SPH2SOL,SOL2SPH,
     5                     IRST   , NOD2SP ,PRINT_FLAG,MAT_PARAM ,
     6                     IXSPS)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE my_alloc_mod
      USE message_mod
      USE r2r_mod
      USE matparam_def_mod
      USE names_and_titles_mod , ONLY : nchartitle
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      "vect01_c.inc"
#include      "com01_c.inc"
#include      "com04_c.inc"
#include      "sphcom.inc"
#include      "units_c.inc"
#include      "param_c.inc"
#include      "scr17_c.inc"
#include      "r2r_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, DIMENSION(KVOISPH,NUMSPH),INTENT(INOUT) :: IXSPS
      INTEGER IDX,ND,
     .        KXSP(NISP,*),IPARG(NPARG,*),DD_IAD(NSPMD+1,*),EADD(*),
     .        IPART(LIPART1,*),IPARTSP(*),CEPSP(*),IXSP(KVOISPH,NUMSPH),
     .        IPM(NPROPMI,*), IGEO(NPROPGI,*),
     .        sph2sol(*), sol2sph(2,*), irst(3,*), nod2sp(*)
      INTEGER, INTENT(IN) :: PRINT_FLAG !< flag to print the element group data
      my_real
     .        PM(NPROPM,*), SPBUF(NSPBUF,NUMSPH)
      TYPE(MATPARAM_STRUCT_) ,DIMENSION(NUMMAT), INTENT(IN) :: MAT_PARAM
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NGR1, NG, MT, MLN, I, P, NEL, MODE, NB,
     .        n, igtyp,jivf,jhbe,ijk,ne1,
     .        issn,iksnod,iorder,iprt,isleep,ieos,nel_prec,iun,ig,ifail, 
     .        work(70000),ngp(nspmd+1),k,j,ii, mx, nfail, ir, ip, stat,
     .        ipartr2r, nod, jale_from_mat, jale_from_prop
      INTEGER, DIMENSION(:,:),ALLOCATABLE :: INUM
      INTEGER, DIMENSION(:),ALLOCATABLE :: INDEX     
      my_real, DIMENSION(:,:),ALLOCATABLE :: RNUM
 
      INTEGER ID
      CHARACTER(LEN=NCHARTITLE)::TITR
      DATA iun/1/
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
      my_real
     .         get_u_geo
      EXTERNAL get_u_geo
C--------------------------------------------------------------
C         BORNAGE DES GROUPES DE MVSIZ
C--------------------------------------------------------------
      CALL my_alloc(inum,13,numsph)
      CALL my_alloc(index,2*numsph)
      CALL my_alloc(rnum,nspbuf,numsph)
      nel = 0
      DO n=1,nd
        nel = nel + eadd(n+1)-eadd(n)
      ENDDO
      ngr1 = ngroup + 1
C
C phase 1 : domain decompostition
C
      idx=idx+nd*(nspmd+1)
      CALL zeroin(1,nd*(nspmd+1),dd_iad(1,nspgroup+1))
C     NSPGROUP = NSPGROUP + ND
      nft = 0
C initialisation dd_iad
      DO n=1,nd
       DO p=1,nspmd+1
         dd_iad(p,nspgroup+n) = 0
       END DO
      END DO
C
      DO n=1,nd
        nel = eadd(n+1)-eadd(n)
        DO i = 1, nel
          index(i) = i
          inum(1,i)=ipartsp(nft+i)
          inum(2,i)=kxsp(1,nft+i)
          inum(3,i)=kxsp(2,nft+i)
          inum(4,i)=kxsp(3,nft+i)
          inum(5,i)=kxsp(4,nft+i)
          inum(6,i)=kxsp(5,nft+i)
          inum(7,i)=kxsp(6,nft+i)
          inum(8,i)=kxsp(7,nft+i)
          inum(9,i)=kxsp(8,nft+i)
c retri spbuf
          DO k=1,nspbuf
           rnum(k,i)=spbuf(k,nft+i)
          END DO
        END DO
 
        DO i = 1, nel
          DO j = 1, kvoisph
            ixsps(j,i) = ixsp(j,nft+i)
          END DO
        END DO
        mode=0
        CALL my_orders( mode, work, cepsp(nft+1), index, nel , 1)
        DO i = 1, nel
          ipartsp(i+nft)=inum(1,index(i))
          kxsp(1,i+nft)=inum(2,index(i))
          kxsp(2,i+nft)=inum(3,index(i))
          kxsp(3,i+nft)=inum(4,index(i))
          kxsp(4,i+nft)=inum(5,index(i))
          kxsp(5,i+nft)=inum(6,index(i))
          kxsp(6,i+nft)=inum(7,index(i))
          kxsp(7,i+nft)=inum(8,index(i))
          kxsp(8,i+nft)=inum(9,index(i))
 
c tri SPBUF
          DO k=1,nspbuf
            spbuf(k,i+nft)=rnum(k,index(i))
          END DO
        END DO
 
        DO i = 1, nel
          DO j = 1, kvoisph
            ixsp(j,i+nft) = ixsps(j,index(i))
          END DO
        END DO
C
        IF(nsphsol/=0)THEN
          DO i=1,nel
            inum(10,i)=sph2sol(nft+i)
            IF(nft+i >= first_sphsol .AND. 
     .         nft+i < first_sphsol+nsphsol)THEN
              inum(11,i)=irst(1,nft+i-first_sphsol+1)
              inum(12,i)=irst(2,nft+i-first_sphsol+1)
              inum(13,i)=irst(3,nft+i-first_sphsol+1)
            END IF
          END DO
          DO i=1,nel
            sph2sol(nft+i) = inum(10,index(i))
C
            IF(nft+i >= first_sphsol .AND. 
     .         nft+i < first_sphsol+nsphsol)THEN
C INDEX(I) < FIRST_SPHSOL <=> internal error
              irst(1,nft+i-first_sphsol+1)=inum(11,index(i))
              irst(2,nft+i-first_sphsol+1)=inum(12,index(i))
              irst(3,nft+i-first_sphsol+1)=inum(13,index(i))
            END IF
          END DO
        END IF
C
        p = cepsp(nft+index(1))
        nb = 1
        DO i = 2, nel
          IF (cepsp(nft+index(i))/=p) THEN
            dd_iad(p+1,nspgroup+n) = nb
            nb = 1
            p = cepsp(nft+index(i))
          ELSE
            nb = nb + 1
          ENDIF
        END DO
        dd_iad(p+1,nspgroup+n) = nb
        DO p = 2, nspmd
          dd_iad(p,nspgroup+n) = dd_iad(p,nspgroup+n)
     .                         + dd_iad(p-1,nspgroup+n)
        END DO
        DO p = nspmd+1,2,-1
          dd_iad(p,nspgroup+n) = dd_iad(p-1,nspgroup+n)+1
        END DO
        dd_iad(1,nspgroup+n) = 1
C
C maj CEP
C
        DO i = 1, nel
          index(i) = cepsp(nft+index(i))          
        END DO
        DO i = 1, nel
          cepsp(nft+i) = index(i)          
        END DO
C
        nft = nft + nel
C
      END DO
C
C Rebuild SOL2SPH, SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
      IF(nsphsol/=0)THEN
        DO n=1,numels8
          sol2sph(1,n)=0
          sol2sph(2,n)=0
        END DO
        n=sph2sol(first_sphsol)
        sol2sph(1,n)=first_sphsol-1
        sol2sph(2,n)=sol2sph(1,n)+1
        DO i=first_sphsol+1,first_sphsol+nsphsol-1
          IF(sph2sol(i)==n)THEN
            sol2sph(2,n)=sol2sph(2,n)+1
          ELSE
            n=sph2sol(i)
            sol2sph(1,n)=i-1
            sol2sph(2,n)=sol2sph(1,n)+1
          END IF
        END DO          
      END IF
C ne pas oublier renumeroter selection th et surface si concerne
C
C
C phase 2 : bornage en groupe de mvsiz
C ngroup est global, iparg est global mais organise en fonction de dd
C
      jale=0
      jeul=0
      jtur=0
      jthe=0
      jivf=0
      jpor=0
C
      issn   =0
      npt    =1
      iksnod =1
      jhbe   =0
C
      DO n=1,nd
       nft = 0
       DO p = 1, nspmd
        ngp(p)=0
        nel = dd_iad(p+1,nspgroup+n)-dd_iad(p,nspgroup+n)
        IF (nel>0) THEN
         nel_prec = dd_iad(p,nspgroup+n)-dd_iad(1,nspgroup+n)
         ngp(p)=ngroup
         ng  = (nel-1)/nvsiz + 1
         DO i=1,ng
C ngroup global
          istrain=0
          ngroup=ngroup+1
          ii = eadd(n)+nft
          iprt  =ipartsp(ii)
          ipartr2r = 0
          IF (nsubdom>0) ipartr2r = tag_part(iprt)          
          mt    =ipart(1,iprt)
          mln   =nint(pm(19,abs(mt)))
          ig    =ipart(2,iprt)
          igtyp = igeo(11,ig)
          isorth= max(igeo(17,ig),min(iun,igeo(2,ig)))
          israt = ipm(3,mt)
          ieos  = ipm(4,mt)
          iorder=int(get_u_geo(5,ig))
          isleep=kxsp(2,ii)
          
          jale_from_mat = nint(pm(72,mt))
          jale_from_prop = igeo(62,ig)
          jale = max(jale_from_mat, jale_from_prop) !if inconsistent, error message was displayed in PART reader
 
          jlag=0
          IF(jale==0.AND.mln/=18)jlag=1
          jeul=0
          IF(jale==2)THEN
            jale=0
            jeul=1
          ENDIF
 
          !ALE UVAR REZONING  (81:MAT,  82:EOS)
          ! n/a with SPH
          iparg(81,ngroup) = 0
          iparg(82,ngroup) = 0
 
 
          IF (jale+jeul/=0) THEN
C           WRITE(ISTDO,*) ' ** ERROR : BAD ANALYSIS TYPE'
C           WRITE(IOUT,*)  ' ** ERROR : BAD ANALYSIS TYPE'
C           WRITE(IOUT,*)  ' ** ERROR : ALE AND EULERIAN ANALYSIS',
C     .                 ' ARE NOT COMPATIBLE WITH SPH.'
C           IERR=IERR+1
            id=igeo(1,ig)
            CALL fretitl2(titr,igeo(npropgi-ltitr+1,ig),ltitr)
            CALL ancmsg(msgid=403,
     .                  msgtype=msgerror,
     .                  anmode=aninfo_blind_1,
     .                  i1=id,
     .                  c1=titr)
          ENDIF
          jtur=nint(pm(70,mt))
          jthe=nint(pm(71,mt))
 
          isph2sol=0
          IF(nsphsol/=0)isph2sol=sph2sol(ii)
C     full geometric non linearities.
          ne1 = min( nvsiz, nel + nel_prec - nft)
          ifail = 0
          IF (mat_param(mt)%NFAIL > 0) ifail = 1
C-------------------------------------------------          
          IF(mln/=14 .AND. mln/=24 .AND. mln/=25 .AND. mln<28) THEN
            nfail = mat_param(mt)%NFAIL
            DO ijk = 1, ne1
              ii = eadd(n)+nft-1+ijk
              mx = ipart(1,ipartsp(ii))
              DO ir = 1,nfail 
                IF (mat_param(mx)%FAIL(ir)%IRUPT == 10) THEN
                  istrain=1
                  GO TO 100
                END IF
              END DO
            END DO
 100        CONTINUE
          ENDIF
C-------------------------------------------------          
          DO ijk = 1, ne1
C
C Attention, penser a sauvegarder KXSP(2) si besoin (backtrack spmd si interface)
C
            kxsp(2,eadd(n)-1+nft+ijk)=
     .                sign(ngroup,isleep)
          END DO
          IF (mt/=0) THEN
            iparg(1,ngroup)=mln
          ELSE
            iparg(1,ngroup)=igtyp
          END IF
          iparg(2,ngroup)=ne1
          iparg(3,ngroup)=eadd(n)-1 + nft
          iparg(4,ngroup)=lbufel+1
          iparg(5,ngroup)=51
          iparg(6,ngroup)=npt
          iparg(7,ngroup) =jale
          IF(isleep==-1.OR.isph2sol/=0)iparg(8,ngroup) =1
          iparg(9,ngroup) =issn
          IF(isleep>0)iparg(10,ngroup)=ne1
          iparg(11,ngroup)=jeul
          iparg(12,ngroup)=jtur
          iparg(13,ngroup)=-abs(jthe)
          iparg(14,ngroup)=jlag
          iparg(18,ngroup)=0    ! NMTV(MLN)-11
          iparg(23,ngroup)=jhbe
          iparg(24,ngroup)=jivf
          iparg(27,ngroup)=jpor
          iparg(28,ngroup)=iksnod
          iparg(32,ngroup)= p-1
          iparg(38,ngroup)=igtyp
          iparg(40,ngroup)=israt
          iparg(42,ngroup)=isorth
          iparg(43,ngroup)=ifail
          iparg(62,ngroup)=ig
          iparg(69,ngroup)=isph2sol
          
C         flag for group of duplicated elements in multidomains
          IF (nsubdom>0) iparg(71,ngroup)= ipartr2r
C         thermal expansion
          IF(ipm(218,mt) > 0 .AND. mln /= 0 .AND. mln /= 13) iparg(49,ngroup)= 1 
C
          IF(mln/=14.AND.mln/=24.AND.mln/=25.AND.mln<28)THEN
            iparg(44,ngroup)= istrain
          ELSEIF(mln>=28)THEN
            istrain=2
            iparg(44,ngroup)=istrain
          ENDIF
C
C         equation of state
          iparg(55,ngroup)= ieos
          nft = nft + ne1
         ENDDO
         ngp(p)=ngroup-ngp(p)
        ENDIF
       ENDDO
C DD_IAD => nb groupes par sous domaine
       ngp(nspmd+1)=0
       DO p = 1, nspmd
         ngp(nspmd+1)=ngp(nspmd+1)+ngp(p)
         dd_iad(p,nspgroup+n)=ngp(p)
       END DO
       dd_iad(nspmd+1,nspgroup+n)=ngp(nspmd+1)
C
      END DO
      nspgroup = nspgroup + nd
C
      ! rebuild NOD2SP after KXSP renumerotation
      nod2sp(1:numnod) = 0
      DO i = 1, numsph
        nod = kxsp(3,i)
        nod2sp(nod) = i
      END DO 
C
      IF(print_flag>6) THEN
          WRITE(iout,1000)
          WRITE(iout,1001)(n,iparg(1,n),iparg(2,n),iparg(3,n)+1,
     +                   iparg(4,n),iparg(5,n),iparg(55,n),
     +                 n=ngr1,ngroup)
      ENDIF
C
 1000 FORMAT(10x,' 3D - SPH CELL GROUPS '/
     +       10x,' -------------------- '/
     +' GROUP   CELL       CELL      FIRST    BUFFER   CELL     IEOS'/
     +'         MATERIAL   NUMBER    CELL     ADDRESS  TYPE     TYPE'/)
 1001 FORMAT(7(1x,i7,1x))
C-----------
      DEALLOCATE(inum)
      DEALLOCATE(index)
      DEALLOCATE(rnum)
C-----------
      RETURN
      END