Commit 600c4cf5 authored by Kochunas, Brendan M.'s avatar Kochunas, Brendan M.
Browse files

Merge branch 'hex_geom_MR' into 'master' 

Modifications for MPACT hex geometry feature

See merge request https://code.ornl.gov/futility/Futility/-/merge_requests/400
parents 3b430733 300eee18

src/Geom.f90

+1 −0
Original line number Diff line number Diff line
@@ -49,6 +49,7 @@ PUBLIC :: OPERATOR(==) 
PUBLIC :: OPERATOR(/=)

PUBLIC :: OPERATOR(.APPROXEQA.)

PUBLIC :: ASSIGNMENT(=)

PUBLIC :: SOFTEQ



INTERFACE newGeom

  !> @copybrief Geom::newGeom_line

src/Geom_CircCyl.f90

+158 −23
Original line number Diff line number Diff line
@@ -222,10 +222,6 @@ ELEMENTAL SUBROUTINE set_CircleType(circ,p,r,angstt,angstp) 
      circ%thetastt=ZERO

      circ%thetastp=TWOPI

    ENDIF

    !circ%cosstt=COS(circ%thetastt)

    !circ%sinstt=SIN(circ%thetastt)

    !circ%cosstp=COS(circ%thetastp)

    !circ%sinstp=SIN(circ%thetastp)

  ENDIF

ENDSUBROUTINE set_CircleType

!
@@ -332,6 +328,7 @@ ENDSUBROUTINE clear_CylinderShellType 
!> @param line the line type thats being tested against the circle

!> @param p1 the first point of intersection (if it exists)

!> @param p2 the second point of intersection (if it exists)

!> @param checkRad logical to check radius line in case of partial circle

!> @note a return code is assigned to @c p1%dim and @c p2%dim indicating the type of

!> intersection. @n

!>  > 0: success; an intersection point was found @n
@@ -339,13 +336,26 @@ ENDSUBROUTINE clear_CylinderShellType 
!>   -1: problem with dummy arguments passed to routine @n

!>   -2: line is not directed toward circle (disjoint) @n

!>   -3: line segment is tangent @n

ELEMENTAL SUBROUTINE intersect_CircleType_and_LineType(circle,line,p1,p2)

ELEMENTAL SUBROUTINE intersect_CircleType_and_LineType(circle,line,p1,p2,checkRad)

  CLASS(CircleType),INTENT(IN) :: circle

  TYPE(LineType),INTENT(IN) :: line

  TYPE(PointType),INTENT(INOUT) :: p1

  TYPE(PointType),INTENT(INOUT) :: p2

  REAL(SRK) :: a,b,c,t1,t2,u(2),w(2),ra,discr

  LOGICAL(SBK),INTENT(IN),OPTIONAL :: checkRad

  LOGICAL(SBK) :: checkRadLine

  REAL(SRK) :: a,b,c,t1,t2,u(2),w(2),ra,discr,theta

  TYPE(LineType),ALLOCATABLE :: rline(:)

  TYPE(PointType),ALLOCATABLE :: rp(:)

  LOGICAL(SBK) :: checkTheta

  REAL(SRK),ALLOCATABLE :: thetaRange(:,:)

  REAL(SRK) :: thetastt,thetastp

  INTEGER(SIK) :: it



  IF(PRESENT(checkRad)) THEN

    checkRadLine=checkRad

  ELSE

    checkRadLine=.FALSE.

  ENDIF

  CALL p1%clear()

  CALL p2%clear()

  p1%dim=-1
