stl_functions_des_mod.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: stl_functions_des                                      !
!  Author: Rahul Garg                                 Date: 24-Oct-13  !
!                                                                      !
!  Purpose: This module containd routines for geometric interaction    !
!  required for STL files.                                             !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      MODULE stl_functions_des

      IMPLICIT NONE

! Use this module only to define functions and subroutines.
      CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ClosestPtPointTriangle                                  !
!  Author: Rahul Garg                                 Date: 24-Oct-13  !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      Subroutine closestptpointtriangle(pointp, points, closest_point)
      USE param1, only: zero, one
      USE discretelement, only: dimn

      IMPLICIT NONE

! points are the three nodes of the triangle
! point p is the sphere center
      double precision, intent(in), dimension(3,3) :: points
      double precision, intent(in), dimension(dimn) :: pointp
      double precision, intent(out), dimension(dimn) ::  closest_point
! Local variables
      double precision, dimension(3) :: pointa, pointb, pointc
      double precision, dimension(dimn) :: ab, ac, ap, bp,cp
      double precision :: d1, d2, d3, d4, vc, v, d5, d6, vb, w, va, denom

      pointa = points(1,:)
      pointb = points(2,:)
      pointc = points(3,:)

      ab = pointb - pointa
      ac = pointc - pointa
      ap = pointp - pointa
      d1 = dot_product(ab, ap)
      d2 = dot_product(ac, ap)

      IF(d1 <= zero .AND. d2 <= zero) then
         closest_point = pointa
         return
      end if

! Check if P in vertex region outside B
      bp = pointp - pointb
      d3 = dot_product(ab, bp);
      d4 = dot_product(ac, bp);
      if (d3 >= zero .and. d4 <= d3) then
         closest_point = pointb
         return
      endif

! Check if P in edge region of AB, if so return projection of P onto AB
      vc = d1*d4 - d3*d2;
      if (vc <= zero .and. d1 >= zero .and. d3 <= zero) then
         v = d1 / (d1 - d3);
         closest_point =  pointa + v * ab;
         return
      end if

! Check if P in vertex region outside C
      cp = pointp - pointc
      d5 = dot_product(ab, cp)
      d6 = dot_product(ac, cp)
      if (d6 >= zero .and. d5 <= d6) then
         closest_point  = pointc
         return
      endif

! Check if P in edge region of AC, if so return projection of P onto AC
      vb = d5*d2 - d1*d6

      if (vb <= zero .and. d2 >= zero .and. d6 <= zero) then
         w = d2 / (d2 - d6)
         closest_point = pointa + w * ac
         return
      end if

! Check if P in edge region of BC, if so return projection of P onto BC
      va = d3*d6 - d5*d4
      if (va <= zero .and.(d4 - d3) >= zero .and. (d5 - d6) >= zero) then
         w = (d4 - d3) / ((d4 - d3) + (d5 - d6))
         closest_point = pointb + w * (pointc - pointb)
         return
      end if


! P inside face region. Compute Q through its barycentric coordinates (u,v,w)
      denom = one / (va + vb + vc)
      v = vb * denom
      w = vc * denom
      closest_point = pointa + ab * v + ac * w;
      return
      end Subroutine closestptpointtriangle



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: intersectLnPlane                                        !
!  Author: Rahul Garg                                 Date: 24-Oct-13  !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      Subroutine intersectlnplane(ref_line, dir_line, ref_plane,       &
         norm_plane, line_param)

      USE discretelement, only: dimn
      USE param1, only: zero

      IMPLICIT NONE

! Reference point and direction of the line
      DOUBLE PRECISION, INTENT(IN) :: REF_LINE(3),  DIR_LINE(3)
! reference point and normal of the plane
      DOUBLE PRECISION, INTENT(IN) :: REF_PLANE(3), NORM_PLANE(3)

! line is parameterized as p = p_ref + t * dir_line, t is line_param
      double precision, intent(out) :: line_param

      !local vars
      double precision :: denom

      denom = dot_product(dir_line, norm_plane)

      if(denom*denom.gt.zero) then
         line_param = dot_product(ref_plane-ref_line, norm_plane)
         line_param = line_param/denom
      endif

      return
      end subroutine intersectlnplane

