Python numerix.MA类代码示例

    def _calcFaceToCellDistances(self):
        tmp = MA.repeat(self.faceCenters[...,numerix.NewAxis,:], 2, 1)
        # array -= masked_array screws up masking for on numpy 1.1

        tmp = tmp - numerix.take(self.cellCenters, self.faceCellIDs, axis=1)
        self.cellToFaceDistanceVectors = tmp
        self.faceToCellDistances = MA.sqrt(MA.sum(tmp * tmp,0))

 def deleteIslands(self):
     cellToCellIDs = self.mesh._cellToCellIDs
     adjVals = numerix.take(self.value, cellToCellIDs)
     adjInterfaceValues = MA.masked_array(adjVals, mask = (adjVals * self.value) > 0)
     masksum = numerix.sum(numerix.logical_not(MA.getmask(adjInterfaceValues)), 0)
     tmp = MA.logical_and(masksum == 4, self.value > 0)
     self.value = numerix.array(MA.where(tmp, -1, self.value))

    def _buildMatrix(self, var, SparseMatrix, boundaryConditions=(), dt=None, equation=None, transientGeomCoeff=None, diffusionGeomCoeff=None):

        oldArray = var.old

        mesh = var.mesh
        NCells = mesh.numberOfCells
        NCellFaces = mesh._maxFacesPerCell

        cellValues = numerix.repeat(oldArray[numerix.newaxis, ...], NCellFaces, axis = 0)

        cellIDs = numerix.repeat(numerix.arange(NCells)[numerix.newaxis, ...], NCellFaces, axis = 0)
        cellToCellIDs = mesh._cellToCellIDs

        if NCells > 0:
            cellToCellIDs = MA.where(MA.getmask(cellToCellIDs), cellIDs, cellToCellIDs)

            adjacentValues = numerix.take(oldArray, cellToCellIDs)

            differences = self._getDifferences(adjacentValues, cellValues, oldArray, cellToCellIDs, mesh)
            differences = MA.filled(differences, 0)

            minsq = numerix.sqrt(numerix.sum(numerix.minimum(differences, numerix.zeros((NCellFaces, NCells), 'l'))**2, axis=0))
            maxsq = numerix.sqrt(numerix.sum(numerix.maximum(differences, numerix.zeros((NCellFaces, NCells), 'l'))**2, axis=0))

            coeff = numerix.array(self._getGeomCoeff(var))

            coeffXdifferences = coeff * ((coeff > 0.) * minsq + (coeff < 0.) * maxsq)
        else:
            coeffXdifferences = 0.

        return (var, SparseMatrix(mesh=var.mesh), -coeffXdifferences * mesh.cellVolumes)

    def _calcValue(self):

        Nfaces = self.mesh.numberOfFaces
        M = self.mesh._maxFacesPerCell
        dim = self.mesh.dim
        cellFaceIDs = self.mesh.cellFaceIDs

        faceNormalAreas = self.distanceVar._levelSetNormals * self.mesh._faceAreas

        cellFaceNormalAreas = numerix.array(MA.filled(numerix.take(faceNormalAreas, cellFaceIDs, axis=-1), 0))
        norms = numerix.array(MA.filled(MA.array(self.mesh._cellNormals), 0))

        alpha = numerix.dot(cellFaceNormalAreas, norms)
        alpha = numerix.where(alpha > 0, alpha, 0)

        alphasum = numerix.sum(alpha, axis=0)
        alphasum += (alphasum < 1e-100) * 1.0
        alpha = alpha / alphasum

        phi = numerix.repeat(self.distanceVar[numerix.newaxis, ...], M, axis=0)
        alpha = numerix.where(phi > 0., 0, alpha)

        volumes = numerix.array(self.mesh.cellVolumes)
        alpha = alpha * volumes * norms

        value = numerix.zeros((dim, Nfaces), 'd')

        vector._putAdd(value, cellFaceIDs, alpha, mask=MA.getmask(MA.array(cellFaceIDs)))

