diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index 8e33672f615..d6447db6638 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -8,7 +8,7 @@ module MethodCellTernaryModule use CellDefnModule use SubcellTriModule, only: SubcellTriType, create_subcell_tri use ParticleModule - use GeomUtilModule, only: area + use GeomUtilModule, only: area, transform use ConstantsModule, only: DZERO, DONE, DTWO, DSIX use GeomUtilModule, only: point_in_polygon implicit none @@ -158,6 +158,8 @@ subroutine apply_mct(this, particle, tmax) real(DP), intent(in) :: tmax ! local integer(I4B) :: i + real(DP) :: x, y, z, xO, yO + real(DP), allocatable :: xs(:), ys(:) ! (Re)allocate type-bound arrays select type (cell => this%cell) @@ -194,10 +196,24 @@ subroutine apply_mct(this, particle, tmax) allocate (this%xvertnext(this%nverts)) allocate (this%yvertnext(this%nverts)) + ! New origin is the corner of the cell's + ! bounding box closest to the old origin + allocate (xs(this%nverts)) + allocate (ys(this%nverts)) + xs = cell%defn%polyvert(1, :) + ys = cell%defn%polyvert(2, :) + xO = xs(minloc(abs(xs), dim=1)) + yO = ys(minloc(abs(ys), dim=1)) + deallocate (xs) + deallocate (ys) + ! Cell vertices do i = 1, this%nverts - this%xvert(i) = cell%defn%polyvert(1, i) - this%yvert(i) = cell%defn%polyvert(2, i) + x = cell%defn%polyvert(1, i) + y = cell%defn%polyvert(2, i) + call transform(x, y, DZERO, x, y, z, xO, yO) + this%xvert(i) = x + this%yvert(i) = y end do ! Top, bottom, and thickness @@ -212,13 +228,21 @@ subroutine apply_mct(this, particle, tmax) this%iprev = cshift(this%iprev, -1) this%xvertnext = cshift(this%xvert, 1) this%yvertnext = cshift(this%yvert, 1) - end select - ! Calculate vertex velocities. - call this%vertvelo() + ! Calculate vertex velocities. + call this%vertvelo() + + ! Transform particle coordinates + call particle%transform(xO, yO) + + ! Track the particle across the cell. + call this%track(particle, 2, tmax) - ! Track the particle across the cell. - call this%track(particle, 2, tmax) + ! Transform particle coordinates back + call particle%transform(xO, yO, invert=.true.) + call particle%reset_transform() + + end select end subroutine apply_mct @@ -403,16 +427,11 @@ subroutine vertvelo(this) sixa = areacell * DSIX wk1 = -(this%xvert * this%yvertnext - this%xvertnext * this%yvert) nvert = size(this%xvert) - if (nvert == 3) then - this%xctr = sum(this%xvert) / 3.0_DP - this%yctr = sum(this%yvert) / 3.0_DP - else - this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa - this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa - end if + this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa + this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa - ! TODO: do we need to check if center is in polygon? - ! only if using the centroid method which misbehaves + ! TODO: can we use some of the terms of the centroid calculation + ! to do a cheap point in polygon check? ! ! allocate(poly(2, nvert)) ! poly(1,:) = this%xvert