@@ -359,6 +369,61 @@ ELEMENTAL SUBROUTINE intersect_CircleType_and_LineType(circle,line,p1,p2) 
    b=w(1)*u(1)+w(2)*u(2)

    c=w(1)*w(1)+w(2)*w(2)-circle%r*circle%r



    checkTheta=.FALSE.

    ! 1. check if it is needed to check partial angled circle

    ! 2. check if target line intersects radial lines of partial circle.

    IF(checkRadLine .AND. ISBETWEEN(ZERO,(circle%thetastp-circle%thetastt),TWOPI,EPSD)) THEN

      thetastt=circle%thetastt

      thetastp=circle%thetastp

      DO WHILE(thetastt >= TWOPI)

        thetastt=thetastt-TWOPI

        thetastp=thetastp-TWOPI

      ENDDO

      DO WHILE(thetastt < 0.0_SRK)

        thetastt=thetastt+TWOPI

        thetastp=thetastp+TWOPI

      ENDDO

      IF(thetastp > TWOPI) THEN

        ALLOCATE(thetaRange(2,2))

        thetaRange(1,1)=thetastt

        thetaRange(2,1)=TWOPI

        thetaRange(1,2)=0.0_SRK

        thetaRange(2,2)=thetastp-TWOPI

      ELSE

        ALLOCATE(thetaRange(2,1))

        thetaRange(1,1)=thetastt

        thetaRange(2,1)=thetastp

      ENDIF

      checkTheta=.TRUE.



      ALLOCATE(rline(2))

      ALLOCATE(rp(2))

      CALL rp(1)%init(DIM=2,X=0.0_SRK,Y=0.0_SRK)

      ! define lines in radial direction

      rp(1)%coord(1)=circle%c%coord(1)+ABS(circle%r)*COS(circle%thetastt)

      rp(1)%coord(2)=circle%c%coord(2)+ABS(circle%r)*SIN(circle%thetastt)

      CALL rline(1)%set(circle%c,rp(1))

      rp(1)%coord(1)=circle%c%coord(1)+ABS(circle%r)*COS(circle%thetastp)

      rp(1)%coord(2)=circle%c%coord(2)+ABS(circle%r)*SIN(circle%thetastp)

      CALL rline(2)%set(circle%c,rp(1))

      ! check if rline intersect the test line, if so add the intersection to rp

      rp(1)=rline(1)%intersect(line)

      IF(rp(1)%dim /= 2) THEN

        rp(1)=rline(2)%intersect(line)

        CALL rp(2)%clear()

      ELSE

        rp(2)=rline(2)%intersect(line)

      ENDIF

      CALL rline(1)%clear()

      CALL rline(2)%clear()

      DEALLOCATE(rline)

      IF(rp(1)%dim /= 2) THEN

        CALL rp(1)%clear()

        CALL rp(2)%clear()

        DEALLOCATE(rp)

      ENDIF

    ENDIF



    IF(c > zero .AND. b > zero) THEN

      p1%dim=-2

      p2%dim=-2
@@ -371,7 +436,6 @@ ELEMENTAL SUBROUTINE intersect_CircleType_and_LineType(circle,line,p1,p2) 
        p2%dim=-2

      ELSEIF(discr .APPROXEQA. zero) THEN

        !Tangent



        ra=one/a

        t1=-ra*b

        IF(t1<ZERO .OR. t1 > ONE ) THEN
