srho3.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/.
!||====================================================================
!||    srho3            ../engine/source/elements/solid/solide/srho3.F
!||--- called by ------------------------------------------------------
!||    ig3duforc3       ../engine/source/elements/ige3d/ig3duforc3.F
!||    s10forc3         ../engine/source/elements/solid/solide10/s10forc3.F
!||    s16forc3         ../engine/source/elements/thickshell/solide16/s16forc3.F
!||    s20forc3         ../engine/source/elements/solid/solide20/s20forc3.F
!||    s4forc3          ../engine/source/elements/solid/solide4/s4forc3.F
!||    s6cforc3         ../engine/source/elements/thickshell/solide6c/s6cforc3.F
!||    s8cforc3         ../engine/source/elements/thickshell/solide8c/s8cforc3.F
!||    s8eforc3         ../engine/source/elements/solid/solide8e/s8eforc3.F
!||    s8sforc3         ../engine/source/elements/solid/solide8s/s8sforc3.F
!||    s8zforc3         ../engine/source/elements/solid/solide8z/s8zforc3.F
!||    scforc3          ../engine/source/elements/thickshell/solidec/scforc3.f
!||    sforc3           ../engine/source/elements/solid/solide/sforc3.F
!||    szforc3          ../engine/source/elements/solid/solidez/szforc3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg           ../engine/source/output/message/message.F
!||    arret            ../engine/source/system/arret.F
!||--- uses       -----------------------------------------------------
!||    ale_mod          ../common_source/modules/ale/ale_mod.F
!||    i22bufbric_mod   ../common_source/modules/interfaces/cut-cell-search_mod.F
!||    message_mod      ../engine/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE srho3(
     1   PM,        VOLO,      RHON,      EINT,
     2   DIVDE,     FLUX,      FLU1,      VOLN,
     3   DVOL,      NGL,       MAT,       OFF,
     4   IS_MAT_BCS,TAG22,     VOLDP,     VOL0DP,
     5   AMU,       OFFG,      NEL,       MTN,
     6   JALE,      ISMSTR,    JEUL,      JLAG)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
      USE i22bufbric_mod
      USE ale_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      "param_c.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com08_c.inc"
#include      "scr06_c.inc"
#include      "impl1_c.inc"
#include      "com01_c.inc"
#include      "com04_c.inc"
#include      "inter22.inc"
#include      "scr05_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: NEL
      INTEGER, INTENT(IN) :: MTN
      INTEGER, INTENT(IN) :: JALE
      INTEGER, INTENT(IN) :: ISMSTR
      INTEGER, INTENT(IN) :: JEUL
      INTEGER, INTENT(IN) :: JLAG
      INTEGER NGL(*), MAT(*), IS_MAT_BCS, IB,NIN,MCELL
      
      my_real pm(npropm,nummat),volo(*), rhon(*),eint(*),flux(6,*), flu1(*),
     .        voln(*), dvol(*),divde(*),off(*),tag22(*),amu(*) ,offg(*)
      DOUBLE PRECISION VOLDP(*),VOL0DP(*),DVDP
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      my_real :: rho0,dvv, e0,vavg,rv,rvp,rhon_old(mvsiz),ddvol,rhoref
      INTEGER :: I, J,COUNT,LIST(MVSIZ),II, MX
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      rhon_old(1:nel)=rhon(1:nel)
      rho0 = zero
      IF(ale%GLOBAL%INCOMP==1 .AND. jeul+jale==1)THEN
        mx        = mat(1)
        rho0      = pm(1,mx)
        e0        = pm(23,mx)
        DO i=1,nel
          dvv     = divde(i)
          rhon(i) = rhon(i)-dvv*rho0
          vavg    = half*(voln(i)+volo(i))
          dvol(i) = vavg*dvv
          eint(i) = eint(i)*volo(i)-e0*dvv*vavg
        ENDDO
      ELSE
        IF(jlag/=0)THEN
          mx      = mat(1)
          rho0    = pm(1,mx)