##         value = numerix.reshape(value, (dim, Nfaces, dim))

        return -value / self.mesh._faceAreas

    def faceCellIDs(self):
        XYids = MA.zeros((2, self.nx, self.ny, self.nz + 1), 'l')
        indices = numerix.indices((self.nx, self.ny, self.nz + 1))
        XYids[1] = indices[0] + (indices[1] + indices[2] * self.ny) * self.nx
        XYids[0] = XYids[1] - self.nx * self.ny
        XYids[0, ..., 0] = XYids[1, ..., 0]
        XYids[1, ..., 0] = MA.masked
        XYids[1, ..., -1] = MA.masked

        XZids = MA.zeros((2, self.nx, self.ny + 1, self.nz), 'l')
        indices = numerix.indices((self.nx, self.ny + 1, self.nz))
        XZids[1] = indices[0] + (indices[1] + indices[2] * self.ny) * self.nx
        XZids[0] = XZids[1] - self.nx
        XZids[0,:, 0,:] = XZids[1,:, 0,:]
        XZids[1,:, 0,:] = MA.masked
        XZids[1,:, -1,:] = MA.masked

        YZids = MA.zeros((2, self.nx + 1, self.ny, self.nz), 'l')
        indices = numerix.indices((self.nx + 1, self.ny, self.nz))
        YZids[1] = indices[0] + (indices[1] + indices[2] * self.ny) * self.nx
        YZids[0] = YZids[1] - 1
        YZids[0, 0] = YZids[1, 0]
        YZids[1, 0] = MA.masked
        YZids[1, -1] = MA.masked

        return MA.concatenate((XYids.swapaxes(1, 3).reshape((2, self.numberOfXYFaces)),
                               XZids.swapaxes(1, 3).reshape((2, self.numberOfXZFaces)),
                               YZids.swapaxes(1, 3).reshape((2, self.numberOfYZFaces))), axis=1)

    def __init__(self, mesh, name = '', value = 0., unit = None, hasOld = 0, narrowBandWidth = 1e+10):
        """
        Creates a `distanceVariable` object.

        :Parameters:
          - `mesh`: The mesh that defines the geometry of this variable.
          - `name`: The name of the variable.
	  - `value`: The initial value.
	  - `unit`: the physical units of the variable
          - `hasOld`: Whether the variable maintains an old value.
          - `narrowBandWidth`: The width of the region about the zero level set
            within which the distance function is evaluated.

        """
        CellVariable.__init__(self, mesh, name = name, value = value, unit = unit, hasOld = hasOld)
        self._markStale()
        self.narrowBandWidth = narrowBandWidth

        self.cellToCellDistances = MA.filled(self.mesh._getCellToCellDistances(), 0)
        self.cellNormals = MA.filled(self.mesh._getCellNormals(), 0)      
        self.cellAreas = MA.filled(self.mesh._getCellAreas(), 0)
##         self.cellToCellDistances = numerix.array(MA.array(self.mesh._getCellToCellDistances()).filled(0))
##         self.cellNormals = numerix.array(MA.array(self.mesh._getCellNormals()).filled(0))       
##         self.cellAreas = numerix.array(MA.array(self.mesh._getCellAreas()).filled(0))
        self.cellToCellIDs = numerix.array(self.mesh._getCellToCellIDsFilled())
        self.adjacentCellIDs = self.mesh._getAdjacentCellIDs()
        self.exteriorFaces = self.mesh.getExteriorFaces()
        self.cellFaceIDs = self.mesh._getCellFaceIDs()

    def _buildMatrix(self, var, SparseMatrix, boundaryCondtions=(), dt=None, equation=None):

        oldArray = var.getOld()

        mesh = var.getMesh()
        NCells = mesh.getNumberOfCells()
        NCellFaces = mesh._getMaxFacesPerCell()

        cellValues = numerix.repeat(oldArray[numerix.newaxis, ...], NCellFaces, axis = 0)
        
        cellIDs = numerix.repeat(numerix.arange(NCells)[numerix.newaxis, ...], NCellFaces, axis = 0)
        cellToCellIDs = mesh._getCellToCellIDs()

        if NCells > 0:
            cellToCellIDs = MA.where(MA.getmask(cellToCellIDs), cellIDs, cellToCellIDs) 

            adjacentValues = numerix.take(oldArray, cellToCellIDs)

            differences = self._getDifferences(adjacentValues, cellValues, oldArray, cellToCellIDs, mesh)
            differences = MA.filled(differences, 0)
            
            minsq = numerix.sqrt(numerix.sum(numerix.minimum(differences, numerix.zeros((NCellFaces, NCells)))**2, axis=0))
            maxsq = numerix.sqrt(numerix.sum(numerix.maximum(differences, numerix.zeros((NCellFaces, NCells)))**2, axis=0))

            coeff = numerix.array(self._getGeomCoeff(mesh))

            coeffXdiffereneces = coeff * ((coeff > 0.) * minsq + (coeff < 0.) * maxsq)
        else:
            coeffXdiffereneces = 0.

        return (SparseMatrix(mesh=var.getMesh()), -coeffXdiffereneces * mesh.getCellVolumes())

    def _getDifferences(self, adjacentValues, cellValues, oldArray, cellToCellIDs, mesh):

        dAP = mesh._cellToCellDistances

