hm_read_sphcel.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/.
!||====================================================================
!||    hm_read_sphcel         ../starter/source/elements/reader/hm_read_sphcel.F
!||--- called by ------------------------------------------------------
!||    lectur                 ../starter/source/starter/lectur.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                 ../starter/source/output/message/message.F
!||    fretitl2               ../starter/source/starter/freform.F
!||    hm_get_intv            ../starter/source/devtools/hm_reader/hm_get_intv.F
!||    hm_option_count        ../starter/source/devtools/hm_reader/hm_option_count.F
!||    hm_option_read_key     ../starter/source/devtools/hm_reader/hm_option_read_key.F
!||    hm_option_start        ../starter/source/devtools/hm_reader/hm_option_start.F
!||    soltosphx4             ../starter/source/elements/sph/soltosph.F
!||    soltosphx8             ../starter/source/elements/sph/soltosph.F
!||    udouble                ../starter/source/system/sysfus.F
!||    usr2sys                ../starter/source/system/sysfus.F
!||--- uses       -----------------------------------------------------
!||    hm_option_read_mod     ../starter/share/modules1/hm_option_read_mod.F
!||    message_mod            ../starter/share/message_module/message_mod.F
!||    reader_old_mod         ../starter/share/modules1/reader_old_mod.f90
!||    submodel_mod           ../starter/share/modules1/submodel_mod.F
!||====================================================================
      SUBROUTINE hm_read_sphcel(ITAB    ,ITABM1  ,IPART ,
     2                  IPARTSP ,IPM     ,IGEO  ,KXSP    ,IXSP    ,
     3                  NOD2SP, RESERVEP ,IXS   ,IPARTS  ,ISOLNOD ,
     4                  SPH2SOL ,SOL2SPH ,IRST  ,X       ,SOL2SPH_TYP,
     5                  LSUBMODEL,SPBUF  ,UNITAB,IPRI    )
C-----------------------------------------------
C   ROUTINE DESCRIPTION :
C   ===================
C   READ /SPHCEL ELEMENTS USING HM_READER
C-----------------------------------------------
C   DUMMY ARGUMENTS DESCRIPTION:
C   ===================
C
C     NAME            DESCRIPTION                         
C
C     IXQ             /QUAD ARRAY : CONNECTIVITY, ID, MID PID
C     ITAB            USER ID OF NODES         
C     ITABM1          REVERSE TAB ITAB
C     IPART           PART ARRAY      
C     IPARTSP          INTERNAL PART ID OF A GIVEN QUAD (INTERNAL ID)   
C     IPM             MATERIAL ARRAY (INTEGER) 
C     IGEO            PROP ARRAY (INTEGER) 
C     LSUBMODEL       SUBMODEL STRUCTURE   
C
C     KXSP(1,*)              :INUTILISE
C     KXSP(2,*)              :NG     :   +/- NO DU GROUPE
C     KXSP(3,*)              :IPRC   :   NO SYSTEME DU NOEUD ASSOCIE
C     KXSP(4,*)              :NVOIS  :   NOMBRE DE VOISINS.
C     KXSP(5,*)              :NVOIS  :   NOMBRE DE CANDIDATS RETENUS PAR LE BUCKET.
C     KXSP(6,*)              :NVOISS :   NOMBRE DE VOISINS DANS LA PARTIE SYMETRIQUE.
C     KXSP(7,*)              :NVOISS :   NOMBRE DE CANDIDATS RETENUS DANS LA PARTIE SYMETRIQUE.
C     KXSP(NISP,*)           :ID     :   ID DE LA CELLULE.
C     UNITAB          UNIT ARRAY 
C     IPRI            PRINT FLAG (in 0.out file)
C-----------------------------------------------
C     IXSP(1:KVOISPH,*)      :IVOIS  :   NOS DES VOISINS. 
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE my_alloc_mod
      USE message_mod
      USE hm_option_read_mod
      USE submodel_mod , ONLY : nsubmod, submodel_data
      USE unitab_mod
      USE names_and_titles_mod , ONLY : nchartitle, ncharkey
      USE reader_old_mod , ONLY : line, kcur, ksphopt, irec, koptad
      USE user_id_mod , ONLY : id_limit
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      "com01_c.inc"
#include      "com04_c.inc"
#include      "units_c.inc"
#include      "sphcom.inc"
#include      "scr17_c.inc"
#include      "param_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),
     .        NOD2SP(*),IPM(NPROPMI,*),IGEO(NPROPGI,*),
     .        ITAB(*),ITABM1(*),IPART(LIPART1,*),IPARTSP(*),
     .        RESERVEP(NBPARTINLET), IXS(NIXS,*), IPARTS(*), ISOLNOD(*),
     .        sph2sol(*), sol2sph(2,*), irst(3,nsphsol),sol2sph_typ(*)
      my_real x(3,*)
      TYPE(submodel_data) LSUBMODEL(NSUBMOD)
      my_real, INTENT(INOUT) :: SPBUF(NSPBUF,NUMSPH)
      TYPE (UNIT_TYPE_),INTENT(IN) ::UNITAB 
      INTEGER, INTENT(IN) :: IPRI
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,N,J,ID,IDS,K,
     .        MID,PID,IPRT,IPIDS,NSPHDIR,
     .        NSPHCEL,NCELL,IDNOD,INOD,IDMAX,KSPHRES,
     .        NBP,IT,NT,NP,NN,ITOPO,STAT,
     .        index_part,uid,iflagunit
      my_real bid,fac_m
      CHARACTER MESS*40
      CHARACTER(LEN=NCHARKEY) :: KEY
      LOGICAL IS_AVAILABLE
 
      LOGICAL :: CHECK_LAW
      INTEGER :: MID_SPH,MID_SOL
      INTEGER :: LAW_SPH,LAW_SOL
      INTEGER :: ERROR_NUM
      INTEGER :: I1,I2,I3,I4,I5
      CHARACTER(LEN=NCHARTITLE) :: C1
      CHARACTER(LEN=NCHARTITLE) :: TITR
      INTEGER :: USER_PART_SPH,USER_PART_SOL
      INTEGER :: USER_MID_SPH,USER_MID_SOL
      LOGICAL, DIMENSION(NPART) :: TAG_PART
      INTEGER, DIMENSION(NPART) :: PART_ID_SPH,PART_ID_SOL
      INTEGER, DIMENSION(:), ALLOCATABLE :: SUB_SPH
      INTEGER, DIMENSION(:), ALLOCATABLE :: TYPE
      real*8, DIMENSION(:), ALLOCATABLE :: hm_mass
      INTEGER, DIMENSION(:), ALLOCATABLE :: UID_SPH
      INTEGER, DIMENSION(:), ALLOCATABLE :: ITAG
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
      INTEGER USR2SYS