C--- due to /INIBRI/EREF        
          IF (tt==zero) THEN
           IF (ismstr==11) THEN
            volo(1:nel)=voln(1:nel)
           ELSEIF(ismstr==1) THEN
            DO i=1,nel
              IF(offg(i)>one) volo(i)=voln(i)
            ENDDO
           END IF
          END IF
          IF (impl_s>0.AND.iline>0) THEN
            DO i=1,nel
              rhon(i) = rho0
              eint(i) = eint(i)*voln(i)
            ENDDO
          ELSE
            IF (mtn /= 115) THEN 
              DO i=1,nel
                IF(offg(i)==zero.AND.voln(i)==one) voln(i)=volo(i)
                dvol(i) = voln(i)-(rho0/rhon(i))*volo(i)
                rhon(i) = rho0*(volo(i)/voln(i))
                eint(i) = eint(i)*volo(i)
              ENDDO
            ELSE
              DO i=1,nel
                IF(offg(i)==zero.AND.voln(i)==one) voln(i)=volo(i)
                dvol(i) = voln(i)-(rhon(i+nel)/rhon(i))*volo(i)
                rhon(i) = rhon(i+nel)*(volo(i)/voln(i))
                eint(i) = eint(i)*volo(i)
              ENDDO            
            ENDIF
          ENDIF
        ELSEIF(jale/=0)THEN                                                                                              
          DO i=1,nel                                                                                                     
           rhon(i)    = rhon(i)/voln(i)                                                                                    
           dvol(i)    = voln(i)-volo(i)+half*dt1*(flu1(i)+flux(1,i)+flux(2,i)+flux(3,i)+flux(4,i)+flux(5,i)+flux(6,i))   
           volo(i)    = voln(i)                                                                                            
         ENDDO                                                                                                             
        ELSEIF(jeul/=0)THEN                                                                                              
         DO i=1,nel                                                                                                      
          rhon(i)     = rhon(i)/voln(i)                                                                                    
          dvol(i)     = half*dt1*(flu1(i)+flux(1,i)+flux(2,i)+flux(3,i)+flux(4,i)+flux(5,i)+flux(6,i))                   
         ENDDO                                                                                                             
        endif!(JLAG/=0) 
        
        !---interface22---!
        IF(int22/=0)THEN
          IF(jeul+jale/=0)THEN      
           nin        = 1
           DO i=1,nel         
             ib       = nint(tag22(i))
             IF(ib==0)cycle
             mcell    = brick_list(nin,ib)%mainID
             ddvol    = brick_list(nin,ib)%POLY(mcell)%DDVOL
             dvol(i)  = dt1 * ddvol
             IF(jeul/=0)THEN
               rhon(i) = rhon(i) * voln(i) / brick_list(nin,ib)%vnew_scell
             ENDIF
             voln(i)  = brick_list(nin,ib)%vnew_scell   
             volo(i)  = brick_list(nin,ib)%vold_scell  
             dvol(i)  = dvol(i) + voln(i)-volo(i)           !USE ALE FORMULATION FOR POLYHEDRA EVEN  WITH EULERIAN MATERIAL             
             volo(i)  = voln(i) 
                           
             brick_list(nin,ib)%vold_scell  = voln(i) !for convection during next cycle : aconve() 
C             write(*,FMT='(A,I10,A,F30.16)') "brick id =", NGL(I), "   mass= ", RHON(I)*VOLN(I)
           enddo!next I
          endif!(JEUL+JALE/=0)        
        endif!INT22 
        !-----------------!        
                                                                                                                     
      endif!IF(ALE%GLOBAL%INCOMP==1 .AND. JEUL+JALE==1)THEN
 
      IF(jale+jeul/=0)THEN
       count=0
       DO i=1,nel
         IF(is_mat_bcs== 1)cycle
         IF(rhon(i)> zero)cycle
         IF(off(i)== zero )cycle
           count       = count + 1
           list(count) = i
       ENDDO
       
       DO ii = 1,count
         i = list(ii)
         CALL ancmsg(msgid=167,anmode=aninfo,i1=ngl(i),r1=rhon(i))
         CALL arret(2)
       ENDDO
      ENDIF
C      
      IF (ismdisp>0.OR.ismstr==11) THEN
C      change DXX,DYY,DZZ by DIVDE(I)=dt1*(DXX(I)+DYY(I)+DZZ(I)) as input
C---------RHON(I) = RHO0*(VOLO(I)/VOLN(I)) just calculated above
C---------VOLN(I)=VOLO(I) for small strain excepting initial stat case
       DO i=1,nel
         dvdp   = divde(i)
         rhon(i) = rhon_old(i) - rhon(i)*dvdp
         rhon(i) = max(rhon(i),em30)
         dvol(i)=voln(i)*dvdp
       ENDDO
      ELSEIF ((ismstr<=4.OR.ismstr==12).AND.jlag>0) THEN
       DO i=1,nel
        IF(offg(i)>one) THEN
         dvdp  = divde(i)
         rhoref = rhon(i)
         rhon(i) = rhon_old(i) - rhoref*dvdp
         rhon(i) = max(rhon(i),em30)
         dvol(i)=voln(i)*dvdp
         IF (ismstr==12) amu(i) =rhon(i)/rhoref - one
        END IF
       ENDDO
      ENDIF
C      
      IF((ale%GLOBAL%INCOMP/=1 .OR. (jeul+jale)/=1).AND.jlag/=0.AND.n2d==0
     .     .AND.impl_s==0.AND.ismstr/=1.AND.ismstr/=3.AND.ismstr/=11)THEN
       IF(iresp==1)THEN
C------VOLDP doesn't change after switching to small strain,modifying VOL0DP to      
        DO i=1,nel
         IF(offg(i)>one) THEN
          dvdp = divde(i)*(volo(i)/voln(i))
          vol0dp(i)=vol0dp(i)-dvdp*voldp(i)
         ELSEIF(offg(i)==zero) THEN
          voldp(i)=volo(i)
         ELSE
          dvdp = voldp(i)-(rho0/rhon_old(i))*vol0dp(i)
          dvol(i) = dvdp
          rhon(i) = rho0*(vol0dp(i)/voldp(i))
         END IF
        ENDDO
        amu(1:nel) = vol0dp(1:nel)/voldp(1:nel) - one
       END IF
      ENDIF
C-----AMU compute for DP in MMAIN: AMU() = RHON()/RHO0 - 1, because RHON*VOLDP=RHO0*VOL0DP     
      
C-----------------------------------------------
      RETURN
      END SUBROUTINE srho3