!......................................................................!
! Subroutine TRI_BOX_OVERLAP                                           !
! Author: J.Musser                                   Date: 10-22-2015  !
!                                                                      !
! Purpose: Determine if a box (DES grid cell) intersects the triangle  !
!    (SLT). Note that the DES grid size is slightly increased to       !
!    capture STLs near the boarder of the cell. Otherwise, some        !
!    collisions could ve over looked.                                  !
!                                                                      !
! Author: Tomas Akenine-Moller                   Accessed: 10-22-2015  !
! REF: http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/    !
!         code/tribox2.txt                                             !
!......................................................................!
      SUBROUTINE tri_box_overlap(pCENTER, pHALFSIZE, pVERTS, pOVERLAP)

      IMPLICIT NONE

      DOUBLE PRECISION, INTENT(IN) :: pCENTER(3), pHALFSIZE(3)
      DOUBLE PRECISION, INTENT(IN) :: pVERTS(3,3)
      LOGICAL, INTENT(OUT) :: pOVERLAP

      DOUBLE PRECISION :: v0(3), v1(3), v2(3)
      DOUBLE PRECISION :: fex, fey, fez
      DOUBLE PRECISION :: normal(3), e0(3), e1(3), e2(3)

      poverlap = .false.

      v0 = pverts(1,:) - pcenter
      v1 = pverts(2,:) - pcenter
      v2 = pverts(3,:) - pcenter

      e0 = v1-v0
      e1 = v2-v1
      e2 = v0-v2

      fex = abs(e0(1))
      fey = abs(e0(2))
      fez = abs(e0(3))

      if(atest_x01(e0(3),e0(2),fez,fey)) return
      if(atest_y02(e0(3),e0(1),fez,fex)) return
      if(atest_z12(e0(2),e0(1),fey,fex)) return

      fex = abs(e1(1))
      fey = abs(e1(2))
      fez = abs(e1(3))

      if(atest_x01(e1(3),e1(2),fez,fey)) return
      if(atest_y02(e1(3),e1(1),fez,fex)) return
      if(atest_z0(e1(2),e1(1),fey,fex)) return

      fex = abs(e2(1))
      fey = abs(e2(2))
      fez = abs(e2(3))

      if(atest_x2(e2(3),e2(2),fez,fey)) return
      if(atest_y1(e2(3),e2(1),fez,fex)) return
      if(atest_z12(e2(2),e2(1),fey,fex)) return

      if(findmin(v0(1),v1(1),v2(1)) > phalfsize(1) .OR. &
         findmax(v0(1),v1(1),v2(1)) <-phalfsize(1)) return

      if(findmin(v0(2),v1(2),v2(2)) > phalfsize(2) .OR. &
         findmax(v0(2),v1(2),v2(2)) <-phalfsize(2)) return

      if(findmin(v0(3),v1(3),v2(3)) > phalfsize(3) .OR. &
         findmax(v0(3),v1(3),v2(3)) <-phalfsize(3)) return


      normal(1) = e0(2)*e1(3)-e0(3)*e1(2)
      normal(2) = e0(3)*e1(1)-e0(1)*e1(3)
      normal(3) = e0(1)*e1(2)-e0(2)*e1(1)

      if(.NOT.planeboxoverlap(normal,v0,phalfsize)) return

      poverlap = .true.

      RETURN

      CONTAINS

!``````````````````````````````````````````````````````````````````````!
! Function: planeBoxOverlap                                            !
! Support function for TRI_BOX_OVERLAP                                 !
!``````````````````````````````````````````````````````````````````````!
      LOGICAL FUNCTION planeboxoverlap(norm, vert, maxbox)

      double precision :: norm(3), vert(3), maxbox(3)

      integer :: lc
      double precision :: vmin(3), vmax(3), v

      do lc=1,3
         v=vert(lc)
         if(norm(lc) > 0.0d0) then
            vmin(lc) = -maxbox(lc) - v
            vmax(lc) =  maxbox(lc) - v
         else
            vmin(lc) = maxbox(lc) - v
            vmax(lc) =-maxbox(lc) - v
         endif
      enddo

      if(dot_product(norm,vmin) > 0.0d0) then
         planeboxoverlap = .false.
         return
      elseif(dot_product(norm,vmax) >= 0.0d0) then
         planeboxoverlap = .true.
         return
      endif

      planeboxoverlap = .false.
      return

      RETURN
      END FUNCTION planeboxoverlap

