pyamg.gallery#

Matrix gallery of model problems.

Functions#

  • poisson() : Poisson problem using Finite Differences

  • linear_elasticity() : Linear Elasticity using Finite Elements

  • stencil_grid() : General stencil generation from 1D, 2D, and 3D

  • diffusion_stencil_2d() : 2D rotated anisotropic FE/FD stencil

pyamg.gallery.advection_2d(grid, theta=0.7853981633974483, l_bdry=1.0, b_bdry=1.0)[source]#

Generate matrix and RHS for upwind FD discretization of 2D advection.

(cos(theta),sin(theta)) dot grad(u) = 0,

with inflow boundaries on the left and bottom of the domain. Assume uniform grid spacing, dx=dy, even for grid[0] != grid[1].

Parameters:
gridtuple

(ny, nx) number of points in y and x

thetafloat, optional

Rotation angle theta in radians defines direction of advection (cos(theta),sin(theta))

l_bdryfloat, array

left boundary value. If float, then constant in-flow boundary value applied. If array, then length of array must be equal to ny = grid[0], and this array defines non-constant boundary value on the left.

b_bdryfloat, array

bottom boundary value. If float, then constant in-flow boundary value applied. If array, then length of array must be equal to nx = grid[1], and this array defines non-constant boundary value on the bottom.

Returns:
Acsr_matrix

Defines 2D FD upwind discretization, with boundary

rhsarray

Defines right-hand-side with boundary contributions

See also

poisson

Examples

>>> from numpy import pi
>>> from pyamg.gallery import advection_2d
>>> A, rhs = advection_2d( (4,4), theta=pi/4)
>>> print(A.toarray().round(4))
    [[ 1.4142  0.      0.     -0.7071  0.      0.      0.      0.      0.    ]
     [-0.7071  1.4142  0.      0.     -0.7071  0.      0.      0.      0.    ]
     [ 0.     -0.7071  1.4142  0.      0.     -0.7071  0.      0.      0.    ]
     [ 0.      0.      0.      1.4142  0.      0.     -0.7071  0.      0.    ]
     [ 0.      0.      0.     -0.7071  1.4142  0.      0.     -0.7071  0.    ]
     [ 0.      0.      0.      0.     -0.7071  1.4142  0.      0.     -0.7071]
     [ 0.      0.      0.      0.      0.      0.      1.4142  0.      0.    ]
     [ 0.      0.      0.      0.      0.      0.     -0.7071  1.4142  0.    ]
     [ 0.      0.      0.      0.      0.      0.      0.     -0.7071  1.4142]]
pyamg.gallery.demo()[source]#

Outline basic demo.

pyamg.gallery.diffusion_stencil_2d(epsilon=1.0, theta=0.0, type='FE')[source]#

Rotated Anisotropic diffusion in 2d of the form.

-div Q A Q^T grad u

Q = [cos(theta) -sin(theta)]

[sin(theta) cos(theta)]

A = [1 0 ]

[0 eps ]

Parameters:
epsilonfloat, optional

Anisotropic diffusion coefficient: -div A grad u, where A = [1 0; 0 epsilon]. The default is isotropic, epsilon=1.0

thetafloat, optional

Rotation angle theta in radians defines -div Q A Q^T grad, where Q = [cos(theta) -sin(theta); sin(theta) cos(theta)].

type{‘FE’,’FD’}

Specifies the discretization as Q1 finite element (FE) or 2nd order finite difference (FD) The default is theta = 0.0

Returns:
stencilnumpy array

A 3x3 diffusion stencil

See also

stencil_grid, poisson

Notes

Not all combinations are supported.

Examples

>>> import scipy as sp
>>> from pyamg.gallery.diffusion import diffusion_stencil_2d
>>> sten = diffusion_stencil_2d(epsilon=0.0001,theta=sp.pi/6,type='FD')
>>> print(sten)
[[-0.2164847 -0.750025   0.2164847]
 [-0.250075   2.0002    -0.250075 ]
 [ 0.2164847 -0.750025  -0.2164847]]

Consider a 2 x 4 grid ([x0, x1] x [y0, y1, y2, y3]). The first dimension of the stencil defines x.

>>> nx, ny = (2, 4)
>>> sten = pyamg.gallery.diffusion_stencil_2d(epsilon=0.1, type='FD')
>>> A = pyamg.gallery.stencil_grid(sten, (nx, ny)).toarray()
>>> print(sten)
    [[-0.  -1.   0. ]
     [-0.1  2.2 -0.1]
     [ 0.  -1.  -0. ]]