C-----------------------------------------------
      DATA MESS /'SPH CONNECTIVITIES DEFINITION           '/
C--------------------------------------------------
C      ALLOCS & INITS
c      use NUMELQ IN PLACE OF NUMELC ( NBADMESH routine is modifying NUMELC )
C--------------------------------------------------
      CALL my_alloc(itag,numnod)
      ALLOCATE (sub_sph(numsph),stat=stat)
      IF (stat /= 0) CALL ancmsg(msgid=268,anmode=aninfo,
     .                             msgtype=msgerror,
     .                             c1='SUB_SPH')
      ALLOCATE (uid_sph(numsph),stat=stat)
      IF (stat /= 0) CALL ancmsg(msgid=268,anmode=aninfo,
     .                             msgtype=msgerror,
     .                             c1='UID_SPH') 
      ALLOCATE (TYPE(numsph),STAT=stat)
      IF (stat /= 0) CALL ancmsg(msgid=268,anmode=aninfo,
     .                             msgtype=msgerror,
     .                             c1='TYPE') 
      ALLOCATE (hm_mass(numsph),stat=stat)
      IF (stat /= 0) CALL ancmsg(msgid=268,anmode=aninfo,
     .                             msgtype=msgerror,
     .                             c1='HM_MASS') 
      sub_sph(1:numsph) = 0
      uid_sph(1:numsph) = 0
      hm_mass(1:numsph) = zero
      TYPE(1:numsph) = 0
      index_part = 1
      uid = -1 
 
C--------------------------------------------------
C      READING QUADS INPUTS IN HM STRUCTURE
C-------------------------------------------------- 
      CALL cpp_sphcel_read(kxsp,nisp,ipartsp,sub_sph,TYPE,hm_mass,uid_sph)
C--------------------------------------------------
      call hm_option_count('/SPHCEL',nsphcel)
      NCELL=0
      idmax=0
C--------------------------------------------------
C      FILL OTHER STRUCTURES + CHECKS
C--------------------------------------------------
      DO i=1,nsphcel
C--------------------------------------------------
C      SUBMODEL OFFSET
C--------------------------------------------------
        IF(sub_sph(i) /= 0)THEN
          IF(uid_sph(i) == 0 .AND. lsubmodel(sub_sph(i))%UID /= 0) 
     .         uid_sph(i) = lsubmodel(sub_sph(i))%UID
        ENDIF
C--------------------------------------------------
C      UNITS
C--------------------------------------------------
        fac_m = one
        IF(uid_sph(i) /= uid )THEN
          uid = uid_sph(i)
          iflagunit = 0          
          DO j=1,unitab%NUNITS                         
            IF (unitab%UNIT_ID(j) == uid) THEN  
              fac_m = unitab%FAC_M(j)                 
              iflagunit = 1                     
            ENDIF                               
          ENDDO                      
          IF (uid/=0.AND.iflagunit==0) THEN
            CALL ancmsg(msgid=643,anmode=aninfo,msgtype=msgerror,
     .                  i1=uid,c1='/SPHCELL')
          ENDIF 
        ENDIF 
        hm_mass(i) = hm_mass(i) * fac_m