@@ -395,13 +459,55 @@ ELEMENTAL SUBROUTINE intersect_CircleType_and_LineType(circle,line,p1,p2) 
          p1=line%p1

          p1%coord(1)=p1%coord(1)+u(1)*t1

          p1%coord(2)=p1%coord(2)+u(2)*t1

          IF(checkTheta) THEN

            ! if the point is on the circle, set %dim to be 2. Otherwise set dummy value.

            theta=ATAN2PI((p1%coord(1)-circle%c%coord(1)),(p1%coord(2)-circle%c%coord(2)))

            p1%dim=-5

            DO it=1,SIZE(thetaRange,DIM=2)

              IF(ISBETWEEN(thetaRange(1,it),theta,thetaRange(2,it),-EPSD)) THEN

                p1%dim=line%p1%dim

                EXIT

              ENDIF

            ENDDO

          ENDIF

        ENDIF

        IF(ZERO <= t2 .AND. t2 <= ONE) THEN

          p2=line%p1

          p2%coord(1)=p2%coord(1)+u(1)*t2

          p2%coord(2)=p2%coord(2)+u(2)*t2

          IF(checkTheta) THEN

            ! if the point is on the circle, set %dim to be 2. Otherwise set dummy value.

            theta=ATAN2PI((p2%coord(1)-circle%c%coord(1)),(p2%coord(2)-circle%c%coord(2)))

            p2%dim=-4

            DO it=1,SIZE(thetaRange,DIM=2)

              IF(ISBETWEEN(thetaRange(1,it),theta,thetaRange(2,it),-EPSD)) THEN

                p2%dim=line%p1%dim

                EXIT

              ENDIF

            ENDDO

          ENDIF

        ENDIF

      ENDIF

    ENDIF

    ! if test line has intersections with radial lines

    IF(ALLOCATED(rp) .AND. .NOT.((p1%dim == 2) .AND. (p2%dim == 2))) THEN

      ! there should be maximum two points intersected with the target line

      ! and this circle geom.

      IF(p1%dim /= 2) THEN

        p1=rp(1)

      ELSEIF(p2%dim /= 2) THEN

        p2=rp(1)

      ENDIF

      IF(rp(2)%dim == 2) THEN

        IF(p1%dim /= 2) THEN

          p1=rp(2)

        ELSEIF(p2%dim /= 2) THEN

          p2=rp(2)

        ENDIF

      ENDIF

      CALL rp(1)%clear()

      CALL rp(2)%clear()

      DEALLOCATE(rp)

    ENDIF

  ENDIF

ENDSUBROUTINE intersect_CircleType_and_LineType
@@ -978,24 +1084,16 @@ ELEMENTAL FUNCTION inside_CircleType(circle,point) RESULT(bool) 
  CLASS(CircleType),INTENT(IN) :: circle

  TYPE(PointType),INTENT(IN) :: point

  LOGICAL(SBK) :: bool

  REAL(SRK) :: x,y,dtheta ,cosstt,cosstp,sinstt,sinstp

  REAL(SRK) :: x,y,dtheta

  bool=.FALSE.

  x=point%coord(1)-circle%c%coord(1)

  y=point%coord(2)-circle%c%coord(2)

  dtheta=circle%thetastp-circle%thetastt

  IF((x*x+y*y) .APPROXLE. circle%r*circle%r) THEN

    IF(dtheta .APPROXLE. PI)THEN

      cosstt=COS(circle%thetastt)

      cosstp=COS(circle%thetastp)

      sinstt=SIN(circle%thetastt)

      sinstp=SIN(circle%thetastp)

      IF(COS(circle%thetastt)*y .APPROXGE. SIN(circle%thetastt)*x) &

          bool=COS(circle%thetastp)*y .APPROXLE. SIN(circle%thetastp)*x

    ELSEIF(dtheta < TWOPI) THEN

      cosstt=COS(circle%thetastt)

      cosstp=COS(circle%thetastp)

      sinstt=SIN(circle%thetastt)

      sinstp=SIN(circle%thetastp)

      bool=(COS(circle%thetastt)*y .APPROXGE. SIN(circle%thetastt)*x) &

          .OR. (COS(circle%thetastp)*y .APPROXLE. SIN(circle%thetastp)*x)

    ELSE
@@ -1008,25 +1106,62 @@ ENDFUNCTION inside_CircleType 
!> @brief Determines whether a point lies within a circle type.

!> @param circle the circle type to test for intersection

!> @param point the first to query

!> @param checkRad logical to check radius line in case of partial circle

!> @returns bool .TRUE. if the point lies on the surface of the circle type.

!>

ELEMENTAL FUNCTION onSurface_CircleType(circle,point) RESULT(bool)