!``````````````````````````````````````````````````````````````````````!
! Function: findMin                                                    !
! Support function for TRI_BOX_OVERLAP                                 !
!``````````````````````````````````````````````````````````````````````!
      DOUBLE PRECISION FUNCTION findmin(x0,x1,x2)

      double precision :: x0,x1,x2

      findmin = x0

      if(x1<findmin) findmin=x1
      if(x2<findmin) findmin=x2

      RETURN
      END FUNCTION findmin

!``````````````````````````````````````````````````````````````````````!
! Function: findMax                                                    !
! Support function for TRI_BOX_OVERLAP                                 !
!``````````````````````````````````````````````````````````````````````!
      DOUBLE PRECISION FUNCTION findmax(x0,x1,x2)

      double precision :: x0,x1,x2

      findmax = x0

      if(x1>findmax) findmax=x1
      if(x2>findmax) findmax=x2

      RETURN
      END FUNCTION findmax

!``````````````````````````````````````````````````````````````````````!
! Function: ATEST_X01                                                  !
! Support function for TRI_BOX_OVERLAP                                 !
!``````````````````````````````````````````````````````````````````````!
      LOGICAL FUNCTION atest_x01(a,b,fa,fb)

      double precision :: a, b, fa, fb
      double precision :: lMin, lMax, p0, p2, rad

      p0 = a*v0(2) - b*v0(3)
      p2 = a*v2(2) - b*v2(3)

      if(p0<p2) then; lmin=p0; lmax=p2
      else; lmin=p2; lmax=p0; endif

      rad=fa*phalfsize(2) + fb*phalfsize(3)
      atest_x01=(lmin>rad .OR. lmax<-rad)

      END FUNCTION atest_x01

!``````````````````````````````````````````````````````````````````````!
! Function: ATEST_X2                                                   !
! Support function for TRI_BOX_OVERLAP                                 !
!``````````````````````````````````````````````````````````````````````!
      LOGICAL FUNCTION atest_x2(a,b,fa,fb)

      double precision :: a, b, fa, fb
      double precision :: lMin, lMax, p0, p1, rad

      p0 = a*v0(2) - b*v0(3)
      p1 = a*v1(2) - b*v1(3)

      if(p0<p1) then; lmin=p0; lmax=p1
      else; lmin=p1; lmax=p0; endif

      rad=fa*phalfsize(2) + fb*phalfsize(3)
      atest_x2=(lmin>rad .OR. lmax<-rad)

      END FUNCTION atest_x2

!``````````````````````````````````````````````````````````````````````!
! Function: ATEST_Y02                                                  !
! Support function for TRI_BOX_OVERLAP                                 !
!``````````````````````````````````````````````````````````````````````!
      LOGICAL FUNCTION atest_y02(a,b,fa,fb)

      double precision :: a, b, fa, fb
      double precision :: lMin, lMax, p0, p2, rad

      p0 = -a*v0(1) + b*v0(3)
      p2 = -a*v2(1) + b*v2(3)

      if(p0<p2) then; lmin=p0; lmax=p2
      else; lmin=p2; lmax=p0; endif

      rad=fa*phalfsize(1) + fb*phalfsize(3)
      atest_y02=(lmin>rad .OR. lmax<-rad)

      END FUNCTION atest_y02

!``````````````````````````````````````````````````````````````````````!
! Function: ATEST_Y1                                                   !
! Support function for TRI_BOX_OVERLAP                                 !
!``````````````````````````````````````````````````````````````````````!
      LOGICAL FUNCTION atest_y1(a,b,fa,fb)

      double precision :: a, b, fa, fb
      double precision :: lMin, lMax, p0, p1, rad

      p0 = -a*v0(1) + b*v0(3)
      p1 = -a*v1(1) + b*v1(3)

      if(p0<p1) then; lmin=p0; lmax=p1
      else; lmin=p1; lmax=p0; endif

      rad=fa*phalfsize(1) + fb*phalfsize(3)
      atest_y1=(lmin>rad .OR. lmax<-rad)

      END FUNCTION atest_y1

