dist_node_seg3n.F File Reference

#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	dist_node_seg3n (dist, dmin, dmax, nod_x, nod_y, nod_z, ax, ay, az, bx, by, bz, cx, cy, cz)

Function/Subroutine Documentation

dist_node_seg3n()

subroutine dist_node_seg3n	(	intent(out)	dist,
		intent(in)	dmin,
		intent(in)	dmax,
		intent(in)	nod_x,
		intent(in)	nod_y,
		intent(in)	nod_z,
		intent(in)	ax,
		intent(in)	ay,
		intent(in)	az,
		intent(in)	bx,
		intent(in)	by,
		intent(in)	bz,
		intent(in)	cx,
		intent(in)	cy,
		intent(in)	cz )

Definition at line 28 of file dist_node_seg3n.F.

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