Utilities

Matrix Utilities

discretize.utils.matutils.mkvc(x, numDims=1)[source]

Creates a vector with the number of dimension specified

e.g.:

a = np.array([1, 2, 3])

mkvc(a, 1).shape
    > (3, )

mkvc(a, 2).shape
    > (3, 1)

mkvc(a, 3).shape
    > (3, 1, 1)
discretize.utils.matutils.sdiag(h)[source]

Sparse diagonal matrix

discretize.utils.matutils.sdInv(M)[source]

Inverse of a sparse diagonal matrix

discretize.utils.matutils.speye(n)[source]

Sparse identity

discretize.utils.matutils.kron3(A, B, C)[source]

Three kron prods

discretize.utils.matutils.spzeros(n1, n2)[source]

a sparse matrix of zeros

discretize.utils.matutils.ddx(n)[source]

Define 1D derivatives, inner, this means we go from n+1 to n

discretize.utils.matutils.av(n)[source]

Define 1D averaging operator from nodes to cell-centers.

discretize.utils.matutils.av_extrap(n)[source]

Define 1D averaging operator from cell-centers to nodes.

discretize.utils.matutils.ndgrid(*args, **kwargs)[source]

Form tensorial grid for 1, 2, or 3 dimensions.

Returns as column vectors by default.

To return as matrix input:

ndgrid(..., vector=False)

The inputs can be a list or separate arguments.

e.g.:

a = np.array([1, 2, 3])
b = np.array([1, 2])

XY = ndgrid(a, b)
    > [[1 1]
       [2 1]
       [3 1]
       [1 2]
       [2 2]
       [3 2]]

X, Y = ndgrid(a, b, vector=False)
    > X = [[1 1]
           [2 2]
           [3 3]]
    > Y = [[1 2]
           [1 2]
           [1 2]]
discretize.utils.matutils.ind2sub(shape, inds)[source]

From the given shape, returns the subscripts of the given index

discretize.utils.matutils.sub2ind(shape, subs)[source]

From the given shape, returns the index of the given subscript

discretize.utils.matutils.getSubArray(A, ind)[source]

subArray

discretize.utils.matutils.inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33, returnMatrix=True)[source]

B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33)

inverts a stack of 3x3 matrices

Input:
A - a11, a12, a13, a21, a22, a23, a31, a32, a33
Output:
B - inverse
discretize.utils.matutils.inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True)[source]

B = inv2X2BlockDiagonal(a11, a12, a21, a22)

Inverts a stack of 2x2 matrices by using the inversion formula

inv(A) = (1/det(A)) * cof(A)^T

Input: A - a11, a12, a21, a22

Output: B - inverse

class discretize.utils.matutils.TensorType(M, tensor)[source]
discretize.utils.matutils.makePropertyTensor(M, tensor)[source]
discretize.utils.matutils.invPropertyTensor(M, tensor, returnMatrix=False)[source]
class discretize.utils.matutils.Zero[source]
transpose()[source]
T
class discretize.utils.matutils.Identity(positive=True)[source]
T
transpose()[source]

Curv Utilities

discretize.utils.curvutils.volTetra(xyz, A, B, C, D)[source]

Returns the volume for tetrahedras volume specified by the indexes A to D.

Parameters:
Return type:

numpy.ndarray

Returns:

V, volume of the tetrahedra

Algorithm https://en.wikipedia.org/wiki/Tetrahedron#Volume

\[ \begin{align}\begin{aligned}V = {1 \over 3} A h\\V = {1 \over 6} | ( a - d ) \cdot ( ( b - d ) ( c - d ) ) |\end{aligned}\end{align} \]
discretize.utils.curvutils.indexCube(nodes, gridSize, n=None)[source]

Returns the index of nodes on the mesh.

Input:
nodes - string of which nodes to return. e.g. ‘ABCD’ gridSize - size of the nodal grid n - number of nodes each i,j,k direction: [ni,nj,nk]
Output:
index - index in the order asked e.g. ‘ABCD’ –> (A,B,C,D)