>>> print(A)
    [[ 2.2 -0.1  0.   0.  -1.   0.   0.   0. ]
     [-0.1  2.2 -0.1  0.   0.  -1.   0.   0. ]
     [ 0.  -0.1  2.2 -0.1  0.   0.  -1.   0. ]
     [ 0.   0.  -0.1  2.2  0.   0.   0.  -1. ]
     [-1.   0.   0.   0.   2.2 -0.1  0.   0. ]
     [ 0.  -1.   0.   0.  -0.1  2.2 -0.1  0. ]
     [ 0.   0.  -1.   0.   0.  -0.1  2.2 -0.1]
     [ 0.   0.   0.  -1.   0.   0.  -0.1  2.2]]
pyamg.gallery.gauge_laplacian(npts, spacing=1.0, beta=0.1)[source]#

Construct a Gauge Laplacian from Quantum Chromodynamics for regular 2D grids.

Note that this function is not written efficiently, but should be fine for N x N grids where N is in the low hundreds.

Parameters:
nptsint

number of pts in x and y directions

spacingfloat

grid spacing between points

betafloat

temperature Note that if beta=0, then we get the typical 5pt Laplacian stencil

Returns:
Acsr matrix

A is Hermitian positive definite for beta > 0.0 A is Symmetric semi-definite for beta = 0.0

References

[1]

MacLachlan, S. and Oosterlee, C., “Algebraic Multigrid Solvers for Complex-Valued Matrices”, Vol. 30, SIAM J. Sci. Comp, 2008

Examples

>>> from pyamg.gallery import gauge_laplacian
>>> A = gauge_laplacian(10)
pyamg.gallery.linear_elasticity(grid, spacing=None, E=100000.0, nu=0.3, format=None)[source]#

Linear elasticity problem with Q1 finite elements on a regular rectangular grid.

Parameters:
gridtuple

length 2 tuple of grid sizes, e.g. (10, 10)

spacingtuple

length 2 tuple of grid spacings, e.g. (1.0, 0.1)

Efloat

Young’s modulus

nufloat

Poisson’s ratio

formatstring

Format of the returned sparse matrix (eg. ‘csr’, ‘bsr’, etc.)

Returns:
Acsr_matrix

FE Q1 stiffness matrix

Barray

rigid body modes

Notes

  • only 2d for now

References

[1]

J. Alberty, C. Carstensen, S. A. Funken, and R. KloseDOI “Matlab implementation of the finite element method in elasticity” Computing, Volume 69, Issue 3 (November 2002) Pages: 239 - 263 http://www.math.hu-berlin.de/~cc/

Examples

>>> from pyamg.gallery import linear_elasticity
>>> A, B = linear_elasticity((4, 4))
pyamg.gallery.linear_elasticity_p1(vertices, elements, E=100000.0, nu=0.3, format=None)[source]#

P1 elements in 2 or 3 dimensions.

Parameters:
verticesarray_like

array of vertices of a triangle or tets

elementsarray_like

array of vertex indices for tri or tet elements

Efloat

Young’s modulus

nufloat

Poisson’s ratio

formatstring

‘csr’, ‘csc’, ‘coo’, ‘bsr’

Returns:
Acsr_matrix

FE Q1 stiffness matrix

Notes

  • works in both 2d and in 3d

References

[1]

J. Alberty, C. Carstensen, S. A. Funken, and R. KloseDOI “Matlab implementation of the finite element method in elasticity” Computing, Volume 69, Issue 3 (November 2002) Pages: 239 - 263 http://www.math.hu-berlin.de/~cc/

Examples

>>> from pyamg.gallery import linear_elasticity_p1
>>> import numpy as np
>>> E = np.array([[0, 1, 2],[1, 3, 2]])
>>> V = np.array([[0.0, 0.0],[1.0, 0.0],[0.0, 1.0],[1.0, 1.0]])
>>> A, B = linear_elasticity_p1(V, E)
pyamg.gallery.load_example(name)[source]#

Load an example problem by name.

Parameters:
namestring (e.g. ‘airfoil’)

Name of the example to load

Notes

Each example is stored in a dictionary with the following keys:
  • ‘A’ : sparse matrix

  • ‘B’ : near-nullspace candidates

  • ‘vertices’ : dense array of nodal coordinates

  • ‘elements’ : dense array of element indices

Current example names are:

airfoil bar helmholtz_2D knot local_disc_galerkin_diffusion recirc_flow unit_cube unit_square

Examples

>>> from pyamg.gallery import load_example
>>> ex = load_example('knot')
pyamg.gallery.poisson(grid, dtype=<class 'float'>, format=None, type='FD')[source]#

Return a sparse matrix for the N-dimensional Poisson problem.

