Define _localNonOverlappingFaceIDs for Gmsh

fipy/matrices/trilinosMatrix.py

@@ -822,6 +822,52 @@ def copy(self):
         copy.matrix = tmp.matrix
         return copy
 
+    @staticmethod
+    def _ArrayAsEpetraVector(arr, emap):
+        """Convert numpy `ndarray` to `Epetra.Vector`
+
+        Parameters
+        ----------
+        arr : array_like
+            Input data, in any form that can be converted to an array.
+        emap : Epetra.Map
+            The map to distribute across.
+
+        Returns
+        -------
+        vec : Epetra.Vector
+            Epetra interpretation of `arr`.
+        dtype : data-type
+            The original data type of `arr`, as :class:`Epetra.Vector` can
+            only handle doubles.
+        """
+        return (Epetra.Vector(emap, numerix.asarray(arr, dtype=float)),
+                arr.dtype)
+
+    @staticmethod
+    def _EpetraVectorAsArray(vec, dtype=None, shape=None):
+        """Convert `Epetra.Vector` to numpy `ndarray`
+
+        Parameters
+        ----------
+        vec : Epetra.Vector
+            Input data.
+        dtype : data-type, optional
+            By default, the data-type is float, as :class:`Epetra.Vector`
+            can only hold doubles.
+        shape : tuple
+            Desired shape of result (default: None leaves shape along)
+
+        Returns
+        -------
+        out : ndarray
+            Array interpretation of `vec`.
+        """
+        arr = numerix.asarray(vec, dtype=dtype)
+        if shape is not None:
+            arr = numerix.reshape(arr, shape=shape)
+        return arr
+
     def _getGhostedValues(self, var):
         """Obtain current ghost values from across processes
 
@@ -830,27 +876,27 @@ def _getGhostedValues(self, var):
         ndarray
             Ghosted values
         """
-        mesh = var.mesh
-        localNonOverlappingCellIDs = mesh._localNonOverlappingCellIDs
+        s = (var.mesh._localNonOverlappingCellIDs,)
 
         ## The following conditional is required because empty indexing is
         ## not altogether functional.  This numpy.empty((0,))[[]] and this
         ## numpy.empty((0,))[...,[]] both work, but this numpy.empty((3,
         ## 0))[...,[]] is broken.
         if var.shape[-1] != 0:
-            s = (Ellipsis, localNonOverlappingCellIDs)
-        else:
-            s = (localNonOverlappingCellIDs,)
+            s = (Ellipsis,) + s
 
-        nonOverlappingVector = Epetra.Vector(self.domainMap,
-                                             var[s].ravel())
+        (nonOverlappingVector,
+         original_dtype) = self._ArrayAsEpetraVector(arr=var[s].ravel(),
+                                                     emap=self.domainMap)
 
         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 self._EpetraVectorAsArray(overlappingVector,
+                                         dtype=original_dtype,
+                                         shape=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():