Define _localNonOverlappingFaceIDs for Gmsh

fipy/matrices/trilinosMatrix.py

@@ -842,15 +842,22 @@ def _getGhostedValues(self, var):
         else:
             s = (localNonOverlappingCellIDs,)
 
+        nonOverlappingVector = var[s].ravel()
+        dtype = nonOverlappingVector.dtype
+
+        # Epetra can only handle doubles. Ugh.
+        nonOverlappingVector = nonOverlappingVector.astype(float)
+
         nonOverlappingVector = Epetra.Vector(self.domainMap,
-                                             var[s].ravel())
+                                             nonOverlappingVector)
 
         overlappingVector = Epetra.Vector(self.colMap)
         overlappingVector.Import(nonOverlappingVector,
                                  Epetra.Import(self.colMap, self.domainMap),
                                  Epetra.Insert)
 
-        return numerix.reshape(numerix.asarray(overlappingVector), var.shape)
+        return numerix.reshape(numerix.asarray(overlappingVector, dtype=dtype),
+                               var.shape)
 
     def put(self, vector, id1, id2, overlapping=False):
         """Insert local overlapping values and coordinates into global

fipy/meshes/cylindricalUniformGrid2D.py

@@ -276,7 +276,7 @@ def _test(self):
         warnings from the solver.
 
             >>> from fipy import *
-            >>> mesh = CylindricalUniformGrid2D(nx=3., ny=3., dx=1., dy=1.)
+            >>> mesh = CylindricalUniformGrid2D(nx=3, ny=3, dx=1., dy=1.)
             >>> var = CellVariable(mesh=mesh)
             >>> DiffusionTerm().solve(var, solver=DummySolver())
 

fipy/meshes/gmshMesh.py

@@ -1387,6 +1387,28 @@ def _localOverlappingCellIDs(self):
         return nx.arange(len(self.mesh.cellGlobalIDs)
                          + len(self.mesh.gCellGlobalIDs))
 
+    @property
+    def _localNonOverlappingFaceIDs(self):
+        """Return the IDs of the local mesh in isolation.
+
+        Does not include the IDs of faces of boundary cells.
+
+        E.g., would return [0, 1, 3, 4, 6, 7, 9, 10, 11, 13, 14, 15]
+        for mesh A
+
+        ```
+            A   ||   B
+        --6---7-----7---8--
+       13   14 15/14 15   16
+        --3---4-----4---5--
+        9   10 11/10 11   12
+        --0---1-----1---2--
+                ||
+        ```
+
+        .. note:: Trivial except for parallel meshes
+        """
+        return nx.arange(self.mesh.numberOfFaces)[..., self._nonOverlappingFaces]
 
 
 
@@ -2221,6 +2243,80 @@ def _test(self):
 
         >>> if parallelComm.procID == 0:
         ...     os.remove(posFile)
+
+
+        Ensure that ghost faces are excluded from accumulating operations
+        (#856).  Six exterior surfaces of :math:`10\times 10\times 10` cube
+        mesh should each have a total area of 100, regardless of
+        partitioning.
+
+        >>> geo = '''
+        ... cellSize = 1.;
+        ...
+        ... Point(1) = {0, 0, 0, cellSize};
+        ... Point(2) = {10, 0, 0, cellSize};
+        ... Point(3) = {10, 10, 0, cellSize};
+        ... Point(4) = {0, 10, 0, cellSize};
+        ...
+        ... Point(11) = {0, 0, 10, cellSize};
+        ... Point(12) = {10, 0, 10, cellSize};
+        ... Point(13) = {10, 10, 10, cellSize};
+        ... Point(14) = {0, 10, 10, cellSize};
+        ...
+        ... Line(1) = {1, 2};
+        ... Line(2) = {2, 3};
+        ... Line(3) = {3, 4};
+        ... Line(4) = {4, 1};
+        ...
+        ... Line(11) = {11, 12};
+        ... Line(12) = {12, 13};
+        ... Line(13) = {13, 14};
+        ... Line(14) = {14, 11};
+        ...
+        ... Line(21) = {1, 11};
+        ... Line(22) = {2, 12};
+        ... Line(23) = {3, 13};
+        ... Line(24) = {4, 14};
+        ...
+        ... Line Loop(1) = {1, 2, 3, 4};
+        ... Line Loop(2) = {11, 12, 13, 14};
+        ... Line Loop(3) = {1, 22, -11, -21};
+        ... Line Loop(4) = {2, 23, -12, -22};
+        ... Line Loop(5) = {3, 24, -13, -23};
+        ... Line Loop(6) = {4, 21, -14, -24};
+        ...
+        ... Plane Surface(1) = {1};
+        ... Plane Surface(2) = {2};
+        ... Plane Surface(3) = {3};
+        ... Plane Surface(4) = {4};
+        ... Plane Surface(5) = {5};
+        ... Plane Surface(6) = {6};
+        ...
+        ... Surface Loop(1) = {1, 2, 3, 4, 5, 6};
+        ...
+        ... Volume(1) = {1};
+        ...
+        ... Physical Surface("bottom") = {1};
+        ... Physical Surface("top") = {2};
+        ... Physical Surface("front") = {3};
+        ... Physical Surface("right") = {4};
+        ... Physical Surface("back") = {5};
+        ... Physical Surface("left") = {6};
+        ...
+        ... Physical Volume("box") = {1};
+        ... '''
+
+        >>> cube = Gmsh3D(geo) # doctest: +GMSH
+
+        >>> for surface in ["bottom", "top", "front", "right", "back", "left"]: # doctest: +GMSH
+        ...     area = (cube._faceAreas * cube.physicalFaces[surface]).sum() # doctest: +GMSH
+        ...     print(surface, numerix.allclose(area, 100)) # doctest: +GMSH
+        bottom True
+        top True
+        front True
+        right True
+        back True
+        left True
         """
 
 class GmshGrid2D(Gmsh2D):