##        adjacentGradient = numerix.take(oldArray.grad, cellToCellIDs)
        adjacentGradient = numerix.take(oldArray.grad, mesh._cellToCellIDs, axis=-1)
        adjacentNormalGradient = numerix.dot(adjacentGradient, mesh._cellNormals)
        adjacentUpValues = cellValues + 2 * dAP * adjacentNormalGradient

        cellIDs = numerix.repeat(numerix.arange(mesh.numberOfCells)[numerix.newaxis, ...],
                mesh._maxFacesPerCell, axis=0)
        cellIDs = MA.masked_array(cellIDs, mask = MA.getmask(mesh._cellToCellIDs))
        cellGradient = numerix.take(oldArray.grad, cellIDs, axis=-1)
        cellNormalGradient = numerix.dot(cellGradient, mesh._cellNormals)
        cellUpValues = adjacentValues - 2 * dAP * cellNormalGradient

        cellLaplacian = (cellUpValues + adjacentValues - 2 * cellValues) / dAP**2

        adjacentLaplacian = (adjacentUpValues + cellValues - 2 * adjacentValues) / dAP**2
        adjacentLaplacian = adjacentLaplacian.filled(0)
        cellLaplacian = cellLaplacian.filled(0)

        mm = numerix.where(cellLaplacian * adjacentLaplacian < 0.,
                           0.,
                           numerix.where(abs(cellLaplacian) > abs(adjacentLaplacian),
                                         adjacentLaplacian,
                                         cellLaplacian))

        return FirstOrderAdvectionTerm._getDifferences(self, adjacentValues, cellValues, oldArray, cellToCellIDs, mesh) -  mm * dAP / 2.

    def __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs, communicator=serialComm, _RepresentationClass=_MeshRepresentation, _TopologyClass=_MeshTopology):
        super(Mesh, self).__init__(communicator=communicator,
                                   _RepresentationClass=_RepresentationClass,
                                   _TopologyClass=_TopologyClass)

        """faceVertexIds and cellFacesIds must be padded with minus ones."""
                                   
        self.vertexCoords = vertexCoords
        self.faceVertexIDs = MA.masked_values(faceVertexIDs, -1)
        self.cellFaceIDs = MA.masked_values(cellFaceIDs, -1)

        self.dim = self.vertexCoords.shape[0]

        if not hasattr(self, "numberOfFaces"):
            self.numberOfFaces = self.faceVertexIDs.shape[-1]
        if not hasattr(self, "numberOfCells"):
            self.numberOfCells = self.cellFaceIDs.shape[-1]
        if not hasattr(self, "globalNumberOfCells"):
            self.globalNumberOfCells = self.numberOfCells
        if not hasattr(self, "globalNumberOfFaces"):
            self.globalNumberOfFaces = self.numberOfFaces

        self.faceCellIDs = self._calcFaceCellIDs() 

        self._setTopology()
        self._setGeometry(scaleLength = 1.)

    def __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs):
        """faceVertexIds and cellFacesIds must be padded with minus ones."""

        self.vertexCoords = vertexCoords
        self.faceVertexIDs = MA.masked_values(faceVertexIDs, -1)
        self.cellFaceIDs = MA.masked_values(cellFaceIDs, -1)

        _CommonMesh.__init__(self)

    def getCellInterfaceAreas(self):
        """
        Returns the length of the interface that crosses the cell

        A simple 1D test:

        >>> from fipy.meshes.grid1D import Grid1D
        >>> mesh = Grid1D(dx = 1., nx = 4)
        >>> distanceVariable = DistanceVariable(mesh = mesh, 
        ...                                     value = (-1.5, -0.5, 0.5, 1.5))
        >>> answer = CellVariable(mesh=mesh, value=(0, 0., 1., 0))
        >>> print numerix.allclose(distanceVariable.getCellInterfaceAreas(), 
        ...                        answer)
        True

        A 2D test case:
        
        >>> from fipy.meshes.grid2D import Grid2D
        >>> from fipy.variables.cellVariable import CellVariable
        >>> mesh = Grid2D(dx = 1., dy = 1., nx = 3, ny = 3)
        >>> distanceVariable = DistanceVariable(mesh = mesh, 
        ...                                     value = (1.5, 0.5, 1.5,
        ...                                              0.5,-0.5, 0.5,
        ...                                              1.5, 0.5, 1.5))
        >>> answer = CellVariable(mesh=mesh,
        ...                       value=(0, 1, 0, 1, 0, 1, 0, 1, 0))
        >>> print numerix.allclose(distanceVariable.getCellInterfaceAreas(), answer)
        True

        Another 2D test case:

        >>> mesh = Grid2D(dx = .5, dy = .5, nx = 2, ny = 2)
        >>> from fipy.variables.cellVariable import CellVariable
        >>> distanceVariable = DistanceVariable(mesh = mesh, 
        ...                                     value = (-0.5, 0.5, 0.5, 1.5))
        >>> answer = CellVariable(mesh=mesh,
        ...                       value=(0, numerix.sqrt(2) / 4,  numerix.sqrt(2) / 4, 0))
        >>> print numerix.allclose(distanceVariable.getCellInterfaceAreas(), 
        ...                        answer)
        True

        Test to check that the circumfrence of a circle is, in fact, 
        :math:`2\pi r`.
	
        >>> mesh = Grid2D(dx = 0.05, dy = 0.05, nx = 20, ny = 20)
        >>> r = 0.25
        >>> x, y = mesh.getCellCenters()
        >>> rad = numerix.sqrt((x - .5)**2 + (y - .5)**2) - r
        >>> distanceVariable = DistanceVariable(mesh = mesh, value = rad)
        >>> print distanceVariable.getCellInterfaceAreas().sum()
        1.57984690073
        """        
        normals = numerix.array(MA.filled(self._getCellInterfaceNormals(), 0))
        areas = numerix.array(MA.filled(self.mesh._getCellAreaProjections(), 0))
        return CellVariable(mesh=self.mesh, 
                            value=numerix.sum(abs(numerix.dot(normals, areas)), axis=0))

 def _calcFaceAreas(self):
     faceVertexIDs = MA.filled(self.faceVertexIDs, -1)
     substitute = numerix.repeat(faceVertexIDs[numerix.newaxis, 0], 
                                 faceVertexIDs.shape[0], axis=0)
     faceVertexIDs = numerix.where(MA.getmaskarray(self.faceVertexIDs), substitute, faceVertexIDs)
     faceVertexCoords = numerix.take(self.vertexCoords, faceVertexIDs, axis=1)
     faceOrigins = numerix.repeat(faceVertexCoords[:,0], faceVertexIDs.shape[0], axis=0)
     faceOrigins = numerix.reshape(faceOrigins, MA.shape(faceVertexCoords))
     faceVertexCoords = faceVertexCoords - faceOrigins
     left = range(faceVertexIDs.shape[0])
     right = left[1:] + [left[0]]
     cross = numerix.sum(numerix.cross(faceVertexCoords, numerix.take(faceVertexCoords, right, 1), axis=0), 1)
     self.faceAreas = numerix.sqrtDot(cross, cross) / 2.

 def _cellToCellDistances(self):
     distances = MA.zeros((2, self.numberOfCells), 'd')
     distances[:] = self.dx
     if self.numberOfCells > 0:
         distances[0, 0] = self.dx / 2.
         distances[1, -1] = self.dx / 2.
     return distances

 def _calcInteriorAndExteriorFaceIDs(self):
     from fipy.variables.faceVariable import FaceVariable
     mask = MA.getmask(self.faceCellIDs[1])
     self.exteriorFaces = FaceVariable(mesh=self, 
                                       value=mask)
     self.interiorFaces = FaceVariable(mesh=self, 
                                       value=numerix.logical_not(mask))

 def _getCellToCellDistances(self):
     distances = MA.zeros((2, self.numberOfCells), "d")
     distances[:] = self.dx
     if self.numberOfCells > 0:
         distances[0, 0] = self.dx / 2.0
         distances[1, -1] = self.dx / 2.0
     return distances