ELEMENTAL FUNCTION onSurface_CircleType(circle,point,checkRad) RESULT(bool)

  CLASS(CircleType),INTENT(IN) :: circle

  TYPE(PointType),INTENT(IN) :: point

  LOGICAL(SBK) :: bool

  LOGICAL(SBK),INTENT(IN),OPTIONAL :: checkRad

  LOGICAL(SBK) :: bool,checkRadLine

  REAL(SRK) :: x,y,theta,theta_shift

  TYPE(LineType) :: line

  TYPE(PointType) :: ep

  REAL(SRK) :: stt,stp



  IF(PRESENT(checkRad)) THEN

    checkRadLine=checkRad

  ELSE

    checkRadLine=.FALSE.

  ENDIF

  bool=.FALSE.

  IF(point%dim == 2 .AND. .NOT.(circle%r .APPROXEQA. ZERO)) THEN

    x=point%coord(1)-circle%c%coord(1)

    y=point%coord(2)-circle%c%coord(2)

    IF((x*x+y*y) .APPROXEQA. circle%r*circle%r) THEN

    IF(SOFTEQR((x*x+y*y),circle%r*circle%r,EPSD*100._SRK)) THEN

      theta=ATAN2PI(x,y)

      ! potential issue with Circle is there is no 'normalization' process or

      ! post processing for theta. Probably needs to do something like:

      ! - Stp always should be larger than Stt.

      ! - Stt start between 0 - 360 deg.

      stt=circle%thetastt

      stp=circle%thetastp

      ! correct angStt to start in [0,360.0)

      stt=MODULO(stt,TWOPI)

      stp=MODULO(stp,TWOPI)

      IF(circle%thetastp .APPROXEQ. TWOPI) stp=TWOPI

      theta_shift=0.0_SRK

      IF(circle%thetastt > circle%thetastp) theta_shift=TWOPI

      IF(stt > stp) theta_shift=TWOPI

      IF(y > 0.0_SRK .OR. ((y .APPROXEQA. 0.0_SRK) .AND. &

          (x .APPROXGE. 0.0_SRK))) theta=theta+theta_shift

      bool=(circle%thetastt .APPROXLE. theta) .AND. &

          (theta .APPROXLE. circle%thetastp+theta_shift)

      bool=(stt .APPROXLE. theta) .AND. &

          (theta .APPROXLE. stp+theta_shift)

    ENDIF

    IF(.NOT.bool .AND. checkRadLine) THEN

      IF(ISBETWEEN(ZERO,(circle%thetastp-circle%thetastt),TWOPI,EPSD)) THEN

        CALL ep%init(DIM=2,X=0.0_SRK,Y=0.0_SRK)

        ep%coord(1)=circle%c%coord(1)+ABS(circle%r)*COS(circle%thetastt)

        ep%coord(2)=circle%c%coord(2)+ABS(circle%r)*SIN(circle%thetastt)

        CALL line%set(circle%c,ep)

        bool= ((line%distance2Point(point)) < EPSD*1000_SRK)

        IF(.NOT.bool) THEN

          ep%coord(1)=circle%c%coord(1)+ABS(circle%r)*COS(circle%thetastp)

          ep%coord(2)=circle%c%coord(2)+ABS(circle%r)*SIN(circle%thetastp)

          CALL line%set(circle%c,ep)

          bool= ((line%distance2Point(point)) < EPSD*1000_SRK)

        ENDIF

        CALL ep%clear()

        CALL line%clear()

      ENDIF

    ENDIF

  ENDIF

ENDFUNCTION onSurface_CircleType

src/Geom_Graph.f90

+2 −2
Original line number Diff line number Diff line
@@ -747,7 +747,7 @@ RECURSIVE SUBROUTINE defineQuadraticEdge_coords(this,coord1,coord2,c0,r) 
    y2=coord2(2)-c0(2)

    r2=x2*x2+y2*y2

    rsq=r*r

    IF((rsq .APPROXEQA. r1) .AND. (rsq .APPROXEQA. r2)) THEN

    IF(SOFTEQR(rsq,r1,EPSD*100._SRK) .AND. SOFTEQR(rsq,r2,EPSD*100._SRK)) THEN

      v1=this%getVertIndex(coord1)

      v2=this%getVertIndex(coord2)

      IF(v1 > 0 .AND. v2 > 0 .AND. v1 /= v2) THEN
