Unverified Commit 8d9b1b01 authored by Graham, Aaron's avatar Graham, Aaron Committed by GitHub
Browse files

Hotfix geom graph combine take2 (#302) 

* Mass rename and interface improvements to Geom_* files

* WIP Add failing BWR test to unit test

* WIP first pass at refactoring combine_GraphType

Everything compiles

* WIP All futility unit tests working

* WIP add generic for Geom_Plane

* WIP fix graphs: add theta filtering for circle-circle intersections

* WIP add comments, fix small bugs, delete old combine_graphtype routine

* Fix bugs in adding circular arcs

* Remove unused variable and debug statements

* Change == to .APPROXEQA. in unit test

* Add RECURSIVE attribute to subroutines

* Add more enumerations
parent 3ebc783f

src/Geom_CircCyl.f90

+73 −20
Original line number Diff line number Diff line
@@ -28,6 +28,7 @@ IMPLICIT NONE 
PRIVATE !Default contents of module to private

!

! List of Public items

PUBLIC :: generateCircle

PUBLIC :: CircleType

PUBLIC :: CylinderType

PUBLIC :: OPERATOR(==)
@@ -53,13 +54,14 @@ TYPE :: CircleType 
    PROCEDURE,PASS :: clear => clear_CircleType

    !> @copybrief Geom_CircCyl::intersect_CircleType_and_LineType

    !> @copydetails Geom_CircCyl::intersect_CircleType_and_LineType

    PROCEDURE,PASS :: intersectLine => intersect_CircleType_and_LineType

    PROCEDURE,PASS,PRIVATE :: intersectLine => intersect_CircleType_and_LineType

    !> @copybrief Geom_CircCyl::intersect_CircleType_and_CircleType

    !> @copydetails Geom_CircCyl::intersect_CircleType_and_CircleType

    PROCEDURE,PASS :: intersectCircle => intersect_CircleType_and_CircleType

    PROCEDURE,PASS,PRIVATE :: intersectCircle => intersect_CircleType_and_CircleType

    !> @copybrief Geom_CircCyl::intersect_ArcCircleType_and_LineType

    !> @copydetails Geom_CircCyl::intersect_ArcCircleType_and_LineType

    PROCEDURE,PASS :: intersectArcLine => intersect_ArcCircleType_and_LineType

    PROCEDURE,PASS,PRIVATE :: intersectArcLine => intersect_ArcCircleType_and_LineType

    GENERIC :: intersect => intersectLine,intersectCircle,intersectArcLine

    !> @copybrief Geom_CircCyl::inside_CircleType

    !> @copydetails Geom_CircCyl::inside_CircleType

    PROCEDURE,PASS :: inside => inside_CircleType
@@ -90,7 +92,8 @@ TYPE :: CylinderType 
    PROCEDURE,PASS :: clear => clear_CylinderType

    !> @copybrief Geom_CircCyl::intersect_CylinderType_and_LineType

    !> @copydetails Geom_CircCyl::intersect_CylinderType_and_LineType

    PROCEDURE,PASS :: intersectLine => intersect_CylinderType_and_LineType

    PROCEDURE,PASS,PRIVATE :: intersectLine => intersect_CylinderType_and_LineType

    GENERIC :: intersect => intersectLine

ENDTYPE CylinderType



!> @brief Generic interface for 'is equal to' operator (==)
@@ -109,6 +112,56 @@ ENDINTERFACE 
CONTAINS

!

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

!> @brief Generates a circle object containing an arc between 2 points

!> @param p1 the first point of the arc

!> @param p2 the second point of the arc

!> @param centroid the centroid of the circle/arc

!> @param radius the radius of the arc

!> @returns circle the resulting @c CircleType object

!>

!> The shortest arc from @c p1 to @c p2 is taken to generate the object.  In the

!> event that a semi-circle is being represented, the arc will start at @c p1 and

!> go in the positive azimuthal direction to @c p2 if @c radius is negative;

!> otherwise it will begin at @c p2 and go in the positive azimuthal direction to @c p1

!>

FUNCTION generateCircle(p1,p2,centroid,radius) RESULT(circle)

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

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

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

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

  TYPE(CircleType) :: circle

  !

  REAL(SRK) :: thetastt,thetastp



  thetastt=ATAN2PI(p1%coord(1)-centroid%coord(1),p1%coord(2)-centroid%coord(2))

  thetastp=ATAN2PI(p2%coord(1)-centroid%coord(1),p2%coord(2)-centroid%coord(2))



  !Handle semi-circles first: in this case, the sign of the radius

  !apparently tells us which semi-circle is actually intended

  IF(ABS(thetastt-thetastp) .APPROXEQ. PI) THEN

    CALL circle%set(centroid,radius,thetastt,thetastp)

  !Now handle other arcs and ensure we're taking the shorter of 2 possibilities

  ELSE

    !Distance between thetastt and thetastp is < PI because point 1

    !is on the x-axis so thetastt is either 0 or PI (and we want the shortest arc)

    IF(p1%coord(2) .APPROXEQ. centroid%coord(2)) THEN

      !Point 2 is above the x-axis

      IF(p2%coord(2) >= p1%coord(2)) THEN

        CALL circle%set(centroid,ABS(radius),thetastp,thetastt)

      !Point 2 is below the x-axis

      ELSE

        CALL circle%set(centroid,ABS(radius),thetastt,thetastp)

      ENDIF

    ELSEIF(p1%coord(2) < centroid%coord(2)) THEN

      CALL circle%set(centroid,ABS(radius),thetastt,thetastp)

    ELSE

      CALL circle%set(centroid,radius,thetastp,thetastt)

    ENDIF

  ENDIF



ENDFUNCTION generateCircle

!

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

!> @brief Sets the values for the CircleType object attributes

!> @param circ the CircleType object to set

!> @param p a point for the center of rotation (must be 2-D)
@@ -121,7 +174,7 @@ ELEMENTAL SUBROUTINE set_CircleType(circ,p,r,angstt,angstp) 
  REAL(SRK),INTENT(IN),OPTIONAL :: angstp

  REAL(SRK) :: dtheta

  CALL circ%clear()

  IF(p%dim == 2 .AND. r > ZERO) THEN

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

    circ%c=p

    circ%r=r

    circ%thetastt=ZERO
@@ -167,7 +220,7 @@ ELEMENTAL SUBROUTINE set_CylinderType(cyl,p,q,r,angstt,angstp) 
  REAL(SRK),INTENT(IN),OPTIONAL :: angstp

  REAL(SRK) :: dtheta

  CALL cyl%clear()

  IF(p%dim == 3 .AND. .NOT.(p .APPROXEQA. q)  .AND. r > ZERO) THEN

  IF(p%dim == 3 .AND. .NOT.(p .APPROXEQA. q)  .AND. .NOT.(r .APPROXEQA. ZERO)) THEN

    CALL cyl%axis%set(p,q)

    cyl%r=r

    cyl%thetastt=ZERO
@@ -218,7 +271,7 @@ ELEMENTAL SUBROUTINE intersect_CircleType_and_LineType(circle,line,p1,p2) 
  p1%dim=-1

  p2%dim=-1

  IF(circle%c%dim == 2 .AND. line%p1%dim == 2 .AND. &

      line%p2%dim == 2 .AND. circle%r > 0.0_SRK) THEN

      line%p2%dim == 2 .AND. .NOT.(circle%r .APPROXEQA. ZERO)) THEN

    u(1)=line%p2%coord(1)-line%p1%coord(1)  !dx

    u(2)=line%p2%coord(2)-line%p1%coord(2)  !dy

    w(1)=line%p1%coord(1)-circle%c%coord(1)
