intanl_qd.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_qd       ../starter/source/fluid/intanl_qd.F
!||--- called by ------------------------------------------------------
!||    mass_fluid_qd   ../starter/source/fluid/mass-fluid_qd.F
!||====================================================================
      SUBROUTINE intanl_qd(X1,  Y1,  Z1,  X2,  Y2,  Z2,  X3,  Y3,  Z3, 
     .                     X4,  Y4,  Z4,  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, x4, y4, z4, 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(5), eta(5), dksi(4), deta(4), r(5),
     .        xls, yls, zls, s(4), v, fln, arg,
     .        d2, l12, l22, l32, l42, lm2
      my_real cs(4), sn(4)
C
      x0=fourth*(x1+x2+x3+x4)
      y0=fourth*(y1+y2+y3+y4)
      z0=fourth*(z1+z2+z3+z4)
C
C SIMPLIFICATION SI SOURCE LOIN DU POINT P
      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=(x4-x3)**2+(y4-y3)**2+(z4-z3)**2
      l42=(x1-x4)**2+(y1-y4)**2+(z1-z4)**2
      lm2=max(l12,l22,l32,l42)
 
      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 COORDONNEES LOCALES
        vx1=x2-x1
        vy1=y2-y1
        vz1=z2-z1
        vx2=x4-x1
        vy2=y4-y1
        vz2=z4-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)=(x4-x0)*vx1+(y4-y0)*vy1+(z4-z0)*vz1
        eta(4)=(x4-x0)*vx2+(y4-y0)*vy2+(z4-z0)*vz2
        ksi(5)=ksi(1)
        eta(5)=eta(1)
 
        dksi(1)=ksi(2)-ksi(1)
        dksi(2)=ksi(3)-ksi(2)
        dksi(3)=ksi(4)-ksi(3)
        dksi(4)=ksi(1)-ksi(4)
        deta(1)=eta(2)-eta(1)
        deta(2)=eta(3)-eta(2)
        deta(3)=eta(4)-eta(3)
        deta(4)=eta(1)-eta(4)
        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)
        r(4)=sqrt((xp-x4)**2+(yp-y4)**2+(zp-z4)**2)
        s(1)=sqrt(l12)
        s(2)=sqrt(l22)
        s(3)=sqrt(l32)
        s(4)=sqrt(l42)
        r(5)=r(1)
            
        DO i=1,4
          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,4
              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,4
           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
           ELSE
              WRITE(6,'(A,E13.5)') '? ARG=',arg
           ENDIF
        ENDDO
        rvlg=-rvlg+zls*rvlh
 
      ENDIF   
C
      RETURN
      END