C--------------------------------------------------
C       STORAGE OF MASS PER PARTICULE
C--------------------------------------------------
        IF ((TYPE(i)==0).AND.(hm_mass(i) > zero)) THEN
          TYPE(i)=1
        endif 
        spbuf(12,i) = hm_mass(i)
        spbuf(13,i) = TYPE(i)
C--------------------------------------------------
C       INTERNAL PART ID
C--------------------------------------------------
        IF( ipart(4,index_part) /= ipartsp(i) )THEN  
          DO j=1,npart                            
            IF(ipart(4,j)== ipartsp(i) ) index_part = j 
          ENDDO  
        ENDIF
        IF(ipart(4,index_part) /= ipartsp(i)) THEN
          CALL ancmsg(msgid=402,
     .                msgtype=msgerror,
     .                anmode=aninfo_blind_1,
     .                c1="SPHCEL",
     .                i1=ipartsp(i),
     .                i2=ipartsp(i),
     .                prmod=msg_cumu)
        ENDIF 
        idnod = kxsp(3,i)
        inod=usr2sys(idnod,itabm1,mess,id)
        kxsp(3,i)=inod
        ncell=ncell+1
        ipartsp(ncell)=index_part
        nod2sp(inod) =ncell
C         meme identifiant que le noeud.
        kxsp(nisp,ncell)=idnod
        idmax=max(idmax,idnod)
C--------------------------------------------------
        IF (kxsp(nisp,i)>id_limit%GLOBAL)THEN
            CALL ancmsg(msgid=509,anmode=aninfo,msgtype=msgerror,
     .                  i1=kxsp(nisp,i),c1=line,c2='/SPHCEL')
        ENDIF
      ENDDO
C-------------------------------------
       IF(ALLOCATED(sub_sph)) DEALLOCATE(sub_sph)
       IF(ALLOCATED(uid_sph)) DEALLOCATE(uid_sph)
       IF(ALLOCATED(type)) DEALLOCATE(type)
       IF(ALLOCATED(hm_mass)) DEALLOCATE(hm_mass)
C-------------------------------------
       first_sphres=ncell+1
       IF(nsphres/=0)THEN
        kcur   =ksphopt
        irec   =koptad(kcur)-1
        inod   =isphres
        nbp = 1
 
        CALL hm_option_start('/SPH/RESERVE')
        DO n=1,nbpartinlet 
          CALL hm_option_read_key(lsubmodel,
     .                         option_id = id,
     .                         keyword2 = key)
          ids=0
          DO j=1,npart
            IF(ipart(4,j)==id) THEN
             IF(igeo(11,ipart(2,j))/=34)THEN
               CALL ancmsg(msgid=1068,
     .                     msgtype=msgerror,
     .                     anmode=aninfo,
     .                     i1=id,
     .                     i2=id)
             ELSE
               ids=j
             END IF
             GOTO 175
            ENDIF
          ENDDO
          CALL ancmsg(msgid=441,
     .                msgtype=msgerror,
     .                anmode=aninfo,
     .                i1=id,
     .                i2=id)
175       CONTINUE
          CALL hm_get_intv('Np',ksphres,is_available,lsubmodel)
C if SPH part not found reserve is emptied to prevent additional errors 
          IF (ids==0) THEN
            nsphres = nsphres - ksphres*nspmd
            numsph = numsph - ksphres*nspmd
            ksphres = 0
          ENDIF 
c store nb of reserve for this part  
          reservep(nbp)=ksphres
          nbp=nbp+1
c KSPHRES by proc  
          ksphres = ksphres*nspmd
          DO j=1,ksphres
            ncell=ncell+1
            ipartsp(ncell)=ids
            inod =inod+1
            kxsp(3,ncell) =inod
            nod2sp(inod)  =ncell
            kxsp(2,ncell)=-1
            idmax=idmax+1
            kxsp(nisp,ncell)=idmax
          ENDDO
        ENDDO
       ENDIF
C-------------------------------------
       first_sphsol=ncell+1
       IF(nsphsol/=0)THEN
C
        kcontact=1