TWO DIMENSIONS:

node(i,j)          node(i,j+1)
     A -------------- B
     |                |
     |    cell(i,j)   |
     |        I       |
     |                |
    D -------------- C
node(i+1,j)        node(i+1,j+1)

THREE DIMENSIONS:

      node(i,j,k+1)       node(i,j+1,k+1)
          E --------------- F
         /|               / |
        / |              /  |
       /  |             /   |
node(i,j,k)         node(i,j+1,k)
     A -------------- B     |
     |    H ----------|---- G
     |   /cell(i,j)   |   /
     |  /     I       |  /
     | /              | /
     D -------------- C
node(i+1,j,k)      node(i+1,j+1,k)
discretize.utils.curvutils.faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True)[source]

function [N] = faceInfo(y,A,B,C,D)

Returns the averaged normal, area, and edge lengths for a given set of faces.

If average option is FALSE then N is a cell array {nA,nB,nC,nD}

Input:
xyz - X,Y,Z vertex vector A,B,C,D - vert index of the face (counter clockwize)
Options:
average - [true]/false, toggles returning all normals or the average
Output:
N - average face normal or {nA,nB,nC,nD} if average = false area - average face area edgeLengths - exact edge Lengths, 4 column vector [AB, BC, CD, DA]

see also testFaceNormal testFaceArea

@author Rowan Cockett

Last modified on: 2013/07/26

Mesh Utilities

discretize.utils.meshutils.exampleLrmGrid(nC, exType)[source]
discretize.utils.meshutils.random_model(shape, seed=None, anisotropy=None, its=100, bounds=None)[source]

Create a random model by convolving a kernel with a uniformly distributed model.

Parameters:
  • shape (tuple) – shape of the model.
  • seed (int) – pick which model to produce, prints the seed if you don’t choose.
  • anisotropy (numpy.ndarray) – this is the (3 x n) blurring kernel that is used.
  • its (int) – number of smoothing iterations
  • bounds (list) – bounds on the model, len(list) == 2
Return type:

numpy.ndarray

Returns:

M, the model

(Source code, png, hires.png, pdf)

../_images/utils-1.png
discretize.utils.meshutils.meshTensor(value)[source]

meshTensor takes a list of numbers and tuples that have the form:

mT = [ float, (cellSize, numCell), (cellSize, numCell, factor) ]

For example, a time domain mesh code needs many time steps at one time:

[(1e-5, 30), (1e-4, 30), 1e-3]

Means take 30 steps at 1e-5 and then 30 more at 1e-4, and then one step of 1e-3.

Tensor meshes can also be created by increase factors:

[(10.0, 5, -1.3), (10.0, 50), (10.0, 5, 1.3)]

When there is a third number in the tuple, it refers to the increase factor, if this number is negative this section of the tensor is flipped right-to-left.

(Source code, png, hires.png, pdf)

../_images/utils-2.png
discretize.utils.meshutils.closestPoints(mesh, pts, gridLoc='CC')[source]

Move a list of points to the closest points on a grid.

Parameters:
  • mesh (BaseMesh) – The mesh
  • pts (numpy.ndarray) – Points to move
  • gridLoc (str) – [‘CC’, ‘N’, ‘Fx’, ‘Fy’, ‘Fz’, ‘Ex’, ‘Ex’, ‘Ey’, ‘Ez’]
Return type:

numpy.ndarray

Returns:

nodeInds

discretize.utils.meshutils.ExtractCoreMesh(xyzlim, mesh, meshType='tensor')[source]

Extracts Core Mesh from Global mesh

Parameters:

This function ouputs:

- actind: corresponding boolean index from global to core
- meshcore: core sdiscretize mesh

Interpolation Utilities

discretize.utils.interputils.interpmat(locs, x, y=None, z=None)[source]

Local interpolation computed for each receiver point in turn

