[WIP] Polyregions

src/2d/flagregions2.f90

@@ -39,9 +39,10 @@ subroutine flagregions2(mx,my,mbuff,xlower,ylower,dx,dy,level,t, &
     real(kind=8), intent(in) :: DOFLAG
 
     ! Locals
-    integer :: i,j,m,i1,i2,j1,j2
-    real(kind=8) :: x_low,y_low,x_hi,y_hi, xupper,yupper
+    integer :: i,j,m,i1,i2,j1,j2,ip0,intersections
+    real(kind=8) :: x,y,xp0,yp0,xp1,yp1,xupper,yupper
     integer, allocatable :: minlevel(:,:), maxlevel(:,:)
+    logical :: isinside
 
     allocate(minlevel(mx,my), maxlevel(mx,my))
 
@@ -72,8 +73,36 @@ subroutine flagregions2(mx,my,mbuff,xlower,ylower,dx,dy,level,t, &
 
         do j=j1,j2
             do i=i1,i2
-                minlevel(i,j) = max(minlevel(i,j), regions(m)%min_level)
-                maxlevel(i,j) = max(maxlevel(i,j), regions(m)%max_level)
+                ! we are in the bounding box of the arbitrary polygon
+                ! now check point(i,j), may be in bounding box but not polygon
+                x = xlower + (i-0.5d0)*dx
+                y = ylower + (j-0.5d0)*dy
+                ! ray casting
+                isinside = .false. !alternated with each intersection
+                do ip0 = 1, regions(m)%num_points
+                !looping over edges of polygon
+                    xp0 = regions(m)%xpoints(ip0)
+                    yp0 = regions(m)%ypoints(ip0)
+                    xp1 = regions(m)%xpoints(ip0+1)
+                    yp1 = regions(m)%ypoints(ip0+1)
+                    !ray casting from point (x,y), horizontal ray to + infinity
+                    if (yp0==yp1) then !dismiss horizontal intersection of any kind
+                        cycle !ray on top of horizontal segement...do not count as intersection
+                    elseif ((y>max(yp0,yp1)).or.(y<min(yp0,yp1)).or.(x>max(xp0,xp1))) then
+                        cycle !no intersection ...ray above/below or to the right of segment
+                    elseif (x<=min(xp0,xp1)) then !intersection. 
+                        !point is to left of segement. ray must intersect
+                        isinside = .not.isinside !intesection...change isinside
+                    elseif (x<=((xp1-xp0)*(y-yp0)/(yp1-yp0) + xp0)) then
+                        !checked actual intersection...x to the left of segment...must intersect
+                        isinside = .not.isinside !intesection...change isinside
+                    endif
+                enddo
+
+                if (isinside) then
+                    minlevel(i,j) = max(minlevel(i,j), regions(m)%min_level)
+                    maxlevel(i,j) = max(maxlevel(i,j), regions(m)%max_level)
+                endif
             enddo
          enddo
     enddo rloop

src/2d/regions_module.f90

@@ -8,18 +8,26 @@ module regions_module
     implicit none
     save
 
-    ! Region type definition
-    type region_type
+    ! poly Region type definition
+    type polyregion_type
         integer :: min_level,max_level
-        real(kind=8) :: x_low,y_low,x_hi,y_hi,t_low,t_hi
-    end type region_type
+        integer :: num_points
+        real(kind=8) :: t_low,t_hi
+        real(kind=8) :: x_low,x_hi,y_low,y_hi
+        real(kind=8), allocatable :: xpoints(:)
+        real(kind=8), allocatable :: ypoints(:)
+
+    end type polyregion_type
+
 
     integer :: num_regions
-    type(region_type), allocatable :: regions(:)
+    type(polyregion_type), allocatable :: regions(:)
 
 contains
 
+
     subroutine set_regions(fname)
+        ! this sets arbitrary polygon regions
 
         use amr_module
 
@@ -30,7 +38,7 @@ subroutine set_regions(fname)
 
         ! Locals
         integer, parameter :: unit = 7
-        integer :: i
+        integer :: i,ip,np
 
         write(parmunit,*) ' '
         write(parmunit,*) '--------------------------------------------'
@@ -48,17 +56,32 @@ subroutine set_regions(fname)
             write(parmunit,*) '  No regions specified for refinement'
 
         else
-            ! Refinement region data
+            ! Refinement poly region data
             allocate(regions(num_regions))
             do i=1,num_regions
-                read(unit,*) regions(i)%min_level, regions(i)%max_level, &
-                             regions(i)%t_low, regions(i)%t_hi, &
-                             regions(i)%x_low, regions(i)%x_hi, &
-                             regions(i)%y_low, regions(i)%y_hi
+                read(unit,*) regions(i)%num_points, &
+                             regions(i)%min_level, regions(i)%max_level, &
+                             regions(i)%t_low, regions(i)%t_hi
+
+                np = regions(i)%num_points
+                allocate(regions(i)%xpoints(np+1))
+                allocate(regions(i)%ypoints(np+1))
+
+                do ip = 1,np
+                    read(unit,*) regions(i)%xpoints(ip), &
+                                 regions(i)%ypoints(ip) 
+                enddo
+                regions(i)%xpoints(np+1) = regions(i)%xpoints(1) 
+                regions(i)%ypoints(np+1) = regions(i)%ypoints(1) 
+                regions(i)%x_low = minval(regions(i)%xpoints)
+                regions(i)%y_low = minval(regions(i)%ypoints)
+                regions(i)%x_hi =  maxval(regions(i)%xpoints)
+                regions(i)%y_hi =  maxval(regions(i)%ypoints)
             enddo
         endif
         close(unit)
 
     end subroutine set_regions
 
+
 end module regions_module

src/python/amrclaw/data.py

@@ -141,8 +141,26 @@ def write(self,out_file='regions.data',data_source='setrun.py'):
         self.data_write(value=len(self.regions),alt_name='num_regions')
         if (self.num_dim == 3) and (len(self.regions) > 0):
             raise NotImplementedError("*** Regions not yet implemented in 3d")
-        for regions in self.regions:
-            self._out_file.write(8*"%g  " % tuple(regions) +"\n")
+        for region in self.regions:
+            if (len(region)==5):
+                #this is a polyregion specification style
+                #[minlevel,maxlevle,t1,tend,[(x1,y1),(),(),(xpts,ypts)]]
+                pregion = region
+                pts = len(pregion[4])
+            else:
+                #this is a rectangular specification...turn into pregion
+                pts = 4
+                x1 = region[4]
+                x2 = region[5]
+                y1 = region[6]
+                y2 = region[7]
+                ptlist = [(x1,y1),(x2,y1),(x2,y2),(x1,y2)]
+                pregion = region[0:4]+[ptlist]
+            #write out pregion
+            self._out_file.write("%g " % pts)
+            self._out_file.write(4*"%g  " % tuple(pregion[0:4]) +"\n")
+            for pt in xrange(pts):
+                self._out_file.write(2*"%g  " % tuple(pregion[4][pt]) +"\n")
         self.close_data_file()
 
 
@@ -173,10 +191,15 @@ def read(self, path, force=False):
         self.regions = []
         for n in xrange(num_regions):
             line = data_file.readline().split()
-            self.regions.append([int(line[0]), int(line[1]),
-                                 float(line[2]), float(line[3]), 
-                                 float(line[4]), float(line[5]),
-                                 float(line[6]), float(line[7])])
+            pts = int(line[0])
+            pregion1 =  [int(line[1]),int(line[2]), float(line[3]), float(line[4])]
+            pregionpts = []
+            for pt in xrange(pts):
+                xypt =  data_file.readline().split()
+                pregionpts.append((float(xypt[0]),float(xypt[1])))
+
+            pregion = pregion1+[pregionpts]
+            self.regions.append(pregion)
 
         data_file.close()