intanl_tg.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/.
!||====================================================================
!||    intanl_tg       ../starter/source/fluid/intanl_tg.f
!||--- called by ------------------------------------------------------
!||    mass_fluid_qd   ../starter/source/fluid/mass-fluid_qd.F
!||    mass_fluid_tg   ../starter/source/fluid/mass-fluid_tg.F
!||====================================================================
      SUBROUTINE intanl_tg(X1  , Y1  , Z1  , X2 , Y2 ,
     .                     Z2  , X3  , Y3  , Z3 , XP ,
     .                     YP  , ZP  , NRX , NRY, NRZ,
     .                     AREA, RVLH, RVLG, JEL, IEL)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IEL, JEL
      my_real
     .        X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3, XP, YP, ZP,
     .        nrx, nry, nrz, area, rvlh, rvlg
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, J
      my_real
     .        VX1, VY1, VZ1, VX2, VY2, VZ2, S1, S12, S2, NR1, NR2, 
     .        x0, y0, z0, ksi(4), eta(4), dksi(3), deta(3), r(4),
     .        xls, yls, zls, s(3), v, fln, arg,
     .        d2, l12, l22, l32, lm2
      my_real cs(3), sn(3)
 
C
      x0=third*(x1+x2+x3)
      y0=third*(y1+y2+y3)
      z0=third*(z1+z2+z3)
C
C SIMPLIFICATION SI SOURCE LOIN DE L'ELEMENT
      d2=(x0-xp)**2+(y0-yp)**2+(z0-zp)**2
      l12=(x2-x1)**2+(y2-y1)**2+(z2-z1)**2
      l22=(x3-x2)**2+(y3-y2)**2+(z3-z2)**2
      l32=(x1-x3)**2+(y1-y3)**2+(z1-z3)**2
      lm2=max(l12,l22,l32)
 
      IF(d2>twenty5*lm2) THEN  
        rvlg=area/sqrt(d2)
        rvlh=area*(nrx*(x0-xp)+nry*(y0-yp)+nrz*(z0-zp))/(d2**three_half)    
      ELSE
C
C COORDONNEES LOCALES
        vx1=x2-x1
        vy1=y2-y1
        vz1=z2-z1
        vx2=x3-x1
        vy2=y3-y1
        vz2=z3-z1
C
        s1=vx1*vx1+vy1*vy1+vz1*vz1
        s12=vx1*vx2+vy1*vy2+vz1*vz2
        nr1=sqrt(s1)
C     
        vx2=vx2-s12/s1*vx1
        vy2=vy2-s12/s1*vy1
        vz2=vz2-s12/s1*vz1
C
        s2=vx2*vx2+vy2*vy2+vz2*vz2
        nr2=sqrt(s2)
        vx1=vx1/nr1
        vy1=vy1/nr1
        vz1=vz1/nr1
        vx2=vx2/nr2
        vy2=vy2/nr2
        vz2=vz2/nr2
C
        xls=(xp-x0)*vx1+(yp-y0)*vy1+(zp-z0)*vz1
        yls=(xp-x0)*vx2+(yp-y0)*vy2+(zp-z0)*vz2
        zls=(xp-x0)*nrx+(yp-y0)*nry+(zp-z0)*nrz
 
        ksi(1)=(x1-x0)*vx1+(y1-y0)*vy1+(z1-z0)*vz1
        eta(1)=(x1-x0)*vx2+(y1-y0)*vy2+(z1-z0)*vz2
        ksi(2)=(x2-x0)*vx1+(y2-y0)*vy1+(z2-z0)*vz1
        eta(2)=(x2-x0)*vx2+(y2-y0)*vy2+(z2-z0)*vz2
        ksi(3)=(x3-x0)*vx1+(y3-y0)*vy1+(z3-z0)*vz1
        eta(3)=(x3-x0)*vx2+(y3-y0)*vy2+(z3-z0)*vz2
        ksi(4)=ksi(1)
        eta(4)=eta(1)
 
        dksi(1)=ksi(2)-ksi(1)
        dksi(2)=ksi(3)-ksi(2)
        dksi(3)=ksi(1)-ksi(3)
        deta(1)=eta(2)-eta(1)
        deta(2)=eta(3)-eta(2)
        deta(3)=eta(1)-eta(3)
        r(1)=sqrt((xp-x1)**2+(yp-y1)**2+(zp-z1)**2)
        r(2)=sqrt((xp-x2)**2+(yp-y2)**2+(zp-z2)**2)
        r(3)=sqrt((xp-x3)**2+(yp-y3)**2+(zp-z3)**2)
        s(1)=sqrt(l12)
        s(2)=sqrt(l22)
        s(3)=sqrt(l32)
        r(4)=r(1)
            
        DO i=1,3
          cs(i)=dksi(i)/s(i)
          sn(i)=deta(i)/s(i)
        ENDDO
C
        rvlh=zero
C INTEGRALE DOUBLE COUCHE
        IF (zls/=zero) THEN
           DO i=1,3
              j=i+1
              rvlh=rvlh
     .       +atan((deta(i)*((xls-ksi(i))**2+zls**2)-dksi(i)*(xls-ksi(i))*(yls-eta(i)))/(r(i)*zls*dksi(i)))
     .       -atan((deta(i)*((xls-ksi(j))**2+zls**2)-dksi(i)*(xls-ksi(j))*(yls-eta(j)))/(r(j)*zls*dksi(i)))
           ENDDO
        ENDIF
C
C INTEGRALE SIMPLE COUCHE
        rvlg=zero
        DO i=1,3
           j=i+1
           v=(xls-ksi(i))*sn(i)-(yls-eta(i))*cs(i)
           arg=(r(i)+r(j)-s(i))/(r(i)+r(j)+s(i))
           IF (arg>zero) THEN
              fln=-log(arg)
              rvlg=rvlg+v*fln
           ENDIF
        ENDDO
        rvlg=-rvlg+zls*rvlh
C
      ENDIF
C
      RETURN
      END