@@ -258,12 +311,12 @@ ELEMENTAL SUBROUTINE intersect_CircleType_and_LineType(circle,line,p1,p2) 
        discr=SQRT(discr)

        t1=(-b-discr)*ra

        t2=(-b+discr)*ra

        IF(ZERO < t1 .AND. t1 < ONE) THEN

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

          p1=line%p1

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

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

        ENDIF

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

        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
@@ -300,9 +353,9 @@ ELEMENTAL SUBROUTINE intersect_CircleType_and_CircleType(c1,c2,p1,p2) 
  p2%dim=-1



  IF(c1%c%dim==2 .AND. c2%c%dim==2 .AND. &

      c1%r > 0.0_SRK .AND. c2%r > 0.0_SRK) THEN

      .NOT.(c1%r .APPROXEQA. ZERO) .AND. .NOT.(c2%r .APPROXEQA. ZERO)) THEN

    d=Distance(c1%c,c2%c)

    s=c1%r+c2%r

    s=ABS(c1%r)+ABS(c2%r)

    IF(d > s) THEN

      p1%dim=-2

      p2%dim=-2
@@ -313,8 +366,8 @@ ELEMENTAL SUBROUTINE intersect_CircleType_and_CircleType(c1,c2,p1,p2) 
      v1%coord(1)=v1%coord(1)/m

      v1%coord(2)=v1%coord(2)/m

      CALL p1%init(DIM=2,X=ZERO,Y=ZERO)

      p1%coord(1)=c1%c%coord(1)+v1%coord(1)*c1%r

      p1%coord(2)=c1%c%coord(2)+v1%coord(2)*c1%r

      p1%coord(1)=c1%c%coord(1)+v1%coord(1)*ABS(c1%r)

      p1%coord(2)=c1%c%coord(2)+v1%coord(2)*ABS(c1%r)

      p1%dim=-3

      p2=p1

    ELSE