C
        inod   =firstnod_sphsol-1
        DO n=1,numels8
          sol2sph(1,n)=0
          sol2sph(2,n)=0
          ipids  =ipart(2,iparts(n))
          nsphdir=igeo(37,ipids)
          ids    =igeo(38,ipids)
          IF(nsphdir/=0)THEN
            IF(isolnod(n)==8)THEN
              DO j=1,8
                itag(ixs(1+j,n))=0
              END DO
              nn=0
              DO j=1,8
                IF(itag(ixs(1+j,n))==0)THEN
                   nn=nn+1
                   itag(ixs(1+j,n))=1
                 END IF
              END DO
              IF(nn==4)THEN
                itopo=4
                np=0
                nt=0
                DO it=1,nsphdir
                  nt=nt+it
                  np=np+nt
                END DO
              ELSE
                itopo=8
                np=nsphdir*nsphdir*nsphdir
              END IF
            ELSEIF(isolnod(n)==6)THEN
              itopo=6
              np=0
              nt=0
              DO it=1,nsphdir
                nt=nt+it
              END DO
              np=np+nsphdir*nt
            ELSEIF(isolnod(n)==4)THEN
              itopo=4
              np=0
              nt=0
              DO it=1,nsphdir
                nt=nt+it
                np=np+nt
              END DO
            END IF
C---
C           SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
            sol2sph(1,n)=ncell
            sol2sph(2,n)=ncell+np
            sol2sph_typ(n)=itopo
            DO i=1,np
              sph2sol(ncell+i)=n
            END DO
C---
            IF(itopo==4)THEN
C-- Tetra
              CALL soltosphx4(nsphdir,ncell   ,inod    ,ids    ,idmax   ,
     .                       x      ,ixs(1,n),kxsp    ,ipartsp,nod2sp  ,
     .                       irst   )
            ELSEIF (itopo==8) THEN
C-- Hexa + degenerated penta6
              CALL soltosphx8(nsphdir,ncell   ,inod    ,ids    ,idmax   ,
     .                       x      ,ixs(1,n),kxsp    ,ipartsp,nod2sp  ,
     .                       irst   )
            ENDIF
          ENDIF
        ENDDO
       ENDIF
! ------------------------------------
!   check the material consistency between
!   the solide and the sph (SOL2PSH)
       check_law = .false.
       error_num = 0
       tag_part(1:npart) = .false.
       IF(nsphsol/=0) THEN
            DO i =1,numsph
                n = sph2sol(i) 
                IF(n/=0) THEN
                    mid_sph = ipart(1,ipartsp(i))
                    mid_sol = ipart(1,iparts(n))            
                    law_sph = ipm(2,mid_sph)
                    law_sol = ipm(2,mid_sol)                                 
                    IF(law_sph/=law_sol) THEN
                        check_law = .true.
                        IF( .NOT.tag_part(ipartsp(i)) ) THEN
                            error_num = error_num + 1
                            tag_part(ipartsp(i)) = .true.
                            part_id_sph(error_num) = ipartsp(i)
                            part_id_sol(error_num) = iparts(n)
                        ENDIF
                    ENDIF
                ENDIF
            ENDDO
       ENDIF
    
       IF(check_law) THEN
            DO i=1,error_num
                titr(1:nchartitle) =''
                CALL fretitl2(titr,ipart(lipart1-ltitr+1,part_id_sph(i)),ltitr-1)
                user_part_sph = ipart(4,part_id_sph(i))
                user_part_sol = ipart(4,part_id_sol(i))
                user_mid_sph = ipart(5,part_id_sph(i))
                user_mid_sol = ipart(5,part_id_sol(i))
                CALL ancmsg(msgid=1911,
     .                     msgtype=msgerror,
     .                     anmode=aninfo,
     .                     i1=user_part_sph,c1=titr(1:len_trim(titr)),
     .                     i2=user_mid_sph,i3=user_part_sph,
     .                     i4=user_mid_sol,i5=user_part_sol )
            ENDDO
       ENDIF
C-------------------------------------
C Recherche des ID doubles
C-------------------------------------
      CALL udouble(kxsp(nisp,1),nisp,numsph,mess,0,bid)
C-------------------------------------
C Print
C-------------------------------------
      i1=1
      i2=min0(50,numsph)
C
      IF(ipri>=5) THEN
   90   WRITE (iout,300)
        DO 100 i=i1,i2
         iprt=ipartsp(i)
         mid =ipm(1,ipart(1,iprt))
         pid =igeo(1,ipart(2,iprt))
         WRITE (iout,'(6(I10,1X))') i,kxsp(nisp,i),mid,pid,
     .                             kxsp(3,i),itab(kxsp(3,i))
  100   CONTINUE
        IF(i2==numsph)GOTO 200
        i1=i1+50
        i2=min0(i2+50,numsph)
        GOTO 90
C
 200    CONTINUE
        WRITE (iout,'(A)') 'END OF CELL TRACEBACK'
      ENDIF
 
       DEALLOCATE(itag)
C
 300  FORMAT(/' SPH CELLS             '/
     +        ' ----------------------'/
     + '   LOC-CEL    GLO-CEL      MATER ',
     + '      GEOM    LOC-NOD    GLO-NOD ')
      RETURN
C-------------------------------------
      END