The matrix represents a finite Difference approximation to the Poisson problem on a regular n-dimensional grid with unit grid spacing and Dirichlet boundary conditions.

The last dimension is iterated over first: z, then y, then x. This should be used with np.mgrid() or np.ndenumerate.

Parameters:
gridtuple of integers

grid dimensions e.g. (100,100)

Notes

The matrix is symmetric and positive definite (SPD).

Examples

>>> from pyamg.gallery import poisson
>>> # 4 nodes in one dimension
>>> poisson((4,)).toarray()
array([[ 2., -1.,  0.,  0.],
       [-1.,  2., -1.,  0.],
       [ 0., -1.,  2., -1.],
       [ 0.,  0., -1.,  2.]])
>>> # rectangular two dimensional grid
>>> poisson((2,3)).toarray()
array([[ 4., -1.,  0., -1.,  0.,  0.],
       [-1.,  4., -1.,  0., -1.,  0.],
       [ 0., -1.,  4.,  0.,  0., -1.],
       [-1.,  0.,  0.,  4., -1.,  0.],
       [ 0., -1.,  0., -1.,  4., -1.],
       [ 0.,  0., -1.,  0., -1.,  4.]])
pyamg.gallery.regular_triangle_mesh(nx, ny)[source]#

Construct a regular triangular mesh in the unit square.

Parameters:
nxint

Number of nodes in the x-direction

nyint

Number of nodes in the y-direction

Returns:
Vertarray

nx*ny x 2 vertex list

E2Varray

Nex x 3 element list

Examples

>>> from pyamg.gallery import regular_triangle_mesh
>>> E2V,Vert = regular_triangle_mesh(3, 2)
pyamg.gallery.sprand(m, n, density, format='csr')[source]#

Return a random sparse matrix.

Parameters:
m, nint

shape of the result

densityfloat

target a matrix with nnz(A) = m*n*density, 0<=density<=1

formatstring

sparse matrix format to return, e.g. ‘csr’, ‘coo’, etc.

Returns:
Asparse matrix

m x n sparse matrix

Examples

>>> from pyamg.gallery import sprand
>>> A = sprand(5,5,3/5.0)
pyamg.gallery.stencil_grid(S, grid, dtype=None, format=None)[source]#

Construct a sparse matrix form a local matrix stencil.

Parameters:
Sndarray

matrix stencil stored in N-d array

gridtuple

tuple containing the N grid dimensions

dtype

data type of the result

formatstring

sparse matrix format to return, e.g. “csr”, “coo”, etc.

Returns:
Asparse matrix

Sparse matrix which represents the operator given by applying stencil S at each vertex of a regular grid with given dimensions.

Notes

The grid vertices are enumerated as arange(prod(grid)).reshape(grid). This implies that the last grid dimension cycles fastest, while the first dimension cycles slowest. For example, if grid=(2,3) then the grid vertices are ordered as (0,0), (0,1), (0,2), (1,0), (1,1), (1,2).

This coincides with the ordering used by the NumPy functions ndenumerate() and mgrid().

Examples

>>> from pyamg.gallery import stencil_grid
>>> stencil = [-1,2,-1]  # 1D Poisson stencil
>>> grid = (5,)          # 1D grid with 5 vertices
>>> A = stencil_grid(stencil, grid, dtype=float, format='csr')
>>> A.toarray()
array([[ 2., -1.,  0.,  0.,  0.],
       [-1.,  2., -1.,  0.,  0.],
       [ 0., -1.,  2., -1.,  0.],
       [ 0.,  0., -1.,  2., -1.],
       [ 0.,  0.,  0., -1.,  2.]])
>>> stencil = [[0,-1,0],[-1,4,-1],[0,-1,0]] # 2D Poisson stencil
>>> grid = (3,3)                            # 2D grid with shape 3x3
>>> A = stencil_grid(stencil, grid, dtype=float, format='csr')
>>> A.toarray()
array([[ 4., -1.,  0., -1.,  0.,  0.,  0.,  0.,  0.],
       [-1.,  4., -1.,  0., -1.,  0.,  0.,  0.,  0.],
       [ 0., -1.,  4.,  0.,  0., -1.,  0.,  0.,  0.],
       [-1.,  0.,  0.,  4., -1.,  0., -1.,  0.,  0.],
       [ 0., -1.,  0., -1.,  4., -1.,  0., -1.,  0.],
       [ 0.,  0., -1.,  0., -1.,  4.,  0.,  0., -1.],
       [ 0.,  0.,  0., -1.,  0.,  0.,  4., -1.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  0., -1.,  4., -1.],
       [ 0.,  0.,  0.,  0.,  0., -1.,  0., -1.,  4.]])