Parameters:
  • loc (numpy.ndarray) – Location of points to interpolate to
  • x (numpy.ndarray) – Tensor of 1st dimension of grid.
  • y (numpy.ndarray) – Tensor of 2nd dimension of grid. None by default.
  • z (numpy.ndarray) – Tensor of 3rd dimension of grid. None by default.
Return type:

scipy.sparse.csr_matrix

Returns:

Interpolation matrix

(Source code, png, hires.png, pdf)

../_images/utils-3.png

Testing

discretize.Tests.setupMesh(meshType, nC, nDim)[source]

For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc.

class discretize.Tests.OrderTest(methodName='runTest')[source]

OrderTest is a base class for testing convergence orders with respect to mesh sizes of integral/differential operators.

Mathematical Problem:

Given are an operator A and its discretization A[h]. For a given test function f and h –> 0 we compare:

\[error(h) = \| A[h](f) - A(f) \|_{\infty}\]

Note that you can provide any norm.

Test is passed when estimated rate order of convergence is at least within the specified tolerance of the estimated rate supplied by the user.

Minimal example for a curl operator:

class TestCURL(OrderTest):
    name = "Curl"

    def getError(self):
        # For given Mesh, generate A[h], f and A(f) and return norm of error.


        fun  = lambda x: np.cos(x)  # i (cos(y)) + j (cos(z)) + k (cos(x))
        sol = lambda x: np.sin(x)  # i (sin(z)) + j (sin(x)) + k (sin(y))


        Ex = fun(self.M.gridEx[:, 1])
        Ey = fun(self.M.gridEy[:, 2])
        Ez = fun(self.M.gridEz[:, 0])
        f = np.concatenate((Ex, Ey, Ez))

        Fx = sol(self.M.gridFx[:, 2])
        Fy = sol(self.M.gridFy[:, 0])
        Fz = sol(self.M.gridFz[:, 1])
        Af = np.concatenate((Fx, Fy, Fz))

        # Generate DIV matrix
        Ah = self.M.edgeCurl

        curlE = Ah*E
        err = np.linalg.norm((Ah*f -Af), np.inf)
        return err

    def test_order(self):
        # runs the test
        self.orderTest()

See also: test_operatorOrder.py

name = 'Order Test'
expectedOrders = 2.0
tolerance = 0.85
meshSizes = [4, 8, 16, 32]
meshTypes = ['uniformTensorMesh']
meshDimension = 3
setupMesh(nC)[source]
getError()[source]

For given h, generate A[h], f and A(f) and return norm of error.

orderTest()[source]

For number of cells specified in meshSizes setup mesh, call getError and prints mesh size, error, ratio between current and previous error, and estimated order of convergence.

discretize.Tests.Rosenbrock(x, return_g=True, return_H=True)[source]

Rosenbrock function for testing GaussNewton scheme

discretize.Tests.checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.85, eps=1e-10, ax=None)[source]

Basic derivative check

Compares error decay of 0th and 1st order Taylor approximation at point x0 for a randomized search direction.

Parameters:
  • fctn (callable) – function handle
  • x0 (numpy.ndarray) – point at which to check derivative
  • num (int) – number of times to reduce step length, h
  • plotIt (bool) – if you would like to plot
  • dx (numpy.ndarray) – step direction
  • expectedOrder (int) – The order that you expect the derivative to yield.
  • tolerance (float) – The tolerance on the expected order.
  • eps (float) – What is zero?
Return type:

bool

Returns:

did you pass the test?!

from discretize import Tests, utils
import numpy as np

def simplePass(x):
    return np.sin(x), utils.sdiag(np.cos(x))
Tests.checkDerivative(simplePass, np.random.randn(5))

(Source code, png, hires.png, pdf)

../_images/utils-4.png
discretize.Tests.getQuadratic(A, b, c=0)[source]

Given A, b and c, this returns a quadratic, Q

\[\mathbf{Q( x ) = 0.5 x A x + b x} + c\]