inist3.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/.
!||====================================================================
!||    inist3    ../starter/source/interfaces/inter3d1/inist3.F
!||--- called by ------------------------------------------------------
!||    i12tid3   ../starter/source/interfaces/inter3d1/i12tid3.F
!||    i1tid3    ../starter/source/interfaces/inter3d1/i1tid3.F
!||    i3pen3    ../starter/source/interfaces/inter3d1/i3pen3.F
!||    i5pwr3    ../starter/source/interfaces/inter3d1/i3pen3.F
!||    i6pen3    ../starter/source/interfaces/inter3d1/i6pen3.F
!||====================================================================
      SUBROUTINE inist3(N1,N2,N3,SSC,TTC,IER,ALP,XX1,XX2,XX3,XS1,YS1,ZS1,XC,YC,ZC)
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,INTENT(INOUT)  :: IER
      my_real, INTENT(INOUT) :: n1, n2, n3, ssc, ttc
      my_real, INTENT(IN)    :: alp
      my_real,INTENT(INOUT)  :: xx1(4), xx2(4),xx3(4),xs1,ys1,zs1,xc,yc,zc
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      my_real
     .   h(4), x0, y0, z0, xl1, xl2, xl3, xl4, yy1, yy2, yy3, yy4,
     .   zz1, zz2, zz3, zz4, xi1, xi2, xi3, xi4, yi1, yi2, yi3, yi4,
     .   zi1, zi2, zi3, zi4, xn1, yn1, zn1, xn2, yn2, zn2, xn3, yn3,
     .   zn3, xn4, yn4, zn4, an, area, a12, a23, a34, a41, b12, b23,
     .   b34, b41, ab1, ab2, tp, tm, sp, sm, x1,x2,x3,x4,
     .   y1,y2,y3,y4,z1,z2,z3,z4,xi,yi,zi
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      !SECONDARY node coordinates
      xi = xs1
      yi = ys1
      zi = zs1
      
      !MAIN segment coordinates
      x1 = xx1(1)
      x2 = xx1(2)
      x3 = xx1(3)
      x4 = xx1(4)
      y1 = xx2(1)
      y2 = xx2(2)
      y3 = xx2(3)
      y4 = xx2(4)
      z1 = xx3(1)
      z2 = xx3(2)
      z3 = xx3(3)
      z4 = xx3(4)
      !MAIN segment centroid
      x0 = fourth*(x1+x2+x3+x4)
      y0 = fourth*(y1+y2+y3+y4)
      z0 = fourth*(z1+z2+z3+z4)
      !length from centroid
      xl1 = x1-x0
      xl2 = x2-x0
      xl3 = x3-x0
      xl4 = x4-x0      
      yy1 = y1-y0
      yy2 = y2-y0
      yy3 = y3-y0
      yy4 = y4-y0
      zz1 = z1-z0
      zz2 = z2-z0
      zz3 = z3-z0
      zz4 = z4-z0
      !length from SECONDARY node
      xi1 = x1-xi
      xi2 = x2-xi
      xi3 = x3-xi
      xi4 = x4-xi
      yi1 = y1-yi
      yi2 = y2-yi
      yi3 = y3-yi
      yi4 = y4-yi
      zi1 = z1-zi
      zi2 = z2-zi
      zi3 = z3-zi
      zi4 = z4-zi
      
      !1st mater triangle
      xn1 = yy1*zz2 - yy2*zz1
      yn1 = zz1*xl2 - zz2*xl1
      zn1 = xl1*yy2 - xl2*yy1
      n1  = xn1
      n2  = yn1
      n3  = zn1
      !2nd MAIN triangle
      xn2 = yy2*zz3 - yy3*zz2
      yn2 = zz2*xl3 - zz3*xl2
      zn2 = xl2*yy3 - xl3*yy2
      n1  = n1+xn2
      n2  = n2+yn2
      n3  = n3+zn2
      !3rd MAIN triangle
      xn3 = yy3*zz4 - yy4*zz3
      yn3 = zz3*xl4 - zz4*xl3
      zn3 = xl3*yy4 - xl4*yy3
      n1  = n1+xn3
      n2  = n2+yn3
      n3  = n3+zn3
      !4th MAIN triangle
      xn4 = yy4*zz1 - yy1*zz4
      yn4 = zz4*xl1 - zz1*xl4
      zn4 = xl4*yy1 - xl1*yy4
      n1  = n1+xn4
      n2  = n2+yn4
      n3  = n3+zn4
      !average normal of MAIN segment
      an  = max(em20,sqrt(n1*n1+n2*n2+n3*n3))
      n1  = n1/an
      n2  = n2/an
      n3  = n3/an
      !in case of degenerated segment
      IF(an<=em19) THEN
        ier=-1
        RETURN
      ENDIF
      area=half*an
      a12 = (n1*xn1+n2*yn1+n3*zn1)
      a23 = (n1*xn2+n2*yn2+n3*zn2)
      a34 = (n1*xn3+n2*yn3+n3*zn3)
      a41 = (n1*xn4+n2*yn4+n3*zn4)
 
      !1st triangle based on SECONDARY node
      xn1 = yi1*zi2 - yi2*zi1
      yn1 = zi1*xi2 - zi2*xi1
      zn1 = xi1*yi2 - xi2*yi1
      b12 = (n1*xn1+n2*yn1+n3*zn1)
      !2nd triangle based on SECONDARY node
      xn2 = yi2*zi3 - yi3*zi2
      yn2 = zi2*xi3 - zi3*xi2
      zn2 = xi2*yi3 - xi3*yi2
      b23 = (n1*xn2+n2*yn2+n3*zn2)
      !3rd triangle based on SECONDARY node
      xn3 = yi3*zi4 - yi4*zi3
      yn3 = zi3*xi4 - zi4*xi3
      zn3 = xi3*yi4 - xi4*yi3
      b34 = (n1*xn3+n2*yn3+n3*zn3)
      !4th triangle based on SECONDARY node
      xn4 = yi4*zi1 - yi1*zi4
      yn4 = zi4*xi1 - zi1*xi4
      zn4 = xi4*yi1 - xi1*yi4
      b41 = (n1*xn4+n2*yn4+n3*zn4)
      !SSC
      ab1 = a23*b41
      ab2 = b23*a41
      IF(abs(ab1+ab2)/area>em10)THEN
       ssc = (ab1-ab2)/(ab1+ab2)
      ELSE
       ssc = zero
      ENDIF
      !TTC
      IF(abs(a34/area)>em10)THEN
       ab1 = b12*a34
       ab2 = b34*a12
       ttc = (ab1-ab2)/(ab1+ab2)
      ELSE
       ttc = (b12-a12)/a12
       IF(b23<=zero .AND. b41<=zero)THEN
         IF(-b23/a12<=alp.AND.-b41/a12<=alp)ssc=zero
       ELSEIF(b23<=zero)THEN
         IF(-b23/a12<=alp)ssc=one
       ELSEIF(b41<=zero)THEN
         IF(-b41/a12<=alp)ssc=-one
       ENDIF
      ENDIF
      !lecipari : GAP=TWO*EM02 GAP->FRIGAP(2)->ALP (hardcoded tolerance)
 
      IF(abs(ssc)>one+alp.OR.abs(ttc)>one+alp) THEN
        ier = 1  !SECONDARY node not on MAIN segment
      ELSE
        ier = 0
        IF(abs(ssc)>one)ssc=ssc/abs(ssc)
        IF(abs(ttc)>one)ttc=ttc/abs(ttc)
      ENDIF
 
      tp   = fourth*(one+ttc)
      tm   = fourth*(one-ttc)
      sp   = one+ssc
      sm   = one-ssc
      h(1) = tm*sm
      h(2) = tm*sp
      h(3) = tp*sp
      h(4) = tp*sm
      xc   = h(1)*x1+h(2)*x2+h(3)*x3+h(4)*x4
      yc   = h(1)*y1+h(2)*y2+h(3)*y3+h(4)*y4
      zc   = h(1)*z1+h(2)*z2+h(3)*z3+h(4)*z4
 
      RETURN
      END