@@ -405,7 +458,7 @@ ELEMENTAL SUBROUTINE intersect_ArcCircleType_and_LineType(circle,line,p1,p2,p3,p 
  p3%dim=-1

  p4%dim=-1

  IF(circle%c%dim == 2 .AND. line%p1%dim == 2 .AND. &

      line%p2%dim == 2 .AND. circle%r > 0.0_SRK) THEN

      line%p2%dim == 2 .AND. .NOT.(circle%r .APPROXEQA. ZERO)) THEN

    u(1)=line%p2%coord(1)-line%p1%coord(1)  !dx

    u(2)=line%p2%coord(2)-line%p1%coord(2)  !dy

    w(1)=line%p1%coord(1)-circle%c%coord(1)
@@ -470,10 +523,10 @@ ELEMENTAL SUBROUTINE intersect_ArcCircleType_and_LineType(circle,line,p1,p2,p3,p 
        !and ends of arc

        arcLine%p1=circle%c

        c0=circle%c%coord

        c1(1)=c0(1)+circle%r*COS(circle%thetastt)

        c1(2)=c0(2)+circle%r*SIN(circle%thetastt)

        c2(1)=c0(1)+circle%r*COS(circle%thetastp)

        c2(2)=c0(2)+circle%r*SIN(circle%thetastp)

        c1(1)=c0(1)+ABS(circle%r)*COS(circle%thetastt)

        c1(2)=c0(2)+ABS(circle%r)*SIN(circle%thetastt)

        c2(1)=c0(1)+ABS(circle%r)*COS(circle%thetastp)

        c2(2)=c0(2)+ABS(circle%r)*SIN(circle%thetastp)

        CALL arcLine%p2%init(DIM=2,X=c1(1),Y=c1(2))

        p3=arcLine%intersect(line)

        arcLine%p2%coord=c2
@@ -585,7 +638,7 @@ ELEMENTAL SUBROUTINE intersect_CylinderType_and_LineType(cyl,line,p1,p2) 
  p2%dim=-1

  !Check for valid input

  IF(cyl%axis%p1%dim == 3 .AND. cyl%axis%p2%dim == 3 .AND. &

      line%p1%dim == 3 .AND. line%p2%dim == 3 .AND. cyl%r > zero) THEN

      line%p1%dim == 3 .AND. line%p2%dim == 3 .AND. .NOT.(cyl%r .APPROXEQA. ZERO)) THEN

    p1%dim=0

    p2%dim=0
@@ -796,7 +849,7 @@ ELEMENTAL FUNCTION onSurface_CircleType(circle,point) RESULT(bool) 
  LOGICAL(SBK) :: bool

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

  bool=.FALSE.

  IF(point%dim == 2 .AND. circle%r > 0.0_SRK) THEN

  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