@@ -773,7 +773,7 @@ RECURSIVE SUBROUTINE defineQuadraticEdge_coords(this,coord1,coord2,c0,r) 
        r2=(y2-y1)

        r2=r2*r2

        d=SQRT(r1+r2)

        IF(d .APPROXEQA. 2.0_SRK*ABS(r)) &

        IF(SOFTEQR(d,2.0_SRK*ABS(r),EPSD*100._SRK)) &

            this%quadEdges(3,v1,v2)=r !sign of r indicates which half

                                           !of semi-circle, all other cases

                                           !traverse shorter arc between points

src/Geom_Line.f90

+4 −0
Original line number Diff line number Diff line
@@ -425,6 +425,8 @@ ELEMENTAL FUNCTION distance2D_to_point(a,b,c) RESULT(d2) 
    ELSE

      !c projects onto ab

      d2=(ac(1)*ac(1)+ac(2)*ac(2))-e*e/f

      ! prevent negative d2 casued by floating point error

      d2=MAX(0.0_SRK,d2)

    ENDIF

  ENDIF

ENDFUNCTION distance2D_to_point
@@ -465,6 +467,8 @@ ELEMENTAL FUNCTION distance3D_to_point(a,b,c) RESULT(d2) 
    ELSE

      !c projects onto ab

      d2=(ac(1)*ac(1)+ac(2)*ac(2)+ac(3)*ac(3))-e*e/f

      ! prevent negative d2 casued by floating point error

      d2=MAX(0.0_SRK,d2)

    ENDIF

  ENDIF

ENDFUNCTION distance3D_to_point

src/Geom_Points.f90

+31 −0
Original line number Diff line number Diff line
@@ -36,6 +36,7 @@ PUBLIC :: OPERATOR(==) 
PUBLIC :: OPERATOR(/=)

PUBLIC :: OPERATOR(.APPROXEQA.)

PUBLIC :: ASSIGNMENT(=)

PUBLIC :: SOFTEQ



INTEGER(SIK),PARAMETER :: MAX_COORD_STR_LEN=128
@@ -174,6 +175,15 @@ INTERFACE ASSIGNMENT(=) 
  !> @copydetails GeomPoints::assign_PointType

  MODULE PROCEDURE assign_PointType

ENDINTERFACE



!> @brief Generic interface for SOFTEQ

!>

!> Compares two points with SOFTEQ

INTERFACE SOFTEQ

  !> @copybrief GeomPoints::SoftEQ_PointType

  !> @copydetails GeomPoints::SoftEQ_PointType

  MODULE PROCEDURE SoftEQ_PointType

ENDINTERFACE SOFTEQ

!

!===============================================================================

CONTAINS
@@ -430,6 +440,27 @@ ELEMENTAL FUNCTION coord2str_PointType(p,fmt_string) RESULT(cstring) 
ENDFUNCTION coord2str_PointType

!

!-------------------------------------------------------------------------------

!> @brief a wrapper to perform SOFTEQ of two points

!> @param p0 first testing point

!> @param p1 second testing point

!> @param tol tolerance parameter for SOFTEQ

!> @returns @c bool the boolean result of the operation

!>

ELEMENTAL FUNCTION SoftEQ_PointType(p0,p1,tol) RESULT(bool)

  CLASS(PointType),INTENT(IN) :: p0,p1

  REAL(SRK),INTENT(IN) :: tol

  LOGICAL(SBK) :: bool

  INTEGER(SIK) :: i

  bool=.FALSE.

  IF((p0%dim == p1%dim) .AND. (p0%dim > 0)) THEN

    DO i=1,p0%dim

      bool=SOFTEQ(p0%coord(i),p1%coord(i),tol)

      IF(.NOT.bool) EXIT

    ENDDO

  ENDIF

ENDFUNCTION SoftEQ_PointType

!

!-------------------------------------------------------------------------------

!> @brief Computes the distance between two points

!> @param p0 the first point

!> @param p1 the second point