sfint3.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/.
!||====================================================================
!||    sfint3        ../engine/source/elements/solid/solide/sfint3.F
!||--- called by ------------------------------------------------------
!||    sforc3        ../engine/source/elements/solid/solide/sforc3.f
!||    szforc3       ../engine/source/elements/solid/solidez/szforc3.F
!||--- uses       -----------------------------------------------------
!||    ale_mod       ../common_source/modules/ale/ale_mod.F
!||    alefvm_mod    ../common_source/modules/ale/alefvm_mod.F
!||    element_mod   ../common_source/modules/elements/element_mod.F90
!||====================================================================
      SUBROUTINE sfint3(
     1   SIG,     PX1,     PX2,     PX3,
     2   PX4,     PY1,     PY2,     PY3,
     3   PY4,     PZ1,     PZ2,     PZ3,
     4   PZ4,     PX5,     PX6,     PX7,
     5   PX8,     PY5,     PY6,     PY7,
     6   PY8,     PZ5,     PZ6,     PZ7,
     7   PZ8,     F11,     F21,     F31,
     8   F12,     F22,     F32,     F13,
     9   F23,     F33,     F14,     F24,
     A   F34,     F15,     F25,     F35,
     B   F16,     F26,     F36,     F17,
     C   F27,     F37,     F18,     F28,
     D   F38,     VOL,     QVIS,    N1X,
     E   N2X,     N3X,     N4X,     N5X,
     F   N6X,     N1Y,     N2Y,     N3Y,
     G   N4Y,     N5Y,     N6Y,     N1Z,
     H   N2Z,     N3Z,     N4Z,     N5Z,
     I   N6Z,     DFE,     Pbak,    IXS,
     J   NEL,     NFT,     JALE,    JEUL,
     K   SVIS)
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
      !-------------------------------------------------!
      ! ISFINT=1 : volume integration for SIG_total     !
      ! ISFINT=2 : volume integration for SIG_dev only  !
      ! ISFINT=3 : surface integration for SIG_total    !
      !-------------------------------------------------!  
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE alefvm_mod , only:alefvm_param
      USE ale_mod
      use element_mod , only : nixs
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      "inter22.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: NFT
      INTEGER, INTENT(IN) :: JALE
      INTEGER, INTENT(IN) :: JEUL
      INTEGER IXS(NIXS,*),NEL
      my_real
     .   SIG(NEL,6),
     .   PX1(*), PX2(*), PX3(*), PX4(*),  
     .   PY1(*), PY2(*), PY3(*), PY4(*),  
     .   PZ1(*), PZ2(*), PZ3(*), PZ4(*),  
     .   PX5(*), PX6(*), PX7(*), PX8(*),  
     .   PY5(*), PY6(*), PY7(*), PY8(*),  
     .   PZ5(*), PZ6(*), PZ7(*), PZ8(*),  
     .   F11(*),F21(*),F31(*),F12(*),F22(*),F32(*),
     .   F13(*),F23(*),F33(*),F14(*),F24(*),F34(*),
     .   F15(*),F25(*),F35(*),F16(*),F26(*),F36(*),
     .   F17(*),F27(*),F37(*),F18(*),F28(*),F38(*),
     .   VOL(*),QVIS(*)
      my_real
     .   N1X(*), N2X(*), N3X(*), N4X(*), N5X(*), N6X(*),
     .   N1Y(*), N2Y(*), N3Y(*), N4Y(*), N5Y(*), N6Y(*),
     .   n1z(*), n2z(*), n3z(*), n4z(*), n5z(*), n6z(*),
     .   dfe(mvsiz,3),pbak(mvsiz),
     .   tx,ty,tz
      my_real, DIMENSION(MVSIZ,6), INTENT(INOUT) :: svis
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real S1(MVSIZ), S2(MVSIZ), S3(MVSIZ),S4(MVSIZ), S5(MVSIZ), S6(MVSIZ),FINT,QVIS_LOC,VOL_LOC,P(MVSIZ)
C-----------------------------------------------
 
      !-------------------------------------------------!
      ! ISFINT=3 : surface integration for SIG_total    !
      !-------------------------------------------------!  
      IF(ALE%GLOBAL%ISFINT==3 .AND. jale+jeul/=0)THEN
        !------------------------------!
        !    TOTAL STRESS TENSOR       ! 
        !------------------------------!  
        DO i=1,nel                         
          qvis_loc = qvis(i)                 
          s1(i)=sig(i,1)+svis(i,1)-qvis_loc  
          s2(i)=sig(i,2)+svis(i,2)-qvis_loc  
          s3(i)=sig(i,3)+svis(i,3)-qvis_loc  
          s4(i)=sig(i,4)+svis(i,4)           
          s5(i)=sig(i,5)+svis(i,5)           
          s6(i)=sig(i,6)+svis(i,6)           
        ENDDO  
        !------------------------------!
        !          PRESSURE            ! 
        !------------------------------!  
        DO i=1,nel                         
          p(i)=third*(s1(i)+s2(i)+s3(i))    
          pbak(i)=one_over_8*p(i)                       
        ENDDO          
        !------------------------------!
        !          DIV(SIGMA)          ! 
        !------------------------------!  
        ! Integral ( div(sigma) ) dV = Integral (sigma:n) dS  
        DO i=1,nel 
          ! node-1 : faces 1,4,6
          tx=n1x(i)+n4x(i)+n6x(i)
          ty=n1y(i)+n4y(i)+n6y(i)
          tz=n1z(i)+n4z(i)+n6z(i)                    
          f11(i)=f11(i)-one_over_8*(s1(i)*tx + s4(i)*ty + s6(i)*tz)
             !print *, "ID,FINT(11)=", IXS(11,I+NFT), -ONE_OVER_8*(S1(I)*TX + S4(I)*TY + S6(I)*TZ)
          f21(i)=f21(i)-one_over_8*(s4(i)*tx + s2(i)*ty + s5(i)*tz)
          f31(i)=f31(i)-one_over_8*(s6(i)*tx + s5(i)*ty + s3(i)*tz)                    
          ! node-2 : faces 1,4,5
          tx=n1x(i)+n4x(i)+n5x(i)
          ty=n1y(i)+n4y(i)+n5y(i)
          tz=n1z(i)+n4z(i)+n5z(i)                    
          f12(i)=f12(i)-one_over_8*(s1(i)*tx + s4(i)*ty + s6(i)*tz)
          f22(i)=f22(i)-one_over_8*(s4(i)*tx + s2(i)*ty + s5(i)*tz)
          f32(i)=f32(i)-one_over_8*(s6(i)*tx + s5(i)*ty + s3(i)*tz)          
          ! node-3 : faces 1,2,5
          tx=n1x(i)+n2x(i)+n5x(i)
          ty=n1y(i)+n2y(i)+n5y(i)
          tz=n1z(i)+n2z(i)+n5z(i)                    
          f13(i)=f13(i)-one_over_8*(s1(i)*tx + s4(i)*ty + s6(i)*tz)
          f23(i)=f23(i)-one_over_8*(s4(i)*tx + s2(i)*ty + s5(i)*tz)
          f33(i)=f33(i)-one_over_8*(s6(i)*tx + s5(i)*ty + s3(i)*tz)         
          ! node-4 : faces 1,2,6
          tx=n1x(i)+n2x(i)+n6x(i)
          ty=n1y(i)+n2y(i)+n6y(i)
          tz=n1z(i)+n2z(i)+n6z(i)                    
          f14(i)=f14(i)-one_over_8*(s1(i)*tx + s4(i)*ty + s6(i)*tz)
          f24(i)=f24(i)-one_over_8*(s4(i)*tx + s2(i)*ty + s5(i)*tz)
          f34(i)=f34(i)-one_over_8*(s6(i)*tx + s5(i)*ty + s3(i)*tz)   
          ! node-5 : faces 3,4,6
          tx=n3x(i)+n4x(i)+n6x(i)
          ty=n3y(i)+n4y(i)+n6y(i)
          tz=n3z(i)+n4z(i)+n6z(i)                    
          f15(i)=f15(i)-one_over_8*(s1(i)*tx + s4(i)*ty + s6(i)*tz)
          f25(i)=f25(i)-one_over_8*(s4(i)*tx + s2(i)*ty + s5(i)*tz)
          f35(i)=f35(i)-one_over_8*(s6(i)*tx + s5(i)*ty + s3(i)*tz) 
          ! node-6 : faces 3,4,5
          tx=n3x(i)+n4x(i)+n5x(i)
          ty=n3y(i)+n4y(i)+n5y(i)
          tz=n3z(i)+n4z(i)+n5z(i)                    
          f16(i)=f16(i)-one_over_8*(s1(i)*tx + s4(i)*ty + s6(i)*tz)
          f26(i)=f26(i)-one_over_8*(s4(i)*tx + s2(i)*ty + s5(i)*tz)
          f36(i)=f36(i)-one_over_8*(s6(i)*tx + s5(i)*ty + s3(i)*tz)                           
          ! node-7 : faces 2,3,5
          tx=n2x(i)+n3x(i)+n5x(i)
          ty=n2y(i)+n3y(i)+n5y(i)
          tz=n2z(i)+n3z(i)+n5z(i)                    
          f17(i)=f17(i)-one_over_8*(s1(i)*tx + s4(i)*ty + s6(i)*tz)
          f27(i)=f27(i)-one_over_8*(s4(i)*tx + s2(i)*ty + s5(i)*tz)
          f37(i)=f37(i)-one_over_8*(s6(i)*tx + s5(i)*ty + s3(i)*tz)     
          ! node-8 : faces 2,3,6
          tx=n2x(i)+n3x(i)+n6x(i)
          ty=n2y(i)+n3y(i)+n6y(i)
          tz=n2z(i)+n3z(i)+n6z(i)                    
          f18(i)=f18(i)-one_over_8*(s1(i)*tx + s4(i)*ty + s6(i)*tz)
          f28(i)=f28(i)-one_over_8*(s4(i)*tx + s2(i)*ty + s5(i)*tz)
          f38(i)=f38(i)-one_over_8*(s6(i)*tx + s5(i)*ty + s3(i)*tz)                              
        ENDDO  
        RETURN                          
      ENDIF
 
      !-------------------------------------------------!
      ! ISFINT=1 : volume integration for SIG_total     !
      ! ISFINT=2 : volume integration for SIG_dev only  !
      !-------------------------------------------------!       
      IF((ale%GLOBAL%ICAA==1.OR.ale%GLOBAL%ISFINT==2) .AND. jale+jeul/=0)THEN
        DO i=1,nel
         qvis_loc = qvis(i)
         vol_loc = vol(i) 
         s1(i)=sig(i,1)+svis(i,1)-qvis_loc
         s2(i)=sig(i,2)+svis(i,2)-qvis_loc
         s3(i)=sig(i,3)+svis(i,3)-qvis_loc
         p(i)=(s1(i)+s2(i)+s3(i))/three
         s1(i)=(s1(i)-p(i))*vol_loc
         s2(i)=(s2(i)-p(i))*vol_loc
         s3(i)=(s3(i)-p(i))*vol_loc
         s4(i)=(sig(i,4)+svis(i,4))*vol_loc
         s5(i)=(sig(i,5)+svis(i,5))*vol_loc
         s6(i)=(sig(i,6)+svis(i,6))*vol_loc
         p(i)=p(i)*0.125 
        ENDDO
        DO i=1,nel
          pbak(i) = p(i)
        ENDDO
        !--------------------------------------------------------------!
        ! fint from spherical tensor : FVM                             !
        !--------------------------------------------------------------!
        DO i=1,nel
          f11(i)=f11(i)-p(i)*(n1x(i)              +n4x(i)       +n6x(i))
          f21(i)=f21(i)-p(i)*(n1y(i)              +n4y(i)       +n6y(i))
          f31(i)=f31(i)-p(i)*(n1z(i)              +n4z(i)       +n6z(i))
          f12(i)=f12(i)-p(i)*(n1x(i)              +n4x(i)+n5x(i)       )
          f22(i)=f22(i)-p(i)*(n1y(i)              +n4y(i)+n5y(i)       )
          f32(i)=f32(i)-p(i)*(n1z(i)              +n4z(i)+n5z(i)       )
          f13(i)=f13(i)-p(i)*(n1x(i)+n2x(i)              +n5x(i)       )
          f23(i)=f23(i)-p(i)*(n1y(i)+n2y(i)              +n5y(i)       )
          f33(i)=f33(i)-p(i)*(n1z(i)+n2z(i)              +n5z(i)       )
          f14(i)=f14(i)-p(i)*(n1x(i)+n2x(i)                     +n6x(i))
          f24(i)=f24(i)-p(i)*(n1y(i)+n2y(i)                     +n6y(i))
          f34(i)=f34(i)-p(i)*(n1z(i)+n2z(i)                     +n6z(i))
          f15(i)=f15(i)-p(i)*(             +n3x(i)+n4x(i)       +n6x(i))
          f25(i)=f25(i)-p(i)*(             +n3y(i)+n4y(i)       +n6y(i))
          f35(i)=f35(i)-p(i)*(             +n3z(i)+n4z(i)       +n6z(i))
          f16(i)=f16(i)-p(i)*(             +n3x(i)+n4x(i)+n5x(i)       )
          f26(i)=f26(i)-p(i)*(             +n3y(i)+n4y(i)+n5y(i)       )
          f36(i)=f36(i)-p(i)*(             +n3z(i)+n4z(i)+n5z(i)       )
          f17(i)=f17(i)-p(i)*(      +n2x(i)+n3x(i)       +n5x(i)       )
          f27(i)=f27(i)-p(i)*(      +n2y(i)+n3y(i)       +n5y(i)       )
          f37(i)=f37(i)-p(i)*(      +n2z(i)+n3z(i)       +n5z(i)       )
          f18(i)=f18(i)-p(i)*(      +n2x(i)+n3x(i)              +n6x(i))
          f28(i)=f28(i)-p(i)*(      +n2y(i)+n3y(i)              +n6y(i))
          f38(i)=f38(i)-p(i)*(      +n2z(i)+n3z(i)              +n6z(i)) 
        ENDDO
      ELSE !IF(ICAA == 0 .OR. JALE+JEUL==0 )
        DO i=1,nel
         qvis_loc = qvis(i)
         vol_loc = vol(i) 
         s1(i)=(sig(i,1)+svis(i,1)-qvis_loc)*vol_loc
         s2(i)=(sig(i,2)+svis(i,2)-qvis_loc)*vol_loc
         s3(i)=(sig(i,3)+svis(i,3)-qvis_loc)*vol_loc
         s4(i)=(sig(i,4)+svis(i,4))*vol_loc
         s5(i)=(sig(i,5)+svis(i,5))*vol_loc
         s6(i)=(sig(i,6)+svis(i,6))*vol_loc 
        ENDDO
        IF(alefvm_param%IEnabled/=0)THEN
          DO i=1,nel
            pbak(i) =0.125 * (s1(i)+s2(i)+s3(i))/three/vol(i)
          ENDDO
        ENDIF
      ENDIF      
      IF(int22==0)THEN !no yet deviatoric stress contribution with inter22 to simplify (priority : ditching)
        !--------------------------------------------------------------!
        ! FEM internal forces                                          !
        !  /CAA or /INTER/TYPE22 : deviatoric stress tensor            !
        !  otherwise             : total stress tensor                 !      
        !--------------------------------------------------------------!
        IF(jeul==0 .OR. (jeul==1.AND.integ8==0))THEN
         DO i=1,nel
          !shape functions hypothesis : phi1=-phi7, phi2=-phi8, phi3=-phi5, phi4=-phi6
          !---nodes 1-7 ---!
          fint=s1(i)*px1(i)+s4(i)*py1(i)+s6(i)*pz1(i)
          f11(i)=f11(i)-fint
          f17(i)=f17(i)+fint
          fint=s2(i)*py1(i)+s4(i)*px1(i)+s5(i)*pz1(i)
          f21(i)=f21(i)-fint
          f27(i)=f27(i)+fint
          fint=s3(i)*pz1(i)+s6(i)*px1(i)+s5(i)*py1(i)
          f31(i)=f31(i)-fint
          f37(i)=f37(i)+fint
          !---nodes 2-8 ---!          
          fint=s1(i)*px2(i)+s4(i)*py2(i)+s6(i)*pz2(i)
          f12(i)=f12(i)-fint
          f18(i)=f18(i)+fint
          fint=s2(i)*py2(i)+s4(i)*px2(i)+s5(i)*pz2(i)
          f22(i)=f22(i)-fint
          f28(i)=f28(i)+fint
          fint=s3(i)*pz2(i)+s6(i)*px2(i)+s5(i)*py2(i)
          f32(i)=f32(i)-fint
          f38(i)=f38(i)+fint
          !---nodes 3-5 ---!
          fint=s1(i)*px3(i)+s4(i)*py3(i)+s6(i)*pz3(i)
          f13(i)=f13(i)-fint
          f15(i)=f15(i)+fint
          fint=s2(i)*py3(i)+s4(i)*px3(i)+s5(i)*pz3(i)
          f23(i)=f23(i)-fint
          f25(i)=f25(i)+fint
          fint=s3(i)*pz3(i)+s6(i)*px3(i)+s5(i)*py3(i)
          f33(i)=f33(i)-fint
          f35(i)=f35(i)+fint
          !---nodes 4-6 ---!
          fint=s1(i)*px4(i)+s4(i)*py4(i)+s6(i)*pz4(i)
          f14(i)=f14(i)-fint
          f16(i)=f16(i)+fint
          fint=s2(i)*py4(i)+s4(i)*px4(i)+s5(i)*pz4(i)
          f24(i)=f24(i)-fint
          f26(i)=f26(i)+fint
          fint=s3(i)*pz4(i)+s6(i)*px4(i)+s5(i)*py4(i)
          f34(i)=f34(i)-fint
          f36(i)=f36(i)+fint
         ENDDO
        ELSE
        !--------------------------------------------------------------!
        ! FEM internal forces  + INTEG8                                !
        !  /CAA or /INTER/TYPE22 : deviatoric stress tensor            !
        !  otherwise             : total stress tensor                 !   
        !--------------------------------------------------------------!
        !shape functions are different with hidden flag INTEG8 in /ANALY control card     
         DO i=1,nel
          f11(i)=f11(i)-(s1(i)*px1(i)+s4(i)*py1(i)+s6(i)*pz1(i))
          f21(i)=f21(i)-(s2(i)*py1(i)+s4(i)*px1(i)+s5(i)*pz1(i))
          f31(i)=f31(i)-(s3(i)*pz1(i)+s6(i)*px1(i)+s5(i)*py1(i))
          f12(i)=f12(i)-(s1(i)*px2(i)+s4(i)*py2(i)+s6(i)*pz2(i))
          f22(i)=f22(i)-(s2(i)*py2(i)+s4(i)*px2(i)+s5(i)*pz2(i))
          f32(i)=f32(i)-(s3(i)*pz2(i)+s6(i)*px2(i)+s5(i)*py2(i))
          f13(i)=f13(i)-(s1(i)*px3(i)+s4(i)*py3(i)+s6(i)*pz3(i))
          f23(i)=f23(i)-(s2(i)*py3(i)+s4(i)*px3(i)+s5(i)*pz3(i))
          f33(i)=f33(i)-(s3(i)*pz3(i)+s6(i)*px3(i)+s5(i)*py3(i))
          f14(i)=f14(i)-(s1(i)*px4(i)+s4(i)*py4(i)+s6(i)*pz4(i))
          f24(i)=f24(i)-(s2(i)*py4(i)+s4(i)*px4(i)+s5(i)*pz4(i))
          f34(i)=f34(i)-(s3(i)*pz4(i)+s6(i)*px4(i)+s5(i)*py4(i))
          f15(i)=f15(i)-(s1(i)*px5(i)+s4(i)*py5(i)+s6(i)*pz5(i))
          f25(i)=f25(i)-(s2(i)*py5(i)+s4(i)*px5(i)+s5(i)*pz5(i))
          f35(i)=f35(i)-(s3(i)*pz5(i)+s6(i)*px5(i)+s5(i)*py5(i))
          f16(i)=f16(i)-(s1(i)*px6(i)+s4(i)*py6(i)+s6(i)*pz6(i))
          f26(i)=f26(i)-(s2(i)*py6(i)+s4(i)*px6(i)+s5(i)*pz6(i))
          f36(i)=f36(i)-(s3(i)*pz6(i)+s6(i)*px6(i)+s5(i)*py6(i))
          f17(i)=f17(i)-(s1(i)*px7(i)+s4(i)*py7(i)+s6(i)*pz7(i))
          f27(i)=f27(i)-(s2(i)*py7(i)+s4(i)*px7(i)+s5(i)*pz7(i))
          f37(i)=f37(i)-(s3(i)*pz7(i)+s6(i)*px7(i)+s5(i)*py7(i))
          f18(i)=f18(i)-(s1(i)*px8(i)+s4(i)*py8(i)+s6(i)*pz8(i))
          f28(i)=f28(i)-(s2(i)*py8(i)+s4(i)*px8(i)+s5(i)*pz8(i))
          f38(i)=f38(i)-(s3(i)*pz8(i)+s6(i)*px8(i)+s5(i)*py8(i))
         ENDDO
        ENDIF
      ENDIF  
      
      RETURN
      END