src/Geom_Graph.f90

+1411 −968

File changed.

Preview size limit exceeded, changes collapsed.

src/Geom_Line.f90

+2 −2
Original line number Diff line number Diff line
@@ -53,8 +53,8 @@ TYPE :: LineType 
    PROCEDURE,PASS :: midPoint => midPoint_LineType

    !> @copybrief GeomLine::intersect_LineType_and_LineType

    !> @copydetails GeomLine::intersect_LineType_and_LineType

    PROCEDURE,PASS :: intersect => intersect_LineType_and_LineType

    PROCEDURE,PASS :: intersectLine => intersect_LineType_and_LineType

    PROCEDURE,PASS,PRIVATE :: intersectLine => intersect_LineType_and_LineType

    GENERIC :: intersect => intersectLine

    !> @copybrief GeomLine::distance_LineType_to_PointType

    !> @copydetails GeomLine::distance_LineType_to_PointType

    PROCEDURE,PASS :: distance2Point => distance_LineType_to_PointType

src/Geom_Plane.f90

+2 −1
Original line number Diff line number Diff line
@@ -44,7 +44,8 @@ TYPE :: PlaneType 
    PROCEDURE,PASS :: clear => clear_PlaneType

    !> @copybrief GeomPlane::intersect_PlaneType_and_LineType

    !> @copydetails GeomPlane::intersect_PlaneType_and_LineType

    PROCEDURE,PASS :: intersectLine => intersect_PlaneType_and_LineType

    PROCEDURE,PASS,PRIVATE :: intersectLine => intersect_PlaneType_and_LineType

    GENERIC :: intersect => intersectLine

ENDTYPE PlaneType



!> @brief Generic interface for 'is equal to' operator (==)

src/Geom_Poly.f90

+40 −39
Original line number Diff line number Diff line
@@ -111,6 +111,7 @@ TYPE :: PolygonType 
    !> @copybrief Geom_Poly::intersectPoly_PolygonType

    !> @copydetails Geom_Poly::intersectPoly_PolygonType

    PROCEDURE,PASS :: intersectPoly => intersectPoly_PolygonType

    GENERIC :: intersect => intersectLine,intersectPoly

    !> @copybrief Geom_Poly::getRadius_PolygonType

    !> @copydetails Geom_Poly::getRadius_PolygonType

    PROCEDURE,PASS :: getRadius => getRadius_PolygonType
@@ -592,7 +593,7 @@ ENDFUNCTION getOuterRadius_PolygonType 
!> @param bool The logical result of this operation.  TRUE if the point is on

!>        the surface.

!>

PURE RECURSIVE FUNCTION onSurface_PolygonType(thisPoly,point,incSubReg) RESULT(bool)

RECURSIVE FUNCTION onSurface_PolygonType(thisPoly,point,incSubReg) RESULT(bool)

  CLASS(PolygonType),INTENT(IN) :: thisPoly

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

  LOGICAL(SBK),INTENT(IN),OPTIONAL :: incSubReg
@@ -690,7 +691,7 @@ ENDFUNCTION onSurface_PolygonType 
!> @param isSub Optional logical of whether or not it is a recursive subregion check.

!> @param bool The logical result of this operation.  TRUE if the point is inside.

!>

PURE RECURSIVE FUNCTION point_inside_PolygonType(thisPoly,point,isSub) RESULT(bool)

RECURSIVE FUNCTION point_inside_PolygonType(thisPoly,point,isSub) RESULT(bool)

  CLASS(PolygonType),INTENT(IN) :: thisPoly

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

  LOGICAL(SBK),INTENT(IN),OPTIONAL :: isSub