fipy/meshes/nonUniformGrid1D.py

@@ -47,6 +47,9 @@ def __init__(self, dx=1., nx=None, overlap=2,
             'overlap': overlap
         }
 
+        if self.args['nx'] is None:
+            self.args['nx'] = len(self.args['dx'])
+
         builder.buildGridData([dx], [nx], overlap, communicator)
 
         ([self.dx],

fipy/meshes/nonUniformGrid2D.py

@@ -33,6 +33,12 @@ def __init__(self, dx=1., dy=1., nx=None, ny=None, overlap=2, communicator=paral
             'overlap': overlap
         }
 
+        if self.args['nx'] is None:
+            self.args['nx'] = len(self.args['dx'])
+
+        if self.args['ny'] is None:
+            self.args['ny'] = len(self.args['dy'])
+
         builder.buildGridData([dx, dy], [nx, ny], overlap, communicator)
 
         ([self.dx, self.dy],
@@ -244,6 +250,27 @@ def _test(self):
             >>> print(min(m.y) == 5.5) # doctest: +PROCESSOR_2_OF_3
             True
 
+        Ensure that ghost faces are excluded from accumulating operations
+        (#856).  Four exterior surfaces of :math:`10\times 10` square mesh
+        should each have a total area of 10, regardless of partitioning.
+
+            >>> square = NonUniformGrid2D(nx=10, dx=1., ny=10, dy=1.)
+
+            >>> area = (square._faceAreas * square.facesBottom).sum()
+            >>> print(numerix.allclose(area, 10))
+            True
+
+            >>> area = (square._faceAreas * square.facesTop).sum()
+            >>> print(numerix.allclose(area, 10))
+            True
+
+            >>> area = (square._faceAreas * square.facesLeft).sum()
+            >>> print(numerix.allclose(area, 10))
+            True
+
+            >>> area = (square._faceAreas * square.facesRight).sum()
+            >>> print(numerix.allclose(area, 10))
+            True
         """
 
 def _test():

fipy/meshes/nonUniformGrid3D.py

@@ -44,6 +44,15 @@ def __init__(self, dx = 1., dy = 1., dz = 1., nx = None, ny = None, nz = None, o
             'overlap': overlap,
         }
 
+        if self.args['nx'] is None:
+            self.args['nx'] = len(self.args['dx'])
+
+        if self.args['ny'] is None:
+            self.args['ny'] = len(self.args['dy'])
+
+        if self.args['nz'] is None:
+            self.args['nz'] = len(self.args['dz'])
+
         builder.buildGridData([dx, dy, dz], [nx, ny, nz], overlap,
                               communicator)
 
@@ -367,6 +376,36 @@ def _test(self):
             >>> print(min(m.z) == 5.5) # doctest: +PROCESSOR_2_OF_3
             True
 
+        Ensure that ghost faces are excluded from accumulating operations
+        (#856).  Four exterior surfaces of :math:`10\times 10\times 10`
+        cube mesh should each have a total area of 100, regardless of
+        partitioning.
+
+            >>> cube = NonUniformGrid3D(nx=10, dx=1., ny=10, dy=1., nz=10, dz=1.)
+
+            >>> area = (cube._faceAreas * cube.facesBottom).sum()
+            >>> print(numerix.allclose(area, 100))
+            True
+
+            >>> area = (cube._faceAreas * cube.facesTop).sum()
+            >>> print(numerix.allclose(area, 100))
+            True
+
+            >>> area = (cube._faceAreas * cube.facesLeft).sum()
+            >>> print(numerix.allclose(area, 100))
+            True
+
+            >>> area = (cube._faceAreas * cube.facesRight).sum()
+            >>> print(numerix.allclose(area, 100))
+            True
+
+            >>> area = (cube._faceAreas * cube.facesFront).sum()
+            >>> print(numerix.allclose(area, 100))
+            True
+
+            >>> area = (cube._faceAreas * cube.facesBack).sum()
+            >>> print(numerix.allclose(area, 100))
+            True
         """
 
 def _test():

fipy/meshes/topologies/abstractTopology.py

@@ -146,6 +146,87 @@ def _globalOverlappingFaceIDs(self):
         """
         return numerix.arange(self.mesh.numberOfFaces)
 
+    @property
+    def _cellProcID(self):
+        """
+        Return the processor that "owns" each cell.
+        Ownership goes to the lowest `procID` of a neighboring cell.
+        Boundary faces try not to be duplicated
+
+        E.g., would return [0, 0, 1, 0, 0, 1] for mesh A
+                       and [0, 1, 1, 0, 1, 1] for mesh B
+
+        ```
+            A        B
+        ------------------
+        | 4 | 5 || 6 | 7 |
+        ------------------
+        | 0 | 1 || 2 | 3 |
+        ------------------
+        ```
+
+        .. note:: Trivial except for parallel meshes
+        """
+        from fipy.variables.cellVariable import CellVariable
+
+        procID = CellVariable(mesh=self.mesh,
+                              value=self.mesh.communicator.procID)
+        procID._updateGhosts()
+
+        return procID
+
+    @property
+    def _ownedFaceIDs(self):
+        """
+        Return the local face IDs of the mesh "owned" by the current
+        processor.  Ownership goes to the lowest `procID` of a neighboring
+        cell.  Boundary faces try not to be duplicated
+
+        E.g., would return [0, 1, 3, 4, 6, 7, 9, 10, 11, 13, 14, 15] for mesh A
+                       and [1, 2, 4, 5, 7, 8, 11, 12, 15, 16] for mesh B
+
+        ```
+            A   ||   B
+        --6---7-----7---8--
+       13   14 15/14 15   16
+        --3---4-----4---5--
+        9   10 11/10 11   12
+        --0---1-----1---2--
+                ||
+        ```
+
+        .. note:: Trivial except for parallel meshes
+        """
+        minproc = numerix.take(self._cellProcID,
+                               self.mesh.faceCellIDs).min(axis=0)
+
+        return numerix.where(minproc == self.mesh.communicator.procID)[0]
+
+    @property
+    def _nonOverlappingFaces(self):
+        """Return mask of faces on local mesh.
+
+        False for faces only belonging to ghost cells.
+
+        E.g., would return [True, True, False, True, True, False, True,
+        True, False, True, True, True, False, True, True, True, False]
+        for mesh A
+
+        ```
+            A   ||   B
+        --6---7----8------
+       13   14  15  16   |
+        --3---4----5------
+        9   10  11  12   |
+        --0---1----2------
+                ||
+        ```
+        """
+        return (numerix.in1dMA(self.mesh.faceCellIDs[0],
+                               self._localNonOverlappingCellIDs).filled(False)
+                | numerix.in1dMA(self.mesh.faceCellIDs[1],
+                                 self._localNonOverlappingCellIDs).filled(False))
+
     @property
     def _localNonOverlappingFaceIDs(self):
         """Return the IDs of the local mesh in isolation.