!``````````````````````````````````````````````````````````````````````!
! Function: ATEST_Z12                                                  !
! Support function for TRI_BOX_OVERLAP                                 !
!``````````````````````````````````````````````````````````````````````!
      LOGICAL FUNCTION atest_z12(a,b,fa,fb)

      double precision :: a, b, fa, fb
      double precision :: lMin, lMax, p1, p2, rad

      p1 = a*v1(1) - b*v1(2)
      p2 = a*v2(1) - b*v2(2)

      if(p2<p1) then; lmin=p2; lmax=p1
      else; lmin=p1; lmax=p2; endif

      rad=fa*phalfsize(1) + fb*phalfsize(2)
      atest_z12=(lmin>rad .OR. lmax<-rad)

      END FUNCTION atest_z12

!``````````````````````````````````````````````````````````````````````!
! Function: ATEST_Z0                                                   !
! Support function for TRI_BOX_OVERLAP                                 !
!``````````````````````````````````````````````````````````````````````!
      LOGICAL FUNCTION atest_z0(a,b,fa,fb)

      double precision :: a, b, fa, fb
      double precision :: lMin, lMax, p0, p1, rad

      p0 = a*v0(1) - b*v0(2)
      p1 = a*v1(1) - b*v1(2)

      if(p0<p1) then; lmin=p0; lmax=p1
      else; lmin=p1; lmax=p0; endif

      rad=fa*phalfsize(1) + fb*phalfsize(2)
      atest_z0=(lmin>rad .OR. lmax<-rad)

      END FUNCTION atest_z0

      END SUBROUTINE tri_box_overlap


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  SUBROUTINE: CHECK_IF_PARTICLE_OVERLAPS_STL                          !
!  Authors: Rahul Garg                               Date: 21-Mar-2014 !
!                                                                      !
!  Purpose: This subroutine is special written to check if a particle  !
!          overlaps any of the STL faces. The routine exits on         !
!          detecting an overlap. It is called after initial            !
!          generation of lattice configuration to remove out of domain !
!          particles                                                   !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE check_if_particle_overlaps_stl(POS, fI, fJ, fK, REMOVE)

      USE compar
      USE constant
      USE cutcell
      USE desgrid
      USE discretelement, ONLY: dimn, xe, yn, zt, max_radius
      USE functions
      USE geometry
      USE indices
      USE param1
      USE run
      USE stl

      Implicit none

      DOUBLE PRECISION, INTENT(IN) :: POS(dimn)
      INTEGER, INTENT(IN) :: fI, fJ, fK
      LOGICAL, INTENT(OUT) :: REMOVE

! Integers mapping the fluid cell corners to DES Grid indices.
      INTEGER :: I1, I2, J1, J2, K1, K2

      INTEGER I, J, K, IJK, NF, LC

      DOUBLE PRECISION :: LINE_T
      DOUBLE PRECISION :: RADSQ, DIST(3)

      double precision :: vv(3)

      remove = .true.

      i1 = iofpos(xe(fi-1))
      i2 = iofpos(xe( fi ))

      j1 = jofpos(yn(fj-1))
      j2 = jofpos(yn( fj ))

      k1 = kofpos(zt(fk-1))
      k2 = kofpos(zt( fk ))

      radsq = (1.05d0*max_radius)**2

      DO k = k1, k2
      DO j = j1, j2
      DO i = i1, i2

         IF(.NOT.dg_is_on_mype_plus1layers(i,j,k)) cycle

         ijk = dg_funijk(i,j,k)

! The point is on the non-fluid side of the plane if t>0
         DO lc = 1, facets_at_dg(ijk)%COUNT
            nf = facets_at_dg(ijk)%ID(lc)

            vv = vertex(1,:,nf)
            CALL intersectlnplane(pos, norm_face(:,nf), &
               vv, norm_face(:,nf), line_t)

! Orthogonal projection puts the point outside of the domain or less than
! one particle radius to the facet.
            dist = line_t*norm_face(:,nf)
            IF(line_t > zero .OR. dot_product(dist,dist)<=radsq) RETURN

         ENDDO
      ENDDO
      ENDDO
      ENDDO

      remove = .false.

      RETURN
      END SUBROUTINE check_if_particle_overlaps_stl

      END MODULE stl_functions_des