@@ -844,12 +845,12 @@ FUNCTION polygon_inside_PolygonType(thisPoly,thatPoly) RESULT(bool) 
            IF(ALL(thatPoly%quad2edge /= j)) THEN

              CALL tmpLine%set(thatPoly%vert(thatPoly%edge(1,j)), &

                  thatPoly%vert(thatPoly%edge(2,j)))

              p1=Lines(i)%intersectLine(tmpLine)

              p1=Lines(i)%intersect(tmpLine)

            ENDIF

          ELSE

            CALL tmpLine%set(thatPoly%vert(thatPoly%edge(1,j)), &

                thatPoly%vert(thatPoly%edge(2,j)))

            p1=Lines(i)%intersectLine(tmpLine)

            p1=Lines(i)%intersect(tmpLine)

          ENDIF

          !Check if there was an intersection

          IF(p1%dim > 0) THEN
@@ -870,7 +871,7 @@ FUNCTION polygon_inside_PolygonType(thisPoly,thatPoly) RESULT(bool) 
        !Check for line-circle intersections

        DO j=1,thatPoly%nQuadEdge

          CALL createArcFromQuad(thatPoly,j,tmpCirc)

          CALL tmpCirc%intersectLine(Lines(i),p1,p2)

          CALL tmpCirc%intersect(Lines(i),p1,p2)



          IF(p1%dim == -3) p1%dim=2 !Include tangent points
@@ -911,14 +912,14 @@ FUNCTION polygon_inside_PolygonType(thisPoly,thatPoly) RESULT(bool) 
              IF(ALL(thatPoly%quad2edge /= i)) THEN

                CALL tmpLine%set(thatPoly%vert(thatPoly%edge(1,j)), &

                    thatPoly%vert(thatPoly%edge(2,j)))

                CALL Circs(i)%intersectLine(tmpLine,p1,p2)

                CALL Circs(i)%intersect(tmpLine,p1,p2)

              ELSE

                CYCLE

              ENDIF

            ELSE

              CALL tmpLine%set(thatPoly%vert(thatPoly%edge(1,j)), &

                  thatPoly%vert(thatPoly%edge(2,j)))

              CALL Circs(i)%intersectLine(tmpLine,p1,p2)

              CALL Circs(i)%intersect(tmpLine,p1,p2)

            ENDIF



            IF(p1%dim == -3) p1%dim=2 !Include tangent points
@@ -952,7 +953,7 @@ FUNCTION polygon_inside_PolygonType(thisPoly,thatPoly) RESULT(bool) 
          !Loop over quadratic edges

          DO j=1,thatPoly%nQuadEdge

            CALL createArcFromQuad(thatPoly,j,tmpCirc)

            CALL tmpCirc%intersectCircle(Circs(i),p1,p2)

            CALL tmpCirc%intersect(Circs(i),p1,p2)



            IF(p1%dim == -3) p1%dim=2 !Include tangent points
@@ -1042,7 +1043,7 @@ ENDFUNCTION polygon_inside_PolygonType 
!> @param line The line type to intersect with the polygon

!> @param bool The logical result of the operation

!>

PURE FUNCTION doesLineIntersect_PolygonType(thisPolygon,line) RESULT(bool)

FUNCTION doesLineIntersect_PolygonType(thisPolygon,line) RESULT(bool)

  CLASS(PolygonType),INTENT(IN) :: thisPolygon

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

  LOGICAL(SBK) :: bool
@@ -1072,7 +1073,7 @@ PURE FUNCTION doesLineIntersect_PolygonType(thisPolygon,line) RESULT(bool) 
  !Intersect circles first if necessary

  ipoint=1

  Quad: DO i=1,thisPolygon%nQuadEdge

    CALL circles(i)%intersectArcLine(line,tmppoints(ipoint), &

    CALL circles(i)%intersect(line,tmppoints(ipoint), &

        tmppoints(ipoint+1),tmppoints(ipoint+2),tmppoints(ipoint+3))

    !Check if the points are on the circle

    DO j=0,3
@@ -1097,7 +1098,7 @@ PURE FUNCTION doesLineIntersect_PolygonType(thisPolygon,line) RESULT(bool) 
      DO i=1,SIZE(lines)

        !Intersect the line

        IF(ALL(i /= thisPolygon%quad2edge)) THEN

          tmppoints(ipoint)=lines(i)%intersectLine(line)

          tmppoints(ipoint)=lines(i)%intersect(line)

          IF(tmppoints(ipoint)%dim > 0) THEN

            bool=.TRUE.

            EXIT
@@ -1107,7 +1108,7 @@ PURE FUNCTION doesLineIntersect_PolygonType(thisPolygon,line) RESULT(bool) 
      ENDDO

    ELSE

      DO i=1,SIZE(lines)

        tmppoints(ipoint)=lines(i)%intersectLine(line)

        tmppoints(ipoint)=lines(i)%intersect(line)

        IF(tmppoints(ipoint)%dim > 0) THEN

          bool=.TRUE.

          EXIT
@@ -1152,7 +1153,7 @@ ENDFUNCTION doesLineIntersect_PolygonType 
!> @param thatPoly The second polygon type to intersect with the first

!> @param bool The logical result of the operation

!>

PURE FUNCTION doesPolyIntersect_PolygonType(thisPoly,thatPoly) RESULT(bool)

FUNCTION doesPolyIntersect_PolygonType(thisPoly,thatPoly) RESULT(bool)

  CLASS(PolygonType),INTENT(IN) :: thisPoly

  TYPE(PolygonType),INTENT(IN) :: thatPoly

  LOGICAL(SBK) :: bool
@@ -1198,12 +1199,12 @@ PURE FUNCTION doesPolyIntersect_PolygonType(thisPoly,thatPoly) RESULT(bool) 
          IF(ALL(thatPoly%quad2edge /= j)) THEN

            CALL tmpLine%set(thatPoly%vert(thatPoly%edge(1,j)), &

                thatPoly%vert(thatPoly%edge(2,j)))

            p1=Lines(i)%intersectLine(tmpLine)

            p1=Lines(i)%intersect(tmpLine)

          ENDIF

        ELSE

          CALL tmpLine%set(thatPoly%vert(thatPoly%edge(1,j)), &

              thatPoly%vert(thatPoly%edge(2,j)))

          p1=Lines(i)%intersectLine(tmpLine)

          p1=Lines(i)%intersect(tmpLine)

        ENDIF

        !Check if there was an intersection

        IF(p1%dim > 0) THEN
@@ -1224,7 +1225,7 @@ PURE FUNCTION doesPolyIntersect_PolygonType(thisPoly,thatPoly) RESULT(bool) 
      !Check for line-circle intersections

      DO j=1,thatPoly%nQuadEdge

        CALL createArcFromQuad(thatPoly,j,tmpCirc)

        CALL tmpCirc%intersectLine(Lines(i),p1,p2)

        CALL tmpCirc%intersect(Lines(i),p1,p2)



        IF(p1%dim == -3) p1%dim=2 !Include tangent points
@@ -1260,14 +1261,14 @@ PURE FUNCTION doesPolyIntersect_PolygonType(thisPoly,thatPoly) RESULT(bool) 
            IF(ALL(thatPoly%quad2edge /= j)) THEN

              CALL tmpLine%set(thatPoly%vert(thatPoly%edge(1,j)), &

                  thatPoly%vert(thatPoly%edge(2,j)))

              CALL Circs(i)%intersectLine(tmpLine,p1,p2)

              CALL Circs(i)%intersect(tmpLine,p1,p2)

            ELSE

              CYCLE

            ENDIF

          ELSE

            CALL tmpLine%set(thatPoly%vert(thatPoly%edge(1,j)), &

                thatPoly%vert(thatPoly%edge(2,j)))

            CALL Circs(i)%intersectLine(tmpLine,p1,p2)

            CALL Circs(i)%intersect(tmpLine,p1,p2)

          ENDIF
@@ -1297,7 +1298,7 @@ PURE FUNCTION doesPolyIntersect_PolygonType(thisPoly,thatPoly) RESULT(bool) 
        !Loop over quadratic edges

        DO j=1,thatPoly%nQuadEdge

          CALL createArcFromQuad(thatPoly,j,tmpCirc)

          CALL tmpCirc%intersectCircle(Circs(i),p1,p2)

          CALL tmpCirc%intersect(Circs(i),p1,p2)



          IF(p1%dim == -3) p1%dim=2 !Include tangent points
@@ -1357,7 +1358,7 @@ ENDFUNCTION doesPolyIntersect_PolygonType 
!> @brief

!> @param

!>

PURE SUBROUTINE intersectLine_PolygonType(thisPolygon,line,points)

SUBROUTINE intersectLine_PolygonType(thisPolygon,line,points)

  CLASS(PolygonType),INTENT(IN) :: thisPolygon

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

  TYPE(PointType),ALLOCATABLE,INTENT(INOUT) :: points(:)
@@ -1409,7 +1410,7 @@ PURE SUBROUTINE intersectLine_PolygonType(thisPolygon,line,points) 
    npoints=0

    IF(nlines > 0) THEN

      DO i=1,nlines

        p1=lines(i)%intersectLine(line)

        p1=lines(i)%intersect(line)

        IF(p1%dim == 2) THEN

          npoints=npoints+1

          tmpPoints(npoints)=p1
@@ -1426,7 +1427,7 @@ PURE SUBROUTINE intersectLine_PolygonType(thisPolygon,line,points) 
    !Test against arcs

    IF(narcs > 0) THEN

      DO i=1,narcs

        CALL circles(i)%intersectLine(line,p1,p2)

        CALL circles(i)%intersect(line,p1,p2)

        IF(p1%dim == -3) p1%dim=2 !Include tangent points



        !If line segment end points are on circle, intersections
@@ -1501,7 +1502,7 @@ ENDSUBROUTINE intersectLine_PolygonType 
!> @brief

!> @param

!>

PURE SUBROUTINE intersectPoly_PolygonType(thisPoly,thatPoly,points)

SUBROUTINE intersectPoly_PolygonType(thisPoly,thatPoly,points)

  CLASS(PolygonType),INTENT(IN) :: thisPoly

  TYPE(PolygonType),INTENT(IN) :: thatPoly

  TYPE(PointType),ALLOCATABLE,INTENT(INOUT) :: points(:)
@@ -1557,7 +1558,7 @@ PURE SUBROUTINE intersectPoly_PolygonType(thisPoly,thatPoly,points) 
      ENDDO

      !All lines with the those circles

      DO j=1,thatPoly%nQuadEdge

        CALL thoseCircs(j)%intersectArcLine(theseLines(i),tmppoints(ipoint), &

        CALL thoseCircs(j)%intersect(theseLines(i),tmppoints(ipoint), &

            tmppoints(ipoint+1),tmppoints(ipoint+2),tmppoints(ipoint+3))

        !Check if the points are on the circle

        DO k=0,3
@@ -1576,7 +1577,7 @@ PURE SUBROUTINE intersectPoly_PolygonType(thisPoly,thatPoly,points) 
    !TheseCircs

    DO i=1,thisPoly%nQuadEdge

      DO j=1,thatPoly%nVert

        CALL theseCircs(i)%intersectArcLine(thoseLines(j),tmppoints(ipoint), &

        CALL theseCircs(i)%intersect(thoseLines(j),tmppoints(ipoint), &

            tmppoints(ipoint+1),tmppoints(ipoint+2),tmppoints(ipoint+3))

        !Check if the points are on the circle

        DO k=0,3
@@ -1592,7 +1593,7 @@ PURE SUBROUTINE intersectPoly_PolygonType(thisPoly,thatPoly,points) 
      ENDDO

      !ThoseCircs

      DO j=1,thatPoly%nQuadEdge

        CALL theseCircs(i)%intersectCircle(thoseCircs(j),tmppoints(ipoint), &

        CALL theseCircs(i)%intersect(thoseCircs(j),tmppoints(ipoint), &

            tmppoints(ipoint+1))

        ipoint=ipoint+2

      ENDDO
@@ -1765,7 +1766,7 @@ RECURSIVE SUBROUTINE generateGraph_PolygonType(thisPoly,g,incSubReg) 
      r=thisPoly%quadEdge(3,i)

      v0=thisPoly%edge(1,j)

      v1=thisPoly%edge(2,j)

      CALL g%defineQuadraticEdge(verts(:,v0),verts(:,v1),c,r)

      CALL g%defineEdge(verts(:,v0),verts(:,v1),c,r)

    ENDDO

    DO i=1,thisPoly%nVert

      IF(.NOT.isQuadEdge(i)) THEN
@@ -1816,10 +1817,10 @@ SUBROUTINE Polygonize_Circle(circle,polygon) 
      CALL g%insertVertex(v1)

      CALL g%insertVertex(v2)

      CALL g%insertVertex(v3)

      CALL g%defineQuadraticEdge(v0,v1,c,r)

      CALL g%defineQuadraticEdge(v1,v2,c,r)

      CALL g%defineQuadraticEdge(v2,v3,c,r)

      CALL g%defineQuadraticEdge(v0,v3,c,r)

      CALL g%defineEdge(v0,v1,c,r)

      CALL g%defineEdge(v1,v2,c,r)

      CALL g%defineEdge(v2,v3,c,r)

      CALL g%defineEdge(v0,v3,c,r)

    ELSE

      v0(1)=c(1); v0(2)=c(2)

      v1(1)=c(1)+r*COS(circle%thetastt)
@@ -1832,7 +1833,7 @@ SUBROUTINE Polygonize_Circle(circle,polygon) 
      CALL g%insertVertex(v2)

      CALL g%defineEdge(v0,v1)

      CALL g%defineEdge(v0,v2)

      CALL g%defineQuadraticEdge(v1,v2,c,r)

      CALL g%defineEdge(v1,v2,c,r)

    ENDIF



    CALL polygon%set(g)
@@ -1868,10 +1869,10 @@ SUBROUTINE Polygonize_Cylinder(cylinder,polygon) 
      CALL g%insertVertex(v1)

      CALL g%insertVertex(v2)

      CALL g%insertVertex(v3)

      CALL g%defineQuadraticEdge(v0,v1,c(1:2),r)

      CALL g%defineQuadraticEdge(v1,v2,c(1:2),r)

      CALL g%defineQuadraticEdge(v2,v3,c(1:2),r)

      CALL g%defineQuadraticEdge(v0,v3,c(1:2),r)

      CALL g%defineEdge(v0,v1,c(1:2),r)

      CALL g%defineEdge(v1,v2,c(1:2),r)

      CALL g%defineEdge(v2,v3,c(1:2),r)

      CALL g%defineEdge(v0,v3,c(1:2),r)

    ELSE

      v0(1)=c(1); v0(2)=c(2)

      v1(1)=c(1)+r*COS(cylinder%thetastt)
@@ -1884,7 +1885,7 @@ SUBROUTINE Polygonize_Cylinder(cylinder,polygon) 
      CALL g%insertVertex(v2)

      CALL g%defineEdge(v0,v1)

      CALL g%defineEdge(v0,v2)

      CALL g%defineQuadraticEdge(v1,v2,c(1:2),r)

      CALL g%defineEdge(v1,v2,c(1:2),r)

    ENDIF



    CALL polygon%set(g)
@@ -1966,7 +1967,7 @@ ENDSUBROUTINE Polygonize_ABBox 
!> @param p2 The second polygon type to compare

!> @param bool The logical result of the operation

!>

PURE RECURSIVE FUNCTION isequal_PolygonType(p1,p2) RESULT(bool)

RECURSIVE FUNCTION isequal_PolygonType(p1,p2) RESULT(bool)

  TYPE(PolygonType),INTENT(IN) :: p1

  TYPE(PolygonType),INTENT(IN) :: p2

  LOGICAL(SBK) :: bool
@@ -2008,7 +2009,7 @@ ENDFUNCTION isequal_PolygonType 
!> @param iquad The desired quadratic edge index to convert to a circle type

!> @param circle The circle type created from the polygon's quadratic edge

!>

PURE SUBROUTINE createArcFromQuad(thisPoly,iquad,circle)

SUBROUTINE createArcFromQuad(thisPoly,iquad,circle)

  CLASS(PolygonType),INTENT(IN) :: thisPoly

  INTEGER(SIK),INTENT(IN) :: iquad

  TYPE(CircleType),INTENT(INOUT) :: circle