From 2cdf4ee152d5b2f6a658ae05200e7f0c893e1add Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 01:02:43 -0500
Subject: [PATCH 01/22] Add `bilinearintegral` function

---
 src/python/geoclaw/topotools.py | 41 +++++++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index a8b8db08d..0f477c16d 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1593,6 +1593,47 @@ def make_shoreline_xy(self, sea_level=0):
             shoreline_xy = numpy.vstack((shoreline_xy, \
                            numpy.array([numpy.nan,numpy.nan]), c.allsegs[0][k]))
         return shoreline_xy
+       
+    def bilinearintegral(self, domain, topodomain, corners):
+        r"""
+        Suppose corners are (x1, y1), (x1, y2), (x2, y1), (x1, y2). Get a function,
+        f(xi, eta) = a * (xi - x1) / (x2 - x1) + b * (eta - y1) / (y2 - y1)
+                     + c * (xi - x1) * (eta - y1) / {(x2 - x1) * (y2 - y1)} + d,
+        to represent topo of topomain by bilinear interpolation with four corner
+        points. Then integrate f(xi, eta) on the rectangular domain which is the
+        intersection of domain and topodomain.
+        
+        :Input:
+         - *domain* (ndarray(:))
+         - *topodomain* (ndarray(:))
+         - *corners* (ndarray(:, :)) Four corner points of the topodomain.
+
+        :Output:
+         - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
+           the intersection of domain and topodomain.
+        """
+    
+        # Set boundary for integral
+        bound = self.intersection(domain, topodomain)[2]
+        
+        # Find terms for the integration
+        deltax = topodomain[1] - topodomain[0]
+        deltay = topodomain[3] - topodomain[2]
+        area = (bound[3] - bound[2]) * (bound[1] - bound[0])
+        sumxi = (bound[1] + bound[0] - 2.0 * topodomain[0]) / deltax
+        sumeta = (bound[3] + bound[2] - 2.0 * topodomain[2]) / deltay
+        
+        # Find coefficients of the function, f(xi, eta)
+        a = corners[1][0] - corners[0][0]
+        b = corners[0][1] - corners[0][0]
+        c = corners[1][1] - corners[1][0] - corners[0][1] + corners[0][0]
+        d = corners[0][0]
+        
+        bilinearintegral = (0.5 * (a * sumxi + b * sumeta) +
+                            0.25 * c * sumxi * sumeta + d) * area
+    
+        return bilinearintegral
+   
 
 
 

From c64127e00a1a3df4dcdfbc0d35abafafabdf8258 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 01:05:51 -0500
Subject: [PATCH 02/22] Update `bilinearintegral_s` function

---
 src/python/geoclaw/topotools.py | 52 ++++++++++++++++++++++++++++++++-
 1 file changed, 51 insertions(+), 1 deletion(-)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 0f477c16d..138707caf 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1633,7 +1633,57 @@ def bilinearintegral(self, domain, topodomain, corners):
                             0.25 * c * sumxi * sumeta + d) * area
     
         return bilinearintegral
-   
+
+    def bilinearintegral_s(self, domain, topodomain, corners):
+        r"""
+        Suppose corners are (x1, y1), (x1, y2), (x2, y1), (x1, y2). Get a function,
+        f(theta, phi) = a * (theta - x1) + b * (phi - y1)
+                        + c * (theta - x1) * (phi - y1) + d,
+        where theta is longitute and phi is latitude, to represent the topo of
+        topomain by bilinear interpolation with four corner points. Then integrate
+        f(theta, phi) on the domain which is the intersection of the domain
+        and the topodomain on the sphere.
+        
+        :Input:
+         - *domain* (ndarray(:))
+         - *topodomain* (ndarray(:))
+         - *corners* (ndarray(:, :)) Four corner points of the topodomain.
+
+        :Output:
+         - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
+           the intersection of domain and topodomain.
+        """
+    
+        #Set boundary parameter
+        bound = self.intersection(domain, topodomain)[2]
+
+        # Find terms for the integration
+        xdiffhi = bound[1] - topodomain[0]
+        xdifflow = bound[0] - topodomain[0]
+        ydiffhi = bound[3] - topodomain[2]
+        ydifflow = bound[2] - topodomain[2]
+        xdiff2 = 0.5 * (xdiffhi**2 - xdifflow**2)
+        intdx = bound[1] - bound[0]
+        deltax = topodomain[1] - topodomain[0]
+        deltay = topodomain[3] - topodomain[2]
+        
+        d2r = DEG2RAD; r2d = RAD2DEG
+
+        cbsinint = (r2d * numpy.cos(d2r * bound[3]) + ydiffhi * numpy.sin(d2r * bound[3]))
+        - (r2d * numpy.cos(d2r * bound[2]) + ydifflow * numpy.sin(d2r * bound[2]))
+
+        adsinint = r2d * (numpy.sin(d2r * bound[3]) - numpy.sin(d2r * bound[2]))
+
+        # Find coefficients of the function, f(theta, phi)
+        a = (corners[1][0] - corners[0][0]) / deltax
+        b = (corners[0][1] - corners[0][0]) / deltay
+        c = (corners[1][1] - corners[1][0] - corners[0][1] + corners[0][0]) / (deltax * deltay)
+        d = corners[0][0]
+
+        bilinearintegral_s = ((a * xdiff2 + d * intdx) * adsinint +
+                              r2d * (c * xdiff2 + b * intdx) * cbsinint) * (Rearth * d2r)**2
+
+        return bilinearintegral_s
 
 
 

From 7e5133d35b61c649cae2de3693a124b9c3795c3c Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 11:55:21 -0500
Subject: [PATCH 03/22] Add `topointegral` function

---
 src/python/geoclaw/topotools.py | 101 ++++++++++++++++++++++++++++++++
 1 file changed, 101 insertions(+)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 138707caf..96c830e89 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1685,6 +1685,107 @@ def bilinearintegral_s(self, domain, topodomain, corners):
 
         return bilinearintegral_s
 
+    def topointegral(self, domain, topoparam, z, intmethod, coordinate_system):
+        r"""
+        topointegral integrates a surface over a rectangular region which is
+        the intersection of the domain and topo which is defined by "topoparam".
+        The surface integrated is defined by a function gotten by bilinear
+        interpolation through corners of the Cartesian grid.
+        
+        :Input:
+         - *domain* (ndarray(:))
+         - *topoparam* (ndarray(:)) topoparam includes topo's x and y coordinates'
+           minimum and maximum, and topo's grid size, dx and dy.
+         - *z* (ndarray(:, :)) Topo data of domain represented by topoparam.
+         - *intmethod* (int) The method of integral.
+         - *coordinate_system* (int) The type of coordinate system.
+
+        :Output:
+         - *theintegral* (float) The integral of a surface over a rectangular
+           region which is the intersection of the domain and "topoparam".
+        """
+        
+        theintegral = 0.0
+        
+        # Set topo's parameter
+        x1 = topoparam[0]; x2 = topoparam[1]
+        y1 = topoparam[2]; y2 = topoparam[3]
+        topodx = topoparam[4]; topody = topoparam[5]
+        topomx = numpy.array(z).shape[0]
+        topomy = numpy.array(z).shape[1]
+        
+        # Find the intersection of the domain and topo
+        bound = self.intersection(domain, topoparam)[2]
+        xlow = bound[0]; xhi = bound[1]
+        ylow = bound[2]; yhi = bound[3]
+        
+        if intmethod == 1:
+            
+            # Find topo's start and end points for integral
+            # The x coodinate of the topo's start point
+            istart = -1
+            for i in range(topomx):
+                if (x1 + i * topodx) > xlow:
+                    istart = i
+                    break
+            
+            # The y coodinate of the topo's start point
+            jstart = -1
+            for j in range(topomy):
+                if (y1 + j * topody)  > ylow:
+                    jstart = j
+                    break
+            
+            # The x coodinate of the topo's end point
+            iend = topomx
+            for m in range(topomx):
+                if (x1 + m * topodx) >= xhi:
+                    iend = m
+                    break
+            
+            # The y coodinate of the topo's end point
+            jend = topomy
+            for n in range(topomy):
+                if (y1 + n * topody) >= yhi:
+                    jend = n
+                    break
+                            
+            # Prevent overflow
+            jstart = max(jstart, 1)
+            istart = max(istart, 1)
+            jend = min(jend, topomy - 1)
+            iend = min(iend, topomx - 1)
+            
+            # Topo integral
+            for jj in range(jstart, jend + 1):
+                yint1 = y1 + (jj - 1.0) * topody
+                yint2 = y1 + jj * topody
+
+                for ii in range(istart, iend + 1):
+                    xint1 = x1 + (ii - 1.0) * topodx
+                    xint2 = x1 + ii * topodx
+
+                    # Four corners of topo's grid
+                    corners = numpy.empty((2, 2))
+                    corners[0][0] = z[ii-1][topomy - jj]
+                    corners[0][1] = z[ii-1][topomy - 1 - jj]
+                    corners[1][0] = z[ii][topomy - jj]
+                    corners[1][1] = z[ii][topomy - 1 - jj]
+
+                    # Parameters of topo's grid integral
+                    topointparam = [xint1, xint2, yint1, yint2]
+
+                    if coordinate_system == 1:
+                        theintegral += self.bilinearintegral(domain, topointparam, corners)
+                    elif coordinate_system == 2:
+                        theintegral += self.bilinearintegral_s(domain, topointparam, corners)
+                    else:
+                         print("TOPOINTEGRAL: coordinate_system error")
+
+        else:
+            print("TOPOINTEGRAL: only intmethod = 1,2 is supported")
+
+        return theintegral
 
 
 # Define convenience dictionary of URLs for some online DEMs in netCDF form:

From 76a44655d207a1e87373b80a41f249a997b286a9 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 11:58:03 -0500
Subject: [PATCH 04/22] Add `intersection` function

---
 src/python/geoclaw/topotools.py | 42 +++++++++++++++++++++++++++++++++
 1 file changed, 42 insertions(+)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 96c830e89..ddff60827 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1787,6 +1787,48 @@ def topointegral(self, domain, topoparam, z, intmethod, coordinate_system):
 
         return theintegral
 
+    def intersection(self, rect1, rect2):
+        r"""
+        Whether two domains, rect1 and rect2 intersect. If they intersect,
+        find the intersection's boundary and area.
+        
+        :Input:
+         - *rect1* (ndarray(:))
+         - *rect2* (ndarray(:))
+
+        :Output:
+         - *mark* (bool) If rect1 and rect2 intersect, mark = True. Otherwise,
+           mark = False.
+         - *area* (float) The area of the intersection of rect1 and rect2.
+         - *bound* (ndarray(:)) The boundary of the intersection of rect1 and rect2.
+        """
+        
+        # Set boundary for rect1
+        x1_low = rect1[0]; x1hi = rect1[1]
+        y1_low = rect1[2]; y1hi = rect1[3]
+        
+        # Set boundary for rect2
+        x2low = rect2[0]; x2hi = rect2[1]
+        y2low = rect2[2]; y2hi = rect2[3]
+        
+        # Boundary of the intersection part
+        xintlow = max(x1_low, x2low)
+        xinthi = min(x1hi, x2hi)
+        yintlow = max(y1_low, y2low)
+        yinthi = min(y1hi, y2hi)
+
+        # Whether rect1 and rect2 intersect
+        if xinthi > xintlow and yinthi > yintlow:
+            area = (xinthi - xintlow) * (yinthi - yintlow)
+            mark = True
+            bound = [xintlow, xinthi, yintlow, yinthi]
+        else:
+            area = 0.0
+            mark = False
+            bound = [-1, -1, -1, -1]
+
+        return mark, area, bound
+       
 
 # Define convenience dictionary of URLs for some online DEMs in netCDF form:
 remote_topo_urls = {}

From d88cfba7f6cb497928c0899ba11a2c1ece4c9314 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 12:08:41 -0500
Subject: [PATCH 05/22] Add `recurintegral` function

---
 src/python/geoclaw/topotools.py | 61 +++++++++++++++++++++++++++++++++
 1 file changed, 61 insertions(+)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index ddff60827..080f5ac77 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1828,6 +1828,67 @@ def intersection(self, rect1, rect2):
             bound = [-1, -1, -1, -1]
 
         return mark, area, bound
+
+    def recurintegral(self, domain, mtopoorder, mtopofiles, m, topo):
+        r"""
+        Compute the integral of topo over the rectangle domain define by "domain".
+        Find the index of the topo whose resolution is m by mtopoorder. Use this
+        index to find corrosponding topo object in the list "topo".
+
+        The main call to this recursive function has corners of a grid cell for the
+        rectangle and m = 1 in order to compute the integral over the cell
+        using all topo objects.
+
+        The recursive strategy is to first compute the integral using only the topo
+        object of m+1 resolution. Then apply corrections due to adding m resolution
+        topo object.
+         
+        Corrections are needed if the new topo object intersects the grid cell.
+        Let the intersection be represented by mark2[2]. Two corrections are needed.
+        First, subtract out the integral over the rectangle "mark2[2]" computed using
+        topo objects mtopoorder(mtopofiles) to mtopoorder(m+1), and then adding in
+        the integral over this same region using the topo object mtopoorder(m).
+        
+        :Input:
+         - *domain* (ndarray(:)) The specific surface where the integral is computed.
+         - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+         - *mtopofiles* (int) The number of the topo objects.
+         - *m* (int) Resolution of the topo objects.
+         - *topo* (list) The list of topo objects.
+
+        :Output:
+         - *integral* (float) Integral.
+        """
+        
+        # The index of the topo object whose resolution is m
+        mfid = mtopoorder[m]
+        
+        # The parameters of the m resolution topo object
+        topoparam_mfid = [topo[mfid].extent[0], topo[mfid].extent[1], topo[mfid].extent[2],
+                          topo[mfid].extent[3], topo[mfid].delta[0], topo[mfid].delta[1]]
+        
+        # Innermost step of recursion reaches this point, only using coarsest topo grid
+        # compute directly
+        if m == mtopofiles - 1:
+            mark1 = self.intersection(domain, topoparam_mfid)
+            if mark1[0]:
+                integral = self.topointegral(mark1[2], topoparam_mfid, topo[mfid].Z, 1, 1)
+            else:
+                integral = 0.0
+        else:
+            int1 = self.recurintegral(domain, mtopoorder, mtopofiles, m+1, topo)
+
+            # Whether new topo object intersects the grid cell
+            mark2 = self.intersection(domain, topoparam_mfid)
+            
+            # The new topo object intersects the grid cell. Corrections are needed
+            if mark2[1] > 0:
+                int2 = self.recurintegral(mark2[2], mtopoorder, mtopofiles, m+1, topo)
+                int3 = self.topointegral(mark2[2], topoparam_mfid, topo[mfid].Z, 1, 1)
+                integral = int1 - int2 + int3
+            else:
+                integral = int1
+        return integral
        
 
 # Define convenience dictionary of URLs for some online DEMs in netCDF form:

From 3f8b55bf783ff500cdbce1612b7eb74c030718dd Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 12:15:52 -0500
Subject: [PATCH 06/22] Add `cellintegral` function

---
 src/python/geoclaw/topotools.py | 46 ++++++++++++++++++++++++++++++++-
 1 file changed, 45 insertions(+), 1 deletion(-)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 080f5ac77..00dbe6a73 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1889,8 +1889,52 @@ def recurintegral(self, domain, mtopoorder, mtopofiles, m, topo):
             else:
                 integral = int1
         return integral
-       
 
+    def cellintegral(self, cell, mtopoorder, mtopofiles, topo):
+        r"""
+        Compute the integral on the surface defined by "cell" with the topo objects list
+        "topo".
+        
+        :Input:
+        - *cell* (ndarray(:)) The boundary of specific surface where the integral
+          is computed.
+        - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+        - *mtopofiles* (int) The number of the topo objects.
+        - *topo* (list) The list of topo objects.
+
+        :Output:
+        - *topoint* (float) The integral on the cell, "cell".
+        """
+        
+        # Initialization
+        topoint = 0.0
+
+        for m in range(mtopofiles):
+        
+            # The index of the topo file object resolution is m
+            mfid = mtopoorder[m]
+            cellarea = (cell[1] - cell[0]) * (cell[3] - cell[2])
+            
+            # The parameters of this topo object
+            topoparam_mfid = [topo[mfid].extent[0], topo[mfid].extent[1], topo[mfid].extent[2],
+                              topo[mfid].extent[3], topo[mfid].delta[0], topo[mfid].delta[1]]
+            
+            # Whether m-th topo interects with cell
+            mark = self.intersection(cell, topoparam_mfid)
+            if mark[0]:
+                
+                # Whether the cell is completely overlapped by m-th topo
+                if mark[1] == cellarea:
+                    topoint = topoint + self.topointegral(mark[2], topoparam_mfid, topo[mfid].Z, 1, 1)
+                    return topoint
+                else:
+                    
+                    # If the cell is not completely overlapped by m-th topo, use recursion
+                    # to compute the integral on the cell
+                    topoint = self.recurintegral(cell, mtopoorder, mtopofiles, m, topo)
+                    return topoint
+        return topoint
+       
 # Define convenience dictionary of URLs for some online DEMs in netCDF form:
 remote_topo_urls = {}
 

From 07f2c06927516ab4b46e28e47e679d4edea17545 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 12:18:02 -0500
Subject: [PATCH 07/22] Add `patch_value` function

---
 src/python/geoclaw/topotools.py | 30 ++++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 00dbe6a73..8c0cce421 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1935,6 +1935,36 @@ def cellintegral(self, cell, mtopoorder, mtopofiles, topo):
                     return topoint
         return topoint
        
+    def patch_value(self, patch, mtopoorder, mtopofiles, topo):
+        r"""
+        Compute every cell's topo density on the specific path given by "patch".
+        
+        :Input:
+        - *patch* (object) Patch.
+        - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+        - *mtopofiles* (int) The number of the topo objects.
+        - *topo* (list) The list of topo objects.
+
+        :Output:
+        - *cell_value* (ndarray(:, :)) Patch's cells' topo density.
+        """
+    
+        x_num = patch.x.shape[0] - 1
+        y_num = patch.y.shape[0] - 1
+        cell_value = numpy.empty((y_num, x_num))
+        
+        for i in range(y_num):
+            for j in range(x_num):
+                
+                # Set cell's parameter
+                cell = [patch.x[j], patch.x[j+1], patch.y[i], patch.y[i+1]]
+                area = (patch.x[j+1] - patch.x[j]) * (patch.y[i+1] - patch.y[i])
+                
+                # Calculate cell's topo density
+                cell_value[i, j] = self.cellintegral(cell, mtopoorder, mtopofiles, topo) / area
+                
+        return cell_value 
+
 # Define convenience dictionary of URLs for some online DEMs in netCDF form:
 remote_topo_urls = {}
 

From 7037bb451a3109579d585c517cc068fda569fd75 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 12:18:38 -0500
Subject: [PATCH 08/22] Add `patch` class

---
 src/python/geoclaw/topotools.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 8c0cce421..a978c9449 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1964,6 +1964,13 @@ def patch_value(self, patch, mtopoorder, mtopofiles, topo):
                 cell_value[i, j] = self.cellintegral(cell, mtopoorder, mtopofiles, topo) / area
                 
         return cell_value 
+    
+    class patch:
+        def __init__(self, x, y, dx, dy):
+            self.x = x
+            self.y = y
+            self.dx = dx
+            self.dy = dy
 
 # Define convenience dictionary of URLs for some online DEMs in netCDF form:
 remote_topo_urls = {}

From 2509931075a8175d092be0b7d48adc03044b90a4 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 12:21:39 -0500
Subject: [PATCH 09/22] from clawpack.geoclaw.data import data

---
 src/python/geoclaw/topotools.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index a978c9449..b72744e0e 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -41,6 +41,7 @@
 import clawpack.geoclaw.data
 import six
 from six.moves import range
+from clawpack.geoclaw.data import Rearth, DEG2RAD, RAD2DEG
 
 # ==============================================================================
 #  Topography Related Functions

From 833a321102ee57276eb4c7221d5eb4d69c11dbc2 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 12:33:09 -0500
Subject: [PATCH 10/22] Fix `patch_value` docstring format

---
 src/python/geoclaw/topotools.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index b72744e0e..c007c4d79 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1941,13 +1941,13 @@ def patch_value(self, patch, mtopoorder, mtopofiles, topo):
         Compute every cell's topo density on the specific path given by "patch".
         
         :Input:
-        - *patch* (object) Patch.
-        - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
-        - *mtopofiles* (int) The number of the topo objects.
-        - *topo* (list) The list of topo objects.
+         - *patch* (object) Patch.
+         - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+         - *mtopofiles* (int) The number of the topo objects.
+         - *topo* (list) The list of topo objects.
 
         :Output:
-        - *cell_value* (ndarray(:, :)) Patch's cells' topo density.
+         - *cell_value* (ndarray(:, :)) Patch's cells' topo density.
         """
     
         x_num = patch.x.shape[0] - 1

From ab347e92d2dd2326ceadea245c24002ef9e6b458 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 12:35:42 -0500
Subject: [PATCH 11/22] Add `test_integral`

---
 tests/test_topotools.py | 95 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 95 insertions(+)

diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index ed07f583f..35f1f5843 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -427,6 +427,101 @@ def plot_kahului():
     plt.savefig(fname)
     print("Created ",fname)
 
+def test_integral(func, mfile, funcflag, plotflag):
+    r"""
+    test_integral is used to test a set of functions which compute 
+    cell integrals and cell density in topotools.py.
+    
+    :Input:
+     - *func* (function).
+     - *mfile* (int) The number of the topo objects.
+     - *funcflag* (bool) Whether the function, "func" is discontinuous.
+     - *plotflag* (bool) Whether make plots.
+    """
+
+    dx = []
+    data = []
+    topo = []
+    mtopoorder = []
+    
+    # Set boundary for coarsest topo
+    xlow = numpy.random.random() * 10
+    xhi = xlow + 3
+    xarray = numpy.linspace(xlow, xhi, 100)
+    ylow = numpy.random.random() * 10
+    yhi = ylow + 3
+    yarray = numpy.linspace(ylow, yhi, 100)
+    
+    # Set topo data randomly
+    for i in range(mfile):
+        if i == 0:
+            
+            # Set the topo data set which covers the whole patch
+            x = numpy.linspace(xlow, xhi, 100)
+            y = numpy.linspace(ylow, yhi, 100)
+            dx.append(x[1] - x[0])
+        else:
+            
+            # Set x coordinate of the topo data set randomly
+            n1 = int(numpy.random.random() * 70)
+            n2 = n1 + 10 + int(numpy.random.random() * (85 - n1))
+            x1 = xarray[n1]; x2 = xarray[n2]
+            dx.append(dx[i-1] - 0.0012)
+            mx = int((x2 - x1) / dx[i])
+            x2 = x1 + dx[i] * (mx - 1)
+            x = numpy.linspace(x1, x2, mx)
+
+            # Set y coordinate of the topo data set randomly
+            m1 = int(numpy.random.random() * 70)
+            m2 = m1 + 10 + int(numpy.random.random() * (85 - m1))
+            y1 = yarray[m1]; y2 = yarray[m2]
+            my = int((y2 - y1) / dx[i])
+            y2 = y1 + dx[i] * (my - 1)
+            y = numpy.linspace(y1, y2, my)
+        
+        # Set Topography objects parameters
+        topo1 = Topography()
+        
+        # Whether the function is discontinuous function
+        if funcflag == True:
+            z = numpy.empty((len(x), len(y)))
+            for m in range(len(x)):
+                for n in range(len(y)):
+                    z[m][n] = func(x[m], y[n])
+            topo1.Z = z
+        else:
+            topo1.x = numpy.flip(y)
+            topo1.y = x
+            topo1.Z = func(topo1.Y, topo1.X)
+        topo1.x = x; topo1.y = y
+        mtopoorder.append(mfile - 1 - i)
+        topo.append(topo1)
+    
+    # Set patch data
+    patch_x = numpy.linspace(xlow + 1, xhi - 1, 5)
+    patch_y = numpy.linspace(ylow + 1, yhi - 1, 4)
+    patch_dx = patch_x[1] - patch_x[0]
+    patch_dy = patch_y[1] - patch_y[0]
+    patch1 = Topography.patch(patch_x, patch_y, patch_dx, patch_dy)
+    
+    # Accurate cell value
+    real_value = numpy.empty((len(patch1.y) - 1, len(patch1.x) - 1))
+    for i in range(len(patch1.y) - 1):
+        for j in range(len(patch1.x) - 1):
+            area = (patch1.x[j+1] - patch1.x[j]) * (patch1.y[i+1] - patch1.y[i])
+            real_value[i][j] = dblquad(func, patch1.y[i], patch1.y[i+1], 
+                                             patch1.x[j], patch1.x[j+1])[0] / float(area)
+    
+    # Cell value calculated by functions
+    calculated_value = Topography().patch_value(patch1, mtopoorder, mfile, topo)
+        
+    # Whether calculated value achieve the expected resolution
+    npt.assert_almost_equal(real_value, calculated_value, decimal=3)
+    
+    # Whether make plots
+    if plotflag == True:
+        plot(patch_x, patch_y, real_value, calculated_value)
+        
 
 if __name__ == "__main__":
     if len(sys.argv) > 1:

From 4cd81d22c0eff6d0aa6f786f18ce91cdc4511f9d Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 12:44:40 -0500
Subject: [PATCH 12/22] Add `test_integral_plot` for `test_integral`

---
 tests/test_topotools.py | 36 ++++++++++++++++++++++++++++++++++--
 1 file changed, 34 insertions(+), 2 deletions(-)

diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index 35f1f5843..96c2abe0c 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -520,8 +520,40 @@ def test_integral(func, mfile, funcflag, plotflag):
     
     # Whether make plots
     if plotflag == True:
-        plot(patch_x, patch_y, real_value, calculated_value)
-        
+        test_integral_plot(patch_x, patch_y, real_value, calculated_value)
+ 
+def test_integral_plot(x, y, accurate, calculated):
+    r"""
+    Plot function for `test_integral`.
+    
+    :Input:
+     - *x* (ndarray(:)).
+     - *y* (ndarray(:)).
+     - *real* (ndarray(:, :)) Accurate cell value.
+     - *calculated* (ndarray(:, :)) Cell value calculated by functions.
+     
+    :Output:
+     - *fig1* (figure) Figure for accurate cell value.
+     - *fig2* (figure) Figure for calculated cell value.
+     - *fig3* (figure) Figure for error.
+    """
+
+    fig1, plot1 = plt.subplots(figsize=(4,3))
+    im1 = plot1.pcolor(x, y, accurate)
+    plot1.set_title("accurate value")
+    plt.colorbar(im1, ax = plot1)
+    
+    fig2, plot2 = plt.subplots(figsize=(4,3))
+    im2 = plot2.pcolor(x, y, calculated)
+    plot2.set_title("calculated value")
+    plt.colorbar(im2, ax = plot2)
+    
+    fig3, plot3 = plt.subplots(figsize=(4,3))
+    im3 = plot3.pcolor(x, y, numpy.abs(calculated - accurate))
+    plot3.set_title("Error")
+    plt.colorbar(im3, ax = plot3)
+    
+    return fig1, fig2, fig3
 
 if __name__ == "__main__":
     if len(sys.argv) > 1:

From 66211e6f4a1cf47afe121b56cfd124ac82f96782 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 12:52:17 -0500
Subject: [PATCH 13/22] Add import

---
 tests/test_topotools.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index 96c2abe0c..ddbf14adc 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -20,7 +20,11 @@
 
 import clawpack.geoclaw.topotools as topotools
 import clawpack.clawutil.data
+import numpy.testing as npt
+import matplotlib.pyplot as plt
 from six.moves import range
+from clawpack.geoclaw.topotools import Topography
+from scipy.integrate import dblquad
 
 # Set local test directory to get local files
 testdir = os.path.dirname(__file__)

From 2aa438bf775fe6404ab6965a9bf22880b7dd2bf8 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Sun, 15 Dec 2019 13:00:00 -0500
Subject: [PATCH 14/22] Update import

---
 tests/test_topotools.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index ddbf14adc..c6e895ce0 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -21,7 +21,6 @@
 import clawpack.geoclaw.topotools as topotools
 import clawpack.clawutil.data
 import numpy.testing as npt
-import matplotlib.pyplot as plt
 from six.moves import range
 from clawpack.geoclaw.topotools import Topography
 from scipy.integrate import dblquad
@@ -542,6 +541,11 @@ def test_integral_plot(x, y, accurate, calculated):
      - *fig3* (figure) Figure for error.
     """
 
+    try:
+        import matplotlib
+    except ImportError:
+        raise nose.SkipTest("Skipping test since matplotlib not found.")
+    
     fig1, plot1 = plt.subplots(figsize=(4,3))
     im1 = plot1.pcolor(x, y, accurate)
     plot1.set_title("accurate value")

From 5f66683abf242c43896aa59683ef85b59c38f9ab Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Fri, 20 Dec 2019 10:40:31 -0500
Subject: [PATCH 15/22] Update the set of functions for cell integral

Update the docstring. Put the  functions outside the Topography class.
---
 src/python/geoclaw/topotools.py | 695 ++++++++++++++++----------------
 1 file changed, 357 insertions(+), 338 deletions(-)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index c007c4d79..8a36d12f6 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1595,383 +1595,402 @@ def make_shoreline_xy(self, sea_level=0):
                            numpy.array([numpy.nan,numpy.nan]), c.allsegs[0][k]))
         return shoreline_xy
        
-    def bilinearintegral(self, domain, topodomain, corners):
-        r"""
-        Suppose corners are (x1, y1), (x1, y2), (x2, y1), (x1, y2). Get a function,
-        f(xi, eta) = a * (xi - x1) / (x2 - x1) + b * (eta - y1) / (y2 - y1)
-                     + c * (xi - x1) * (eta - y1) / {(x2 - x1) * (y2 - y1)} + d,
-        to represent topo of topomain by bilinear interpolation with four corner
-        points. Then integrate f(xi, eta) on the rectangular domain which is the
-        intersection of domain and topodomain.
-        
-        :Input:
-         - *domain* (ndarray(:))
-         - *topodomain* (ndarray(:))
-         - *corners* (ndarray(:, :)) Four corner points of the topodomain.
+    
+def intersection(rect1, rect2):
+    r"""
+    Whether two domains, rect1 and rect2 intersect. If they intersect,
+    find the intersection's boundary and area.
+
+    :Input:
+     - *rect1* (ndarray(:))
+     - *rect2* (ndarray(:))
+
+    :Output:
+     - *mark* (bool) If rect1 and rect2 intersect, mark = True. Otherwise,
+       mark = False.
+     - *area* (float) The area of the intersection of rect1 and rect2.
+     - *bound* (ndarray(:)) The boundary of the intersection of rect1 and rect2.
+    """
+
+    # Set boundary for rect1
+    x1_low = rect1[0]; x1hi = rect1[1]
+    y1_low = rect1[2]; y1hi = rect1[3]
+
+    # Set boundary for rect2
+    x2low = rect2[0]; x2hi = rect2[1]
+    y2low = rect2[2]; y2hi = rect2[3]
+
+    # Boundary of the intersection part
+    xintlow = max(x1_low, x2low)
+    xinthi = min(x1hi, x2hi)
+    yintlow = max(y1_low, y2low)
+    yinthi = min(y1hi, y2hi)
+
+    # Whether rect1 and rect2 intersect
+    if xinthi > xintlow and yinthi > yintlow:
+        area = (xinthi - xintlow) * (yinthi - yintlow)
+        mark = True
+        bound = [xintlow, xinthi, yintlow, yinthi]
+    else:
+        area = 0.0
+        mark = False
+        bound = [-1, -1, -1, -1]
+
+    return mark, area, bound
 
-        :Output:
-         - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
-           the intersection of domain and topodomain.
-        """
     
-        # Set boundary for integral
-        bound = self.intersection(domain, topodomain)[2]
-        
-        # Find terms for the integration
-        deltax = topodomain[1] - topodomain[0]
-        deltay = topodomain[3] - topodomain[2]
-        area = (bound[3] - bound[2]) * (bound[1] - bound[0])
-        sumxi = (bound[1] + bound[0] - 2.0 * topodomain[0]) / deltax
-        sumeta = (bound[3] + bound[2] - 2.0 * topodomain[2]) / deltay
-        
-        # Find coefficients of the function, f(xi, eta)
-        a = corners[1][0] - corners[0][0]
-        b = corners[0][1] - corners[0][0]
-        c = corners[1][1] - corners[1][0] - corners[0][1] + corners[0][0]
-        d = corners[0][0]
-        
-        bilinearintegral = (0.5 * (a * sumxi + b * sumeta) +
-                            0.25 * c * sumxi * sumeta + d) * area
+def bilinearintegral(domain, topodomain, corners):
+    r"""
+    bilinearintegral integrates a surface over a rectangular region which is the
+    intersection of the surface defined by "domain" and the surface defined by
+    "topodomain". The surface integrated is defined by the function obtained by
+    bilinear interpolation through corners of the Cartesian grid.
+    
+    Suppose corners are (x1, y1), (x1, y2), (x2, y1), (x1, y2). Get a function,
+    f(xi, eta) = a * (xi - x1) / (x2 - x1) + b * (eta - y1) / (y2 - y1)
+                 + c * (xi - x1) * (eta - y1) / {(x2 - x1) * (y2 - y1)} + d,
+    to represent topo of topomain by bilinear interpolation with four corner
+    points. Then integrate f(xi, eta) on the rectangular domain which is the
+    intersection of domain and topodomain.
     
-        return bilinearintegral
+    :Input:
+     - *domain* (ndarray(:))
+     - *topodomain* (ndarray(:))
+     - *corners* (ndarray(:, :)) Four corner points of the topodomain.
 
-    def bilinearintegral_s(self, domain, topodomain, corners):
-        r"""
-        Suppose corners are (x1, y1), (x1, y2), (x2, y1), (x1, y2). Get a function,
-        f(theta, phi) = a * (theta - x1) + b * (phi - y1)
-                        + c * (theta - x1) * (phi - y1) + d,
-        where theta is longitute and phi is latitude, to represent the topo of
-        topomain by bilinear interpolation with four corner points. Then integrate
-        f(theta, phi) on the domain which is the intersection of the domain
-        and the topodomain on the sphere.
-        
-        :Input:
-         - *domain* (ndarray(:))
-         - *topodomain* (ndarray(:))
-         - *corners* (ndarray(:, :)) Four corner points of the topodomain.
+    :Output:
+     - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
+       the intersection of domain and topodomain.
+    """
 
-        :Output:
-         - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
-           the intersection of domain and topodomain.
-        """
+    # Set boundary for integral
+    bound = intersection(domain, topodomain)[2]
     
-        #Set boundary parameter
-        bound = self.intersection(domain, topodomain)[2]
-
-        # Find terms for the integration
-        xdiffhi = bound[1] - topodomain[0]
-        xdifflow = bound[0] - topodomain[0]
-        ydiffhi = bound[3] - topodomain[2]
-        ydifflow = bound[2] - topodomain[2]
-        xdiff2 = 0.5 * (xdiffhi**2 - xdifflow**2)
-        intdx = bound[1] - bound[0]
-        deltax = topodomain[1] - topodomain[0]
-        deltay = topodomain[3] - topodomain[2]
-        
-        d2r = DEG2RAD; r2d = RAD2DEG
+    # Find terms for the integration
+    deltax = topodomain[1] - topodomain[0]
+    deltay = topodomain[3] - topodomain[2]
+    area = (bound[3] - bound[2]) * (bound[1] - bound[0])
+    sumxi = (bound[1] + bound[0] - 2.0 * topodomain[0]) / deltax
+    sumeta = (bound[3] + bound[2] - 2.0 * topodomain[2]) / deltay
+    
+    # Find coefficients of the function, f(xi, eta)
+    a = corners[1][0] - corners[0][0]
+    b = corners[0][1] - corners[0][0]
+    c = corners[1][1] - corners[1][0] - corners[0][1] + corners[0][0]
+    d = corners[0][0]
+    
+    bilinearintegral = (0.5 * (a * sumxi + b * sumeta) +
+                        0.25 * c * sumxi * sumeta + d) * area
+
+    return bilinearintegral
+    
+    
+def bilinearintegral_s(domain, topodomain, corners):
+    r"""
+    bilinearintegral_s integrates a surface over a sphere which is the intersection
+    of the surface defined by "domain" and the surface defined by "topodomain".
+    The surface integrated is defined by the function obtained by bilinear
+    interpolation through corners of the polar grid.
+    
+    Suppose corners are (x1, y1), (x1, y2), (x2, y1), (x1, y2). Get a function,
+    f(theta, phi) = a * (theta - x1) + b * (phi - y1)
+                    + c * (theta - x1) * (phi - y1) + d,
+    where theta is longitute and phi is latitude, to represent the topo of
+    topomain by bilinear interpolation with four corner points. Then integrate
+    f(theta, phi) on the domain which is the intersection of the domain
+    and the topodomain on the sphere.
+    
+    :Input:
+     - *domain* (ndarray(:))
+     - *topodomain* (ndarray(:))
+     - *corners* (ndarray(:, :)) Four corner points of the topodomain.
+
+    :Output:
+     - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
+       the intersection of domain and topodomain.
+    """
 
-        cbsinint = (r2d * numpy.cos(d2r * bound[3]) + ydiffhi * numpy.sin(d2r * bound[3]))
-        - (r2d * numpy.cos(d2r * bound[2]) + ydifflow * numpy.sin(d2r * bound[2]))
+    #Set boundary parameter
+    bound = intersection(domain, topodomain)[2]
+
+    # Find terms for the integration
+    xdiffhi = bound[1] - topodomain[0]
+    xdifflow = bound[0] - topodomain[0]
+    ydiffhi = bound[3] - topodomain[2]
+    ydifflow = bound[2] - topodomain[2]
+    xdiff2 = 0.5 * (xdiffhi**2 - xdifflow**2)
+    intdx = bound[1] - bound[0]
+    deltax = topodomain[1] - topodomain[0]
+    deltay = topodomain[3] - topodomain[2]
+    
+    d2r = DEG2RAD; r2d = RAD2DEG
 
-        adsinint = r2d * (numpy.sin(d2r * bound[3]) - numpy.sin(d2r * bound[2]))
+    cbsinint = (r2d * numpy.cos(d2r * bound[3]) + ydiffhi * numpy.sin(d2r * bound[3]))
+    - (r2d * numpy.cos(d2r * bound[2]) + ydifflow * numpy.sin(d2r * bound[2]))
 
-        # Find coefficients of the function, f(theta, phi)
-        a = (corners[1][0] - corners[0][0]) / deltax
-        b = (corners[0][1] - corners[0][0]) / deltay
-        c = (corners[1][1] - corners[1][0] - corners[0][1] + corners[0][0]) / (deltax * deltay)
-        d = corners[0][0]
+    adsinint = r2d * (numpy.sin(d2r * bound[3]) - numpy.sin(d2r * bound[2]))
 
-        bilinearintegral_s = ((a * xdiff2 + d * intdx) * adsinint +
-                              r2d * (c * xdiff2 + b * intdx) * cbsinint) * (Rearth * d2r)**2
+    # Find coefficients of the function, f(theta, phi)
+    a = (corners[1][0] - corners[0][0]) / deltax
+    b = (corners[0][1] - corners[0][0]) / deltay
+    c = (corners[1][1] - corners[1][0] - corners[0][1] + corners[0][0]) / (deltax * deltay)
+    d = corners[0][0]
 
-        return bilinearintegral_s
+    bilinearintegral_s = ((a * xdiff2 + d * intdx) * adsinint +
+                          r2d * (c * xdiff2 + b * intdx) * cbsinint) * (Rearth * d2r)**2
 
-    def topointegral(self, domain, topoparam, z, intmethod, coordinate_system):
-        r"""
-        topointegral integrates a surface over a rectangular region which is
-        the intersection of the domain and topo which is defined by "topoparam".
-        The surface integrated is defined by a function gotten by bilinear
-        interpolation through corners of the Cartesian grid.
+    return bilinearintegral_s
         
-        :Input:
-         - *domain* (ndarray(:))
-         - *topoparam* (ndarray(:)) topoparam includes topo's x and y coordinates'
-           minimum and maximum, and topo's grid size, dx and dy.
-         - *z* (ndarray(:, :)) Topo data of domain represented by topoparam.
-         - *intmethod* (int) The method of integral.
-         - *coordinate_system* (int) The type of coordinate system.
 
-        :Output:
-         - *theintegral* (float) The integral of a surface over a rectangular
-           region which is the intersection of the domain and "topoparam".
-        """
+def topointegral(domain, topoparam, z, intmethod, coordinate_system):
+    r"""
+    topointegral integrates a surface over a rectangular region which is
+    the intersection of the surface defined by "domain" and the surface
+    defined by "topoparam". The surface integrated is defined by the function
+    obtained by bilinear interpolation through nodes of the Cartesian grid.
+    
+    :Input:
+     - *domain* (ndarray(:))
+     - *topoparam* (ndarray(:)) topoparam includes topo's x and y coordinates'
+       minimum and maximum, and topo's grid size, dx and dy.
+     - *z* (ndarray(:, :)) Topo data of domain represented by topoparam.
+     - *intmethod* (int) The method of integral.
+     - *coordinate_system* (int) The type of coordinate system.
+
+    :Output:
+     - *theintegral* (float) The integral of a surface over a rectangular
+       region which is the intersection of the domain and "topoparam".
+    """
+    
+    theintegral = 0.0
+    
+    # Set topo's parameter
+    x1 = topoparam[0]; x2 = topoparam[1]
+    y1 = topoparam[2]; y2 = topoparam[3]
+    topodx = topoparam[4]; topody = topoparam[5]
+    topomx = numpy.array(z).shape[0]
+    topomy = numpy.array(z).shape[1]
+    
+    # Find the intersection of the domain and topo
+    bound = intersection(domain, topoparam)[2]
+    xlow = bound[0]; xhi = bound[1]
+    ylow = bound[2]; yhi = bound[3]
+    
+    if intmethod == 1:
         
-        theintegral = 0.0
+        # Find topo's start and end points for integral
+        # The x coodinate of the topo's start point
+        istart = -1
+        for i in range(topomx):
+            if (x1 + i * topodx) > xlow:
+                istart = i
+                break
         
-        # Set topo's parameter
-        x1 = topoparam[0]; x2 = topoparam[1]
-        y1 = topoparam[2]; y2 = topoparam[3]
-        topodx = topoparam[4]; topody = topoparam[5]
-        topomx = numpy.array(z).shape[0]
-        topomy = numpy.array(z).shape[1]
+        # The y coodinate of the topo's start point
+        jstart = -1
+        for j in range(topomy):
+            if (y1 + j * topody)  > ylow:
+                jstart = j
+                break
         
-        # Find the intersection of the domain and topo
-        bound = self.intersection(domain, topoparam)[2]
-        xlow = bound[0]; xhi = bound[1]
-        ylow = bound[2]; yhi = bound[3]
+        # The x coodinate of the topo's end point
+        iend = topomx
+        for m in range(topomx):
+            if (x1 + m * topodx) >= xhi:
+                iend = m
+                break
         
-        if intmethod == 1:
-            
-            # Find topo's start and end points for integral
-            # The x coodinate of the topo's start point
-            istart = -1
-            for i in range(topomx):
-                if (x1 + i * topodx) > xlow:
-                    istart = i
-                    break
-            
-            # The y coodinate of the topo's start point
-            jstart = -1
-            for j in range(topomy):
-                if (y1 + j * topody)  > ylow:
-                    jstart = j
-                    break
-            
-            # The x coodinate of the topo's end point
-            iend = topomx
-            for m in range(topomx):
-                if (x1 + m * topodx) >= xhi:
-                    iend = m
-                    break
-            
-            # The y coodinate of the topo's end point
-            jend = topomy
-            for n in range(topomy):
-                if (y1 + n * topody) >= yhi:
-                    jend = n
-                    break
-                            
-            # Prevent overflow
-            jstart = max(jstart, 1)
-            istart = max(istart, 1)
-            jend = min(jend, topomy - 1)
-            iend = min(iend, topomx - 1)
-            
-            # Topo integral
-            for jj in range(jstart, jend + 1):
-                yint1 = y1 + (jj - 1.0) * topody
-                yint2 = y1 + jj * topody
-
-                for ii in range(istart, iend + 1):
-                    xint1 = x1 + (ii - 1.0) * topodx
-                    xint2 = x1 + ii * topodx
-
-                    # Four corners of topo's grid
-                    corners = numpy.empty((2, 2))
-                    corners[0][0] = z[ii-1][topomy - jj]
-                    corners[0][1] = z[ii-1][topomy - 1 - jj]
-                    corners[1][0] = z[ii][topomy - jj]
-                    corners[1][1] = z[ii][topomy - 1 - jj]
-
-                    # Parameters of topo's grid integral
-                    topointparam = [xint1, xint2, yint1, yint2]
-
-                    if coordinate_system == 1:
-                        theintegral += self.bilinearintegral(domain, topointparam, corners)
-                    elif coordinate_system == 2:
-                        theintegral += self.bilinearintegral_s(domain, topointparam, corners)
-                    else:
-                         print("TOPOINTEGRAL: coordinate_system error")
-
-        else:
-            print("TOPOINTEGRAL: only intmethod = 1,2 is supported")
+        # The y coodinate of the topo's end point
+        jend = topomy
+        for n in range(topomy):
+            if (y1 + n * topody) >= yhi:
+                jend = n
+                break
+                        
+        # Prevent overflow
+        jstart = max(jstart, 1)
+        istart = max(istart, 1)
+        jend = min(jend, topomy - 1)
+        iend = min(iend, topomx - 1)
+        
+        # Topo integral
+        for jj in range(jstart, jend + 1):
+            yint1 = y1 + (jj - 1.0) * topody
+            yint2 = y1 + jj * topody
+
+            for ii in range(istart, iend + 1):
+                xint1 = x1 + (ii - 1.0) * topodx
+                xint2 = x1 + ii * topodx
+
+                # Four corners of topo's grid
+                corners = numpy.empty((2, 2))
+                corners[0][0] = z[ii-1][topomy - jj]
+                corners[0][1] = z[ii-1][topomy - 1 - jj]
+                corners[1][0] = z[ii][topomy - jj]
+                corners[1][1] = z[ii][topomy - 1 - jj]
+
+                # Parameters of topo's grid integral
+                topointparam = [xint1, xint2, yint1, yint2]
+
+                if coordinate_system == 1:
+                    theintegral += bilinearintegral(domain, topointparam, corners)
+                elif coordinate_system == 2:
+                    theintegral += bilinearintegral_s(domain, topointparam, corners)
+                else:
+                     print("TOPOINTEGRAL: coordinate_system error")
 
-        return theintegral
+    else:
+        print("TOPOINTEGRAL: only intmethod = 1,2 is supported")
+
+    return theintegral
+
+
+def recurintegral(domain, mtopoorder, mtopofiles, m, topo):
+    r"""
+    Compute the integral of topo over the rectangle domain defined by "domain".
+    Find the index of the topo whose resolution is m by "mtopoorder". Use this
+    index to find corresponding topo object in the list "topo".
+
+    The main call to this recursive function has corners of a grid cell for the
+    rectangle and m = 1 in order to compute the integral over the cell using all
+    topo objects.
+
+    The recursive strategy is to first compute the integral using only the topo
+    object of m+1 resolution. Then apply corrections due to adding m resolution
+    topo object.
+     
+    Corrections are needed if the new topo object intersects the grid cell.
+    Let the intersection be represented by mark2[2]. Two corrections are needed.
+    First, subtract out the integral over the rectangle "mark2[2]" computed with
+    topo objects mtopoorder(mtopofiles) to mtopoorder(m+1), and then adding in
+    the integral over this same region using the topo object mtopoorder(m).
+    
+    :Input:
+     - *domain* (ndarray(:)) The specific surface where the integral is computed.
+     - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+     - *mtopofiles* (int) The number of the topo objects.
+     - *m* (int) Resolution of the topo objects.
+     - *topo* (list) The list of topo objects.
 
-    def intersection(self, rect1, rect2):
-        r"""
-        Whether two domains, rect1 and rect2 intersect. If they intersect,
-        find the intersection's boundary and area.
-        
-        :Input:
-         - *rect1* (ndarray(:))
-         - *rect2* (ndarray(:))
+    :Output:
+     - *integral* (float) Integral.
+    """
+    
+    # The index of the topo object whose resolution is m
+    mfid = mtopoorder[m]
+    
+    # The parameters of the m resolution topo object
+    topoparam_mfid = [topo[mfid].extent[0], topo[mfid].extent[1], topo[mfid].extent[2],
+                      topo[mfid].extent[3], topo[mfid].delta[0], topo[mfid].delta[1]]
+    
+    # Innermost step of recursion reaches this point, only using coarsest topo grid
+    # compute directly
+    if m == mtopofiles - 1:
+        mark1 = intersection(domain, topoparam_mfid)
+        if mark1[0]:
+            integral = topointegral(mark1[2], topoparam_mfid, topo[mfid].Z, 1, 1)
+        else:
+            integral = 0.0
+    else:
+        int1 = recurintegral(domain, mtopoorder, mtopofiles, m+1, topo)
 
-        :Output:
-         - *mark* (bool) If rect1 and rect2 intersect, mark = True. Otherwise,
-           mark = False.
-         - *area* (float) The area of the intersection of rect1 and rect2.
-         - *bound* (ndarray(:)) The boundary of the intersection of rect1 and rect2.
-        """
-        
-        # Set boundary for rect1
-        x1_low = rect1[0]; x1hi = rect1[1]
-        y1_low = rect1[2]; y1hi = rect1[3]
-        
-        # Set boundary for rect2
-        x2low = rect2[0]; x2hi = rect2[1]
-        y2low = rect2[2]; y2hi = rect2[3]
+        # Whether new topo object intersects the grid cell
+        mark2 = intersection(domain, topoparam_mfid)
         
-        # Boundary of the intersection part
-        xintlow = max(x1_low, x2low)
-        xinthi = min(x1hi, x2hi)
-        yintlow = max(y1_low, y2low)
-        yinthi = min(y1hi, y2hi)
-
-        # Whether rect1 and rect2 intersect
-        if xinthi > xintlow and yinthi > yintlow:
-            area = (xinthi - xintlow) * (yinthi - yintlow)
-            mark = True
-            bound = [xintlow, xinthi, yintlow, yinthi]
+        # The new topo object intersects the grid cell. Corrections are needed
+        if mark2[1] > 0:
+            int2 = recurintegral(mark2[2], mtopoorder, mtopofiles, m+1, topo)
+            int3 = topointegral(mark2[2], topoparam_mfid, topo[mfid].Z, 1, 1)
+            integral = int1 - int2 + int3
         else:
-            area = 0.0
-            mark = False
-            bound = [-1, -1, -1, -1]
+            integral = int1
+    return integral
+        
 
-        return mark, area, bound
+def cellintegral(cell, mtopoorder, mtopofiles, topo):
+    r"""
+    Compute the integral on the surface defined by "cell" with the topo objects
+    list "topo".
+    
+    :Input:
+    - *cell* (ndarray(:)) The boundary of specific surface where the integral
+      is computed.
+    - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+    - *mtopofiles* (int) The number of the topo objects.
+    - *topo* (list) The list of topo objects.
 
-    def recurintegral(self, domain, mtopoorder, mtopofiles, m, topo):
-        r"""
-        Compute the integral of topo over the rectangle domain define by "domain".
-        Find the index of the topo whose resolution is m by mtopoorder. Use this
-        index to find corrosponding topo object in the list "topo".
-
-        The main call to this recursive function has corners of a grid cell for the
-        rectangle and m = 1 in order to compute the integral over the cell
-        using all topo objects.
-
-        The recursive strategy is to first compute the integral using only the topo
-        object of m+1 resolution. Then apply corrections due to adding m resolution
-        topo object.
-         
-        Corrections are needed if the new topo object intersects the grid cell.
-        Let the intersection be represented by mark2[2]. Two corrections are needed.
-        First, subtract out the integral over the rectangle "mark2[2]" computed using
-        topo objects mtopoorder(mtopofiles) to mtopoorder(m+1), and then adding in
-        the integral over this same region using the topo object mtopoorder(m).
-        
-        :Input:
-         - *domain* (ndarray(:)) The specific surface where the integral is computed.
-         - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
-         - *mtopofiles* (int) The number of the topo objects.
-         - *m* (int) Resolution of the topo objects.
-         - *topo* (list) The list of topo objects.
+    :Output:
+    - *topoint* (float) The integral on the cell, "cell".
+    """
+    
+    # Initialization
+    topoint = 0.0
 
-        :Output:
-         - *integral* (float) Integral.
-        """
-        
-        # The index of the topo object whose resolution is m
+    for m in range(mtopofiles):
+    
+        # The index of the topo file object resolution is m
         mfid = mtopoorder[m]
+        cellarea = (cell[1] - cell[0]) * (cell[3] - cell[2])
         
-        # The parameters of the m resolution topo object
+        # The parameters of this topo object
         topoparam_mfid = [topo[mfid].extent[0], topo[mfid].extent[1], topo[mfid].extent[2],
                           topo[mfid].extent[3], topo[mfid].delta[0], topo[mfid].delta[1]]
         
-        # Innermost step of recursion reaches this point, only using coarsest topo grid
-        # compute directly
-        if m == mtopofiles - 1:
-            mark1 = self.intersection(domain, topoparam_mfid)
-            if mark1[0]:
-                integral = self.topointegral(mark1[2], topoparam_mfid, topo[mfid].Z, 1, 1)
-            else:
-                integral = 0.0
-        else:
-            int1 = self.recurintegral(domain, mtopoorder, mtopofiles, m+1, topo)
-
-            # Whether new topo object intersects the grid cell
-            mark2 = self.intersection(domain, topoparam_mfid)
+        # Whether m-th topo interects with cell
+        mark = intersection(cell, topoparam_mfid)
+        if mark[0]:
             
-            # The new topo object intersects the grid cell. Corrections are needed
-            if mark2[1] > 0:
-                int2 = self.recurintegral(mark2[2], mtopoorder, mtopofiles, m+1, topo)
-                int3 = self.topointegral(mark2[2], topoparam_mfid, topo[mfid].Z, 1, 1)
-                integral = int1 - int2 + int3
+            # Whether the cell is completely overlapped by m-th topo
+            if mark[1] == cellarea:
+                topoint = topoint + topointegral(mark[2], topoparam_mfid, topo[mfid].Z, 1, 1)
+                return topoint
             else:
-                integral = int1
-        return integral
+                
+                # If the cell is not completely overlapped by m-th topo, use recursion
+                # to compute the integral on the cell
+                topoint = recurintegral(cell, mtopoorder, mtopofiles, m, topo)
+                return topoint
+    return topoint
 
-    def cellintegral(self, cell, mtopoorder, mtopofiles, topo):
-        r"""
-        Compute the integral on the surface defined by "cell" with the topo objects list
-        "topo".
-        
-        :Input:
-        - *cell* (ndarray(:)) The boundary of specific surface where the integral
-          is computed.
-        - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
-        - *mtopofiles* (int) The number of the topo objects.
-        - *topo* (list) The list of topo objects.
 
-        :Output:
-        - *topoint* (float) The integral on the cell, "cell".
-        """
-        
-        # Initialization
-        topoint = 0.0
+class patch:
+    def __init__(self, x, y, dx, dy):
+        self.x = x
+        self.y = y
+        self.dx = dx
+        self.dy = dy
 
-        for m in range(mtopofiles):
-        
-            # The index of the topo file object resolution is m
-            mfid = mtopoorder[m]
-            cellarea = (cell[1] - cell[0]) * (cell[3] - cell[2])
+
+def cell_average_patch(patch, mtopoorder, mtopofiles, topo):
+    r"""
+    Compute every cell's average of the specific path given by "patch".
+
+    :Input:
+    - *patch* (object) Patch.
+    - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+    - *mtopofiles* (int) The number of the topo objects.
+    - *topo* (list) The list of topo objects.
+
+    :Output:
+    - *cell_average* (ndarray(:, :)) Patch's every cell average.
+    """
+
+    x_num = patch.x.shape[0] - 1
+    y_num = patch.y.shape[0] - 1
+    cell_average = numpy.empty((y_num, x_num))
+
+    for i in range(y_num):
+        for j in range(x_num):
             
-            # The parameters of this topo object
-            topoparam_mfid = [topo[mfid].extent[0], topo[mfid].extent[1], topo[mfid].extent[2],
-                              topo[mfid].extent[3], topo[mfid].delta[0], topo[mfid].delta[1]]
+            # Set cell's parameter
+            cell = [patch.x[j], patch.x[j+1], patch.y[i], patch.y[i+1]]
+            area = (patch.x[j+1] - patch.x[j]) * (patch.y[i+1] - patch.y[i])
             
-            # Whether m-th topo interects with cell
-            mark = self.intersection(cell, topoparam_mfid)
-            if mark[0]:
-                
-                # Whether the cell is completely overlapped by m-th topo
-                if mark[1] == cellarea:
-                    topoint = topoint + self.topointegral(mark[2], topoparam_mfid, topo[mfid].Z, 1, 1)
-                    return topoint
-                else:
-                    
-                    # If the cell is not completely overlapped by m-th topo, use recursion
-                    # to compute the integral on the cell
-                    topoint = self.recurintegral(cell, mtopoorder, mtopofiles, m, topo)
-                    return topoint
-        return topoint
-       
-    def patch_value(self, patch, mtopoorder, mtopofiles, topo):
-        r"""
-        Compute every cell's topo density on the specific path given by "patch".
-        
-        :Input:
-         - *patch* (object) Patch.
-         - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
-         - *mtopofiles* (int) The number of the topo objects.
-         - *topo* (list) The list of topo objects.
+            # Calculate every cell average
+            cell_average[i, j] = cellintegral(cell, mtopoorder, mtopofiles, topo) / area
+            
+    return cell_average
 
-        :Output:
-         - *cell_value* (ndarray(:, :)) Patch's cells' topo density.
-        """
-    
-        x_num = patch.x.shape[0] - 1
-        y_num = patch.y.shape[0] - 1
-        cell_value = numpy.empty((y_num, x_num))
-        
-        for i in range(y_num):
-            for j in range(x_num):
-                
-                # Set cell's parameter
-                cell = [patch.x[j], patch.x[j+1], patch.y[i], patch.y[i+1]]
-                area = (patch.x[j+1] - patch.x[j]) * (patch.y[i+1] - patch.y[i])
-                
-                # Calculate cell's topo density
-                cell_value[i, j] = self.cellintegral(cell, mtopoorder, mtopofiles, topo) / area
-                
-        return cell_value 
-    
-    class patch:
-        def __init__(self, x, y, dx, dy):
-            self.x = x
-            self.y = y
-            self.dx = dx
-            self.dy = dy
 
 # Define convenience dictionary of URLs for some online DEMs in netCDF form:
 remote_topo_urls = {}

From 90810d7bf2f391748d0e297f4384bde12de31c1a Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Fri, 20 Dec 2019 10:46:09 -0500
Subject: [PATCH 16/22] Update the test function for cell integral

Update the test function with deterministic data. Test the cell integral and cell averages with continuous and discontinuous functions.
---
 tests/test_topotools.py | 102 ++++++++++++++++------------------------
 1 file changed, 41 insertions(+), 61 deletions(-)

diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index c6e895ce0..2b3074b5d 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -430,16 +430,15 @@ def plot_kahului():
     plt.savefig(fname)
     print("Created ",fname)
 
-def test_integral(func, mfile, funcflag, plotflag):
+def integral_topotool(func, mfile, funcflag):
     r"""
-    test_integral is used to test a set of functions which compute 
-    cell integrals and cell density in topotools.py.
+    integral_topotool is used for testing a set of functions which compute
+    cell integrals and cell averages in topotools.py.
     
     :Input:
      - *func* (function).
      - *mfile* (int) The number of the topo objects.
      - *funcflag* (bool) Whether the function, "func" is discontinuous.
-     - *plotflag* (bool) Whether make plots.
     """
 
     dx = []
@@ -448,35 +447,32 @@ def test_integral(func, mfile, funcflag, plotflag):
     mtopoorder = []
     
     # Set boundary for coarsest topo
-    xlow = numpy.random.random() * 10
-    xhi = xlow + 3
-    xarray = numpy.linspace(xlow, xhi, 100)
-    ylow = numpy.random.random() * 10
-    yhi = ylow + 3
-    yarray = numpy.linspace(ylow, yhi, 100)
+    xlow = 5; xhi = 8
+    xarray = numpy.linspace(xlow, xhi, mfile * 5)
+    ylow = 2; yhi = 5
+    yarray = numpy.linspace(ylow, yhi, mfile * 5)
     
-    # Set topo data randomly
+    # Set a set of topo data
     for i in range(mfile):
         if i == 0:
             
             # Set the topo data set which covers the whole patch
-            x = numpy.linspace(xlow, xhi, 100)
-            y = numpy.linspace(ylow, yhi, 100)
+            x = xarray; y = yarray
             dx.append(x[1] - x[0])
         else:
             
-            # Set x coordinate of the topo data set randomly
-            n1 = int(numpy.random.random() * 70)
-            n2 = n1 + 10 + int(numpy.random.random() * (85 - n1))
+            # Set x coordinate of the topo data set
+            n1 = int(3 * i)
+            n2 = int(3 * i + 1.2 * mfile)
             x1 = xarray[n1]; x2 = xarray[n2]
             dx.append(dx[i-1] - 0.0012)
             mx = int((x2 - x1) / dx[i])
             x2 = x1 + dx[i] * (mx - 1)
             x = numpy.linspace(x1, x2, mx)
 
-            # Set y coordinate of the topo data set randomly
-            m1 = int(numpy.random.random() * 70)
-            m2 = m1 + 10 + int(numpy.random.random() * (85 - m1))
+            # Set y coordinate of the topo data set
+            m1 = int(3 * i)
+            m2 = int(3 * i + 1.2 * mfile)
             y1 = yarray[m1]; y2 = yarray[m2]
             my = int((y2 - y1) / dx[i])
             y2 = y1 + dx[i] * (my - 1)
@@ -501,67 +497,51 @@ def test_integral(func, mfile, funcflag, plotflag):
         topo.append(topo1)
     
     # Set patch data
-    patch_x = numpy.linspace(xlow + 1, xhi - 1, 5)
-    patch_y = numpy.linspace(ylow + 1, yhi - 1, 4)
+    patch_x = numpy.linspace(xlow + 0.1, xhi - 0.1, 5)
+    patch_y = numpy.linspace(ylow + 0.1, yhi - 0.1, 4)
     patch_dx = patch_x[1] - patch_x[0]
     patch_dy = patch_y[1] - patch_y[0]
-    patch1 = Topography.patch(patch_x, patch_y, patch_dx, patch_dy)
+    patch1 = topotools.patch(patch_x, patch_y, patch_dx, patch_dy)
     
     # Accurate cell value
     real_value = numpy.empty((len(patch1.y) - 1, len(patch1.x) - 1))
     for i in range(len(patch1.y) - 1):
         for j in range(len(patch1.x) - 1):
             area = (patch1.x[j+1] - patch1.x[j]) * (patch1.y[i+1] - patch1.y[i])
-            real_value[i][j] = dblquad(func, patch1.y[i], patch1.y[i+1], 
-                                             patch1.x[j], patch1.x[j+1])[0] / float(area)
+            real_value[i][j] = dblquad(func=func, a=patch1.y[i], b=patch1.y[i+1],
+                                       gfun=lambda x: patch1.x[j],
+                                       hfun=lambda x: patch1.x[j+1])[0] / float(area)
     
     # Cell value calculated by functions
-    calculated_value = Topography().patch_value(patch1, mtopoorder, mfile, topo)
+    calculated_value = topotools.cell_average_patch(patch1, mtopoorder, mfile, topo)
         
     # Whether calculated value achieve the expected resolution
     npt.assert_almost_equal(real_value, calculated_value, decimal=3)
-    
-    # Whether make plots
-    if plotflag == True:
-        test_integral_plot(patch_x, patch_y, real_value, calculated_value)
- 
-def test_integral_plot(x, y, accurate, calculated):
+    return None
+
+
+def testcontinuous():
     r"""
-    Plot function for `test_integral`.
-    
-    :Input:
-     - *x* (ndarray(:)).
-     - *y* (ndarray(:)).
-     - *real* (ndarray(:, :)) Accurate cell value.
-     - *calculated* (ndarray(:, :)) Cell value calculated by functions.
-     
-    :Output:
-     - *fig1* (figure) Figure for accurate cell value.
-     - *fig2* (figure) Figure for calculated cell value.
-     - *fig3* (figure) Figure for error.
+    testcontinuous is used to test continuous functions with `integral_topotool`.
     """
-
-    try:
-        import matplotlib
-    except ImportError:
-        raise nose.SkipTest("Skipping test since matplotlib not found.")
     
-    fig1, plot1 = plt.subplots(figsize=(4,3))
-    im1 = plot1.pcolor(x, y, accurate)
-    plot1.set_title("accurate value")
-    plt.colorbar(im1, ax = plot1)
+    integral_topotool(lambda x, y: numpy.sin(x**2 + x * y) + x + y, 20, False)
+    integral_topotool(lambda x, y: numpy.exp(0.5 * x + 0.1 * y), 20, False)
+    integral_topotool(lambda x, y: x**2 + x * y + y**2 + 1, 20, False)
+    return None
     
-    fig2, plot2 = plt.subplots(figsize=(4,3))
-    im2 = plot2.pcolor(x, y, calculated)
-    plot2.set_title("calculated value")
-    plt.colorbar(im2, ax = plot2)
     
-    fig3, plot3 = plt.subplots(figsize=(4,3))
-    im3 = plot3.pcolor(x, y, numpy.abs(calculated - accurate))
-    plot3.set_title("Error")
-    plt.colorbar(im3, ax = plot3)
+def testdiscontinuous():
+    r"""
+    testdiscontinuous is used to test discontinuous function with `integral_topotool`.
+    """
     
-    return fig1, fig2, fig3
+    def fun(x, y):
+        for i in range(20):
+            if int(x * 100) % 20 == i:
+                return 9.99 + 0.001 * i
+    integral_topotool(fun, 20, True)
+    return None
 
 if __name__ == "__main__":
     if len(sys.argv) > 1:

From 74681c3029763c566111cda5bb83dc379e1d4ad2 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Thu, 26 Dec 2019 17:10:29 -0500
Subject: [PATCH 17/22] Fix warning error in the test

---
 tests/test_topotools.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index 2b3074b5d..1271e73c9 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -508,9 +508,15 @@ def integral_topotool(func, mfile, funcflag):
     for i in range(len(patch1.y) - 1):
         for j in range(len(patch1.x) - 1):
             area = (patch1.x[j+1] - patch1.x[j]) * (patch1.y[i+1] - patch1.y[i])
-            real_value[i][j] = dblquad(func=func, a=patch1.y[i], b=patch1.y[i+1],
-                                       gfun=lambda x: patch1.x[j],
-                                       hfun=lambda x: patch1.x[j+1])[0] / float(area)
+            if funcflag:
+                options = {'limit':1000}
+                real_value[i][j] = nquad(func=func, ranges=[[patch1.x[j], patch1.x[j+1]],
+                                         [patch1.y[i], patch1.y[i+1]]], args=None,
+                                         opts=[options,options])[0] / float(area)
+            else:
+                real_value[i][j] = dblquad(func=func, a=patch1.y[i], b=patch1.y[i+1],
+                                           gfun=patch1.x[j],
+                                           hfun=patch1.x[j+1])[0] / float(area)
     
     # Cell value calculated by functions
     calculated_value = topotools.cell_average_patch(patch1, mtopoorder, mfile, topo)

From 6b23e33d9abd9843830d2cb743eea04c00410804 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Thu, 26 Dec 2019 17:14:57 -0500
Subject: [PATCH 18/22] Import nquad

---
 tests/test_topotools.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index 1271e73c9..06df834c7 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -23,7 +23,7 @@
 import numpy.testing as npt
 from six.moves import range
 from clawpack.geoclaw.topotools import Topography
-from scipy.integrate import dblquad
+from scipy.integrate import dblquad, nquad
 
 # Set local test directory to get local files
 testdir = os.path.dirname(__file__)
@@ -482,7 +482,7 @@ def integral_topotool(func, mfile, funcflag):
         topo1 = Topography()
         
         # Whether the function is discontinuous function
-        if funcflag == True:
+        if funcflag:
             z = numpy.empty((len(x), len(y)))
             for m in range(len(x)):
                 for n in range(len(y)):

From e5a338967af4e94a08f38c17f81b7f484289c742 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Thu, 26 Dec 2019 17:45:29 -0500
Subject: [PATCH 19/22] Change patch to a dictionary

---
 src/python/geoclaw/topotools.py | 20 ++++++--------------
 1 file changed, 6 insertions(+), 14 deletions(-)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 8a36d12f6..5251ce5ab 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1953,20 +1953,12 @@ def cellintegral(cell, mtopoorder, mtopofiles, topo):
     return topoint
 
 
-class patch:
-    def __init__(self, x, y, dx, dy):
-        self.x = x
-        self.y = y
-        self.dx = dx
-        self.dy = dy
-
-
 def cell_average_patch(patch, mtopoorder, mtopofiles, topo):
     r"""
-    Compute every cell's average of the specific path given by "patch".
+    Compute every cell's average of the specific patch given by "patch".
 
     :Input:
-    - *patch* (object) Patch.
+    - *patch* (Dictionary) Patch.
     - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
     - *mtopofiles* (int) The number of the topo objects.
     - *topo* (list) The list of topo objects.
@@ -1975,16 +1967,16 @@ def cell_average_patch(patch, mtopoorder, mtopofiles, topo):
     - *cell_average* (ndarray(:, :)) Patch's every cell average.
     """
 
-    x_num = patch.x.shape[0] - 1
-    y_num = patch.y.shape[0] - 1
+    x_num = len(patch['x']) - 1
+    y_num = len(patch['y']) - 1
     cell_average = numpy.empty((y_num, x_num))
 
     for i in range(y_num):
         for j in range(x_num):
             
             # Set cell's parameter
-            cell = [patch.x[j], patch.x[j+1], patch.y[i], patch.y[i+1]]
-            area = (patch.x[j+1] - patch.x[j]) * (patch.y[i+1] - patch.y[i])
+            cell = [patch['x'][j], patch['x'][j+1], patch['y'][i], patch['y'][i+1]]
+            area = (patch['x'][j+1] - patch['x'][j]) * (patch['y'][i+1] - patch['y'][i])
             
             # Calculate every cell average
             cell_average[i, j] = cellintegral(cell, mtopoorder, mtopofiles, topo) / area

From 56cd613eb8025520a9cfc692135b0fbc601865eb Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Thu, 26 Dec 2019 17:50:03 -0500
Subject: [PATCH 20/22] Change patch to a dictionary for test function

---
 tests/test_topotools.py | 27 +++++++++++++--------------
 1 file changed, 13 insertions(+), 14 deletions(-)

diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index 06df834c7..036bc6029 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -446,7 +446,7 @@ def integral_topotool(func, mfile, funcflag):
     topo = []
     mtopoorder = []
     
-    # Set boundary for coarsest topo
+    # Set boundary for the coarsest topo
     xlow = 5; xhi = 8
     xarray = numpy.linspace(xlow, xhi, mfile * 5)
     ylow = 2; yhi = 5
@@ -499,27 +499,26 @@ def integral_topotool(func, mfile, funcflag):
     # Set patch data
     patch_x = numpy.linspace(xlow + 0.1, xhi - 0.1, 5)
     patch_y = numpy.linspace(ylow + 0.1, yhi - 0.1, 4)
-    patch_dx = patch_x[1] - patch_x[0]
-    patch_dy = patch_y[1] - patch_y[0]
-    patch1 = topotools.patch(patch_x, patch_y, patch_dx, patch_dy)
+    patch = {'x':patch_x, 'y':patch_y}
     
     # Accurate cell value
-    real_value = numpy.empty((len(patch1.y) - 1, len(patch1.x) - 1))
-    for i in range(len(patch1.y) - 1):
-        for j in range(len(patch1.x) - 1):
-            area = (patch1.x[j+1] - patch1.x[j]) * (patch1.y[i+1] - patch1.y[i])
+    real_value = numpy.empty((len(patch['y']) - 1, len(patch['x']) - 1))
+    for i in range(len(patch['y']) - 1):
+        for j in range(len(patch['x']) - 1):
+            area = (patch['x'][j+1] - patch['x'][j]) * (patch['y'][i+1] - patch['y'][i])
             if funcflag:
                 options = {'limit':1000}
-                real_value[i][j] = nquad(func=func, ranges=[[patch1.x[j], patch1.x[j+1]],
-                                         [patch1.y[i], patch1.y[i+1]]], args=None,
+                real_value[i][j] = nquad(func=func, ranges=[[patch['x'][j], patch['x'][j+1]],
+                                         [patch['y'][i], patch['y'][i+1]]], args=None,
                                          opts=[options,options])[0] / float(area)
             else:
-                real_value[i][j] = dblquad(func=func, a=patch1.y[i], b=patch1.y[i+1],
-                                           gfun=patch1.x[j],
-                                           hfun=patch1.x[j+1])[0] / float(area)
+                real_value[i][j] = dblquad(func=func, a=patch['y'][i], b=patch['y'][i+1],
+                                           gfun=patch['x'][j],
+                                           hfun=patch['x'][j+1])[0] / float(area)
+        
     
     # Cell value calculated by functions
-    calculated_value = topotools.cell_average_patch(patch1, mtopoorder, mfile, topo)
+    calculated_value = topotools.cell_average_patch(patch, mtopoorder, mfile, topo)
         
     # Whether calculated value achieve the expected resolution
     npt.assert_almost_equal(real_value, calculated_value, decimal=3)

From 74740222c2cb4ad38b673fbe530bda046247e9f7 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Mon, 11 Jan 2021 00:31:15 -0500
Subject: [PATCH 21/22] Update topotools.py

---
 src/python/geoclaw/topotools.py | 73 ++++++++++++++-------------------
 1 file changed, 31 insertions(+), 42 deletions(-)

diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 5251ce5ab..3f30e40a1 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -1600,11 +1600,9 @@ def intersection(rect1, rect2):
     r"""
     Whether two domains, rect1 and rect2 intersect. If they intersect,
     find the intersection's boundary and area.
-
     :Input:
      - *rect1* (ndarray(:))
      - *rect2* (ndarray(:))
-
     :Output:
      - *mark* (bool) If rect1 and rect2 intersect, mark = True. Otherwise,
        mark = False.
@@ -1657,7 +1655,6 @@ def bilinearintegral(domain, topodomain, corners):
      - *domain* (ndarray(:))
      - *topodomain* (ndarray(:))
      - *corners* (ndarray(:, :)) Four corner points of the topodomain.
-
     :Output:
      - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
        the intersection of domain and topodomain.
@@ -1665,20 +1662,20 @@ def bilinearintegral(domain, topodomain, corners):
 
     # Set boundary for integral
     bound = intersection(domain, topodomain)[2]
-    
+
     # Find terms for the integration
     deltax = topodomain[1] - topodomain[0]
     deltay = topodomain[3] - topodomain[2]
     area = (bound[3] - bound[2]) * (bound[1] - bound[0])
     sumxi = (bound[1] + bound[0] - 2.0 * topodomain[0]) / deltax
     sumeta = (bound[3] + bound[2] - 2.0 * topodomain[2]) / deltay
-    
+
     # Find coefficients of the function, f(xi, eta)
     a = corners[1][0] - corners[0][0]
     b = corners[0][1] - corners[0][0]
     c = corners[1][1] - corners[1][0] - corners[0][1] + corners[0][0]
     d = corners[0][0]
-    
+
     bilinearintegral = (0.5 * (a * sumxi + b * sumeta) +
                         0.25 * c * sumxi * sumeta + d) * area
 
@@ -1704,7 +1701,6 @@ def bilinearintegral_s(domain, topodomain, corners):
      - *domain* (ndarray(:))
      - *topodomain* (ndarray(:))
      - *corners* (ndarray(:, :)) Four corner points of the topodomain.
-
     :Output:
      - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
        the intersection of domain and topodomain.
@@ -1740,7 +1736,7 @@ def bilinearintegral_s(domain, topodomain, corners):
                           r2d * (c * xdiff2 + b * intdx) * cbsinint) * (Rearth * d2r)**2
 
     return bilinearintegral_s
-        
+
 
 def topointegral(domain, topoparam, z, intmethod, coordinate_system):
     r"""
@@ -1756,78 +1752,77 @@ def topointegral(domain, topoparam, z, intmethod, coordinate_system):
      - *z* (ndarray(:, :)) Topo data of domain represented by topoparam.
      - *intmethod* (int) The method of integral.
      - *coordinate_system* (int) The type of coordinate system.
-
     :Output:
      - *theintegral* (float) The integral of a surface over a rectangular
        region which is the intersection of the domain and "topoparam".
     """
-    
+
     theintegral = 0.0
-    
+
     # Set topo's parameter
     x1 = topoparam[0]; x2 = topoparam[1]
     y1 = topoparam[2]; y2 = topoparam[3]
     topodx = topoparam[4]; topody = topoparam[5]
-    topomx = numpy.array(z).shape[0]
-    topomy = numpy.array(z).shape[1]
-    
+    topomx = numpy.array(z).shape[1]
+    topomy = numpy.array(z).shape[0]
+
     # Find the intersection of the domain and topo
     bound = intersection(domain, topoparam)[2]
     xlow = bound[0]; xhi = bound[1]
     ylow = bound[2]; yhi = bound[3]
-    
+
     if intmethod == 1:
-        
+
         # Find topo's start and end points for integral
         # The x coodinate of the topo's start point
         istart = -1
         for i in range(topomx):
             if (x1 + i * topodx) > xlow:
-                istart = i
+                istart = max(i-1, 0)
                 break
-        
+
         # The y coodinate of the topo's start point
         jstart = -1
         for j in range(topomy):
             if (y1 + j * topody)  > ylow:
-                jstart = j
+                jstart = max(j-1, 0)
                 break
-        
+ 
         # The x coodinate of the topo's end point
         iend = topomx
         for m in range(topomx):
             if (x1 + m * topodx) >= xhi:
                 iend = m
                 break
-        
+   
         # The y coodinate of the topo's end point
         jend = topomy
         for n in range(topomy):
             if (y1 + n * topody) >= yhi:
                 jend = n
                 break
-                        
+
         # Prevent overflow
-        jstart = max(jstart, 1)
-        istart = max(istart, 1)
+        jstart = max(jstart, 0)
+        istart = max(istart, 0)
         jend = min(jend, topomy - 1)
         iend = min(iend, topomx - 1)
-        
+
         # Topo integral
-        for jj in range(jstart, jend + 1):
-            yint1 = y1 + (jj - 1.0) * topody
-            yint2 = y1 + jj * topody
+        for jj in range(jstart, jend):
+            yint1 = y1 + jj * topody
+            yint2 = y1 + (jj + 1) * topody
 
-            for ii in range(istart, iend + 1):
-                xint1 = x1 + (ii - 1.0) * topodx
-                xint2 = x1 + ii * topodx
+            for ii in range(istart, iend):
+                xint1 = x1 + ii * topodx
+                xint2 = x1 + (ii + 1) * topodx
 
                 # Four corners of topo's grid
                 corners = numpy.empty((2, 2))
-                corners[0][0] = z[ii-1][topomy - jj]
-                corners[0][1] = z[ii-1][topomy - 1 - jj]
-                corners[1][0] = z[ii][topomy - jj]
-                corners[1][1] = z[ii][topomy - 1 - jj]
+                corners[0][0] = z[jj][ii]
+                corners[0][1] = z[jj+1][ii]
+                corners[1][0] = z[jj][ii+1]
+                corners[1][1] = z[jj+1][ii+1]
 
                 # Parameters of topo's grid integral
                 topointparam = [xint1, xint2, yint1, yint2]
@@ -1837,7 +1832,7 @@ def topointegral(domain, topoparam, z, intmethod, coordinate_system):
                 elif coordinate_system == 2:
                     theintegral += bilinearintegral_s(domain, topointparam, corners)
                 else:
-                     print("TOPOINTEGRAL: coordinate_system error")
+                    print("TOPOINTEGRAL: coordinate_system error")
 
     else:
         print("TOPOINTEGRAL: only intmethod = 1,2 is supported")
@@ -1850,11 +1845,9 @@ def recurintegral(domain, mtopoorder, mtopofiles, m, topo):
     Compute the integral of topo over the rectangle domain defined by "domain".
     Find the index of the topo whose resolution is m by "mtopoorder". Use this
     index to find corresponding topo object in the list "topo".
-
     The main call to this recursive function has corners of a grid cell for the
     rectangle and m = 1 in order to compute the integral over the cell using all
     topo objects.
-
     The recursive strategy is to first compute the integral using only the topo
     object of m+1 resolution. Then apply corrections due to adding m resolution
     topo object.
@@ -1871,7 +1864,6 @@ def recurintegral(domain, mtopoorder, mtopofiles, m, topo):
      - *mtopofiles* (int) The number of the topo objects.
      - *m* (int) Resolution of the topo objects.
      - *topo* (list) The list of topo objects.
-
     :Output:
      - *integral* (float) Integral.
     """
@@ -1918,7 +1910,6 @@ def cellintegral(cell, mtopoorder, mtopofiles, topo):
     - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
     - *mtopofiles* (int) The number of the topo objects.
     - *topo* (list) The list of topo objects.
-
     :Output:
     - *topoint* (float) The integral on the cell, "cell".
     """
@@ -1956,13 +1947,11 @@ def cellintegral(cell, mtopoorder, mtopofiles, topo):
 def cell_average_patch(patch, mtopoorder, mtopofiles, topo):
     r"""
     Compute every cell's average of the specific patch given by "patch".
-
     :Input:
     - *patch* (Dictionary) Patch.
     - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
     - *mtopofiles* (int) The number of the topo objects.
     - *topo* (list) The list of topo objects.
-
     :Output:
     - *cell_average* (ndarray(:, :)) Patch's every cell average.
     """

From 7254e493eb8cbd3c654d665d0113754e01f9dfa7 Mon Sep 17 00:00:00 2001
From: Runxin Ni <45889179+RunxinNi@users.noreply.github.com>
Date: Mon, 11 Jan 2021 00:48:30 -0500
Subject: [PATCH 22/22] Update test_topotools.py

---
 tests/test_topotools.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index 036bc6029..f043cad64 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -489,9 +489,8 @@ def integral_topotool(func, mfile, funcflag):
                     z[m][n] = func(x[m], y[n])
             topo1.Z = z
         else:
-            topo1.x = numpy.flip(y)
-            topo1.y = x
-            topo1.Z = func(topo1.Y, topo1.X)
+            topo_x, topo_y = numpy.meshgrid(x, y)
+            topo1.Z = func(topo_x, topo_y)
         topo1.x = x; topo1.y = y
         mtopoorder.append(mfile - 1 - i)
         topo.append(topo1)