dist_node_seg3n.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/.
!||====================================================================
!||    dist_node_seg3n     ../engine/source/tools/sensor/dist_node_seg3n.F
!||--- called by ------------------------------------------------------
!||    sensor_dist_surf0   ../engine/source/tools/sensor/sensor_dist_surf0.F
!||====================================================================
        SUBROUTINE dist_node_seg3n (
     .             DIST,DMIN,DMAX,NOD_X,NOD_Y,NOD_Z,
     .             AX,AY,AZ,BX,BY,BZ,CX,CY,CZ)
c-----------------------------------------------------------------------
c Calculates distance between NODE(NOD_X,NOD_Y,NOD_Z) and a 3 node segment (A-B-C)
c returns distance to plane if node projection is inside the segment
c otherwise : check closest distance to edges and segment points
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-----------------------------------------------
      my_real , INTENT(IN)  :: dmin,dmax,nod_x,nod_y,nod_z,
     .                         ax,ay,az,bx,by,bz,cx,cy,cz
      my_real , INTENT(OUT) :: dist
C----------------------------------------------------------
C Local Variables
C----------------------------------------------------------
      my_real :: abx,aby,abz,acx,acy,acz,bcx,bcy,bcz,vecx,vecy,vecz,nx,ny,nz,
     .   ux,uy,uz,vx,vy,vz,wx,wy,wz,qx,qy,qz,qa_x,qa_y,qa_z,qb_x,qb_y,qb_z,
     .   qc_x,qc_y,qc_z,area,a1,a2,sx,sy,sz,alpha,alpha2,d1,d2,d3,len,norm
C=======================================================================
c     vect AB = P1-P2
c     vect AC = P1-P3
      abx = bx - ax
      aby = by - ay
      abz = bz - az
      acx = cx - ax
      acy = cy - ay
      acz = cz - az
      bcx = cx - bx
      bcy = cy - by
      bcz = cz - bz
      ux = nod_x - ax
      uy = nod_y - ay
      uz = nod_z - az
c     normal to segment plane
      nx = aby*acz - abz*acy
      ny = abz*acx - abx*acz
      nz = abx*acy - aby*acx
      norm = sqrt(nx*nx + ny*ny + nz*nz)
c     distance to plane
      dist = abs(nx*ux + ny*uy + nz*uz) / norm
      
      IF (dist > dmax .or. dist < dmin) RETURN  ! distance is outside of valid range
c---------------------------------------------------
c     Check if projection on segment plane is inside
c     Q = projection of NODE on segment plane using normal(NPX,NPY,NPZ)
c---------------------------------------------------
      qx = nod_x - dist*nx / norm
      qy = nod_y - dist*ny / norm
      qz = nod_z - dist*nz / norm
      qa_x = ax - qx
      qa_y = ay - qy
      qa_z = az - qz
      qb_x = bx - qx
      qb_y = by - qy
      qb_z = bz - qz
      qc_x = cx - qx
      qc_y = cy - qy
      qc_z = cz - qz
c
c     Surfaces de sous-triangles
c
      sx = qa_y * qb_z - qa_z * qb_y
      sy = qa_z * qb_x - qa_x * qb_z
      sz = qa_x * qb_y - qa_y * qb_x
      a1 = sqrt(sx*sx + sy*sy + sz*sz) * half
      sx = qc_y * qa_z - qc_z * qa_y
      sy = qc_z * qa_x - qc_x * qa_z
      sz = qc_x * qa_y - qc_y * qa_x
      a2 = sqrt(sx*sx + sy*sy + sz*sz) * half
      area = half*norm   ! total segment area
c
      IF (a1 + a2 > area) THEN
        ! Point projection on plane is outside the segment
        ! calculate distance to segment nodes
        d1 = sqrt(ux**2 + uy**2 + uz**2)
        vx = nod_x - bx
        vy = nod_y - by
        vz = nod_z - bz
        d2   = sqrt(vx**2 + vy**2 + vz**2)
        dist = min(d1, d2)
        wx = nod_x - cx
        wy = nod_y - cy
        wz = nod_z - cz
        d3   = sqrt(wx**2 + wy**2 + wz**2)
        dist = min(dist, d2)
        ! calculate distance to edges
        len   = sqrt(abx**2 + aby**2 + abz**2)
        alpha = (ux*abx + uy*aby + uz*abz) / len
        IF (alpha > 0 .and. alpha < one) THEN
          qx = ax + alpha * abx
          qy = ay + alpha * aby
          qz = az + alpha * abz
          d1 = sqrt((nod_x - qx)**2 + (nod_y - qy)**2 + (nod_z - qz)**2)
          dist = min(dist, d1)
        END IF
c
        len   = sqrt(bcx**2 + bcy**2 + bcz**2)
        alpha = (vx*bcx + vy*bcy + vz*bcz) / len
        IF (alpha > 0 .and. alpha < one) THEN
          qx = bx + alpha * bcx
          qy = by + alpha * bcy
          qz = bz + alpha * bcz
          d2 = sqrt((nod_x - qx)**2 + (nod_y - qy)**2 + (nod_z - qz)**2)
          dist = min(dist, d2)
        END IF
c
        len   = sqrt(acx**2 + acy**2 + acz**2)
        alpha = (ux*acx + uy*acy + uz*acz) / len
        IF (alpha > 0 .and. alpha < one) THEN
          qx = ax + alpha * acx
          qy = ay + alpha * acy
          qz = az + alpha * acz
          d3 = sqrt((nod_x - qx)**2 + (nod_y - qy)**2 + (nod_z - qz)**2)
          dist = min(dist, d3)
        END IF
c        
      END IF
c-----------
      RETURN
      END SUBROUTINE