s4therm-itet1.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/.
!||====================================================================
!||    s4therm_itet1   ../engine/source/elements/solid/solide4/s4therm-itet1.F
!||--- called by ------------------------------------------------------
!||    s10forc3        ../engine/source/elements/solid/solide10/s10forc3.F
!||====================================================================
      SUBROUTINE s4therm_itet1(PM   ,IMAT  ,NC   ,NEL  ,
     .                         XX   ,YY   ,ZZ   ,DT1  ,HEAT ,
     .                         TEMP ,FPHI ,OFFG ,OFF  ,THEACCFACT)
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"
#include      "param_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NEL, NC(MVSIZ,10)
      INTEGER, INTENT(IN) :: IMAT
      my_real, INTENT(IN) :: theaccfact
      my_real, INTENT(IN) :: dt1 
      my_real :: temp(*), heat(*), fphi(mvsiz,10), pm(npropm,*), 
     .    off(*), offg(*)
      double precision
     .    xx(mvsiz,10), yy(mvsiz,10), zz(mvsiz,10)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, J
      my_real
     .   det(nel),
     .   px1(nel), px2(nel), px3(nel), px4(nel),  
     .   py1(nel), py2(nel), py3(nel), py4(nel),  
     .   pz1(nel), pz2(nel), pz3(nel), pz4(nel)
      my_real
     .   b1(nel), b2(nel), b3(nel), b4(nel), 
     .   c1(nel), c2(nel), c3(nel), c4(nel),
     .   d1(nel), d2(nel), d3(nel), d4(nel)
      my_real
     .   x41, y41, z41, x42, y42, z42, x43, y43, z43
      my_real
     .  ca, cb, kc, phix, phiy, phiz, a, d,
     .  temp1, temp2, temp3, temp4, tempel
C
C---+----1
      DO i=1,nel
       x43 = xx(i,4) - xx(i,3)
       y43 = yy(i,4) - yy(i,3)
       z43 = zz(i,4) - zz(i,3)
       x41 = xx(i,4) - xx(i,1)
       y41 = yy(i,4) - yy(i,1)
       z41 = zz(i,4) - zz(i,1)
       x42 = xx(i,4) - xx(i,2)
       y42 = yy(i,4) - yy(i,2)
       z42 = zz(i,4) - zz(i,2)
C
       b1(i) =  y43*z42 - y42*z43
       b2(i) =  y41*z43 - y43*z41
       b3(i) =  y42*z41 - y41*z42
       b4(i) =  -(b1(i) + b2(i) + b3(i))
C
       c1(i) =  z43*x42 - z42*x43
       c2(i) =  z41*x43 - z43*x41
       c3(i) =  z42*x41 - z41*x42
       c4(i) =  -(c1(i) + c2(i) + c3(i))
C
       d1(i) =  x43*y42 - x42*y43
       d2(i) =  x41*y43 - x43*y41
       d3(i) =  x42*y41 - x41*y42
       d4(i) =  -(d1(i) + d2(i) + d3(i))
C
       det(i) = (x41*b1(i) + y41*c1(i) + z41*d1(i))*one_over_6
      ENDDO
 
C---+----1
      DO i=1,nel
        d = one/det(i)/six
        px1(i)=-b1(i)*d
        py1(i)=-c1(i)*d
        pz1(i)=-d1(i)*d
        px2(i)=-b2(i)*d
        py2(i)=-c2(i)*d
        pz2(i)=-d2(i)*d
        px3(i)=-b3(i)*d
        py3(i)=-c3(i)*d
        pz3(i)=-d3(i)*d
        px4(i)=-b4(i)*d
        py4(i)=-c4(i)*d
        pz4(i)=-d4(i)*d
      ENDDO
 
C---+----1
      ca = pm(75,imat)
      cb = pm(76,imat)
C
      DO i=1,nel
         IF(off(i)==zero.OR.offg(i)<=zero) cycle
         temp1 = temp(nc(i,1))
         temp2 = temp(nc(i,2))
         temp3 = temp(nc(i,3))
         temp4 = temp(nc(i,4))
C 
C - flux
C      
         phix =  temp1*px1(i) + temp2*px2(i) + temp3*px3(i) + temp4*px4(i)                
         phiy =  temp1*py1(i) + temp2*py2(i) + temp3*py3(i) + temp4*py4(i)       
         phiz =  temp1*pz1(i) + temp2*pz2(i) + temp3*pz3(i) + temp4*pz4(i)        
C
         tempel = fourth * (temp1 + temp2 + temp3 + temp4)
         kc = (ca + cb*tempel)*dt1*det(i)*theaccfact
         phix = kc*phix
         phiy = kc*phiy
         phiz = kc*phiz
C
C force thermique nodale
C
         a = fourth * heat(i)
         fphi(i,1) = a - (phix*px1(i) + phiy*py1(i) + pz1(i)*phiz)
         fphi(i,2) = a - (phix*px2(i) + phiy*py2(i) + pz2(i)*phiz)
         fphi(i,3) = a - (phix*px3(i) + phiy*py3(i) + pz3(i)*phiz)
         fphi(i,4) = a - (phix*px4(i) + phiy*py4(i) + pz4(i)*phiz)
      ENDDO
C 
      RETURN
      END