pygimli.solver

General physics independent solver interface.

Overview

Functions

H1Norm(u[, S, mesh]) Create (H1) norm for the finite element space.
L2Norm(u[, M, mesh]) Create Lebesgue (L2) norm for the finite element space.
applyBoundaryValues(uB, mesh, uBBC) TODO Documentme.
assembleBC_(bc, mesh, S, rhs, a[, time, …]) Shortcut to apply all boundary conditions.
assembleDirichletBC(S, boundaryPairs[, rhs, …]) Apply Dirichlet boundary condition.
assembleForceVector(mesh, f[, userData]) DEPRECATED use assembleLoadVector instead
assembleLoadVector(mesh, f[, userData]) Create right hand side vector based on the given mesh and load or force values.
assembleNeumannBC(rhs, boundaryPairs[, a, …]) Apply Neumann condition to the system matrix S.
assembleRobinBC(S, boundaryPairs[, rhs, …]) Apply Robin boundary condition.
boundaryNormals(mesh) Collect all boundary outer normal vectors.
cellDataToBoundaryData(mesh, v) TODO Documentme.
cellDataToBoundaryDataMatrix(mesh) TODO Documentme.
cellDataToBoundaryGrad(mesh, v, vGrad) TODO Documentme.
cellDataToCellGrad(mesh, v, CtB) TODO Documentme.
cellDataToCellGrad2(mesh, v) TODO Documentme.
cellToFace(boundary, vec[, harmonic]) DEPRECATED.
cellToFaceArithmetic(boundary, AMM) TODO Documentme.
checkCFL(times, mesh, vMax) Check Courant-Friedrichs-Lewy condition.
crankNicolson(times, theta, S, I, f[, u0, …]) S = constant over time f = constant over time
createFVPostProzessMesh(mesh, u, uDirichlet) Create node based post processing of cell centered FV solutions.
createLoadVector(mesh, f[, userData])
createMassMatrix(mesh[, b]) Create the mass matrix.
createMesh(poly[, quality, area, smooth, …]) Create a mesh for a given geometry polygon.
createStiffnessMatrix(mesh[, a]) Create the Stiffness matrix.
deepcopy(x[, memo, _nil]) Deep copy operation on arbitrary Python objects.
diffusionConvectionKernel(mesh[, a, b, uB, …]) Generate system matrix for diffusion and convection in a velocity field.
div(mesh, v) Return the discrete interpolated divergence field.
divergence(mesh[, F, normMap, order]) Divergence for callable function F((x,y,z)).
findDiffusion(mesh, a, b, c[, nc]) TODO Documentme.
findVelocity(mesh, v, b, c[, nc]) Find velocity for boundary b based on vector field v.
generateBoundaryValue(boundary, arg[, time, …]) Generate a value for the given Boundary.
grad(mesh, u[, r]) Return the discrete interpolated gradient \(\mathbf{v}\) for a given scalar field \(\mathbf{u}\).
greenDiffusion1D(x[, t, a, dim]) Greens function for diffusion operator.
identity(dom[, start, end, scale]) Create identity matrix.
linSolve(A, b[, solver, verbose]) Direct solution after \(\textbf{x}\) using core LinSolver.
linsolve(A, b[, verbose]) DEPRECATED wrong name style
parseArgPairToBoundaryArray(pair, mesh) Parse boundary related pair argument to [ GIMLI::Boundary, value|callable ] list.
parseArgToArray(arg, nDof[, mesh, userData]) Parse array related arguments to create a valid value array.
parseArgToBoundaries(args, mesh) Parse boundary related arguments to create a valid boundary value list: [ GIMLI::Boundary, value|callable ]
parseMapToCellArray(attributeMap, mesh[, …]) Parse a value map to cell attributes.
show([mesh, data]) Mesh and model visualization.
showSparseMatrix(A[, full]) Show the content of a sparse matrix.
solve(mesh, **kwargs) Solve partial differential equation.
solveFiniteElements(mesh[, a, b, f, bc, …]) Solve partial differential equation with Finite Elements.
solveFiniteVolume(mesh[, a, b, f, fn, vel, …]) Solve partial differential equation with Finite Volumes.
triDiagToeplitz(dom, a, l, r[, start, end]) Create tri-diagonal Toeplitz matrix.
unique(a) Return list of unique elements ever seen with preserving order.

Classes

RungeKutta(solver[, verbose]) TODO DOCUMENT ME
WorkSpace

Functions

H1Norm

pygimli.solver.H1Norm(u, S=None, mesh=None)[source]

Create (H1) norm for the finite element space.

Parameters:
u : iterable

Node based value to compute the H1 norm for.

S : Matrix

Stiffness matrix.

mesh : GIMLI::Mesh

Mesh with the FE space to generate S if necessary.

Returns:
ret : float

\(H1(u)\) norm.

L2Norm

pygimli.solver.L2Norm(u, M=None, mesh=None)[source]

Create Lebesgue (L2) norm for the finite element space.

Find the L2 Norm for a solution in the finite element space. \(u\) exact solution \({\bf M}\) Mass matrix, i.e., Finite element identity matrix.

\[\begin{split}L2(f(x)) = || f(x) ||_{L^2} & = (\int |f(x)|^2 \d x)^{1/2} \\ & \approx h (\sum |f(x)|^2 )^{1/2} \\ L2(u) = || u ||_{L^2} & = (\int |u|^2 \d x)^{1/2} \\ & \approx (\sum M (u)) ^{1/2} \\ e_{L2_rel} = \frac{L2(u)}{L2(u)} & = \frac{(\sum M(u))^{1/2}}{(\sum M u)^{1/2}}\end{split}\]

The error for any approximated solution \(u_h\) correlates to the L2 norm of ‘L2Norm(u - u_h, M)’. If you like relative values, you can also normalize this error with ‘L2Norm(u - u_h, M) / L2Norm(u, M)*100’.

Parameters:
u : iterable

Node based value to compute the L2 norm for.

M : Matrix

Mass element matrix.

mesh : GIMLI::Mesh

Mesh with the FE space to generate M if necessary.

Returns:
ret : float

\(L2(u)\) norm.

applyBoundaryValues

pygimli.solver.applyBoundaryValues(uB, mesh, uBBC)[source]

TODO Documentme.

assembleBC_

pygimli.solver.assembleBC_(bc, mesh, S, rhs, a, time=None, userData=None)[source]

Shortcut to apply all boundary conditions.

This is a helper function for the solver call. Shortcut to apply all boundary conditions will only forward to appropriate assemble functions.

assembleDirichletBC

pygimli.solver.assembleDirichletBC(S, boundaryPairs, rhs=None, time=0.0, userData=None, nodePairs=None)[source]

Apply Dirichlet boundary condition.

Apply Dirichlet boundary condition to the system matrix S and rhs vector. The right hand side values for h can be given for each boundary element individually by setting proper boundary pair arguments.

\[u(\textbf{r}, t) = h \quad\text{for}\quad\textbf{r}\quad\text{on}\quad\delta\Omega=\Gamma_{\text{Dirichlet}}\]
Parameters:
S : GIMLI::RSparseMatrix

System matrix of the system equation.

boundaryPair : list()

List of pairs [ GIMLI::Boundary, h ]. The value \(h\) will assigned to the nodes of the boundaries. Later assignment overwrites prior.

\(h\) need to be a scalar value (float or int) or a value generator function that will be executed at runtime. See pygimli.solver.solver.parseArgToBoundaries and Modelling with Boundary Conditions for example syntax,

nodePairs : list()

List of pairs [ nodeID, uD ]. The value uD will assigned to the nodes given there ids. This node value settings will overwrite any prior settings due to boundaryPair.

rhs : GIMLI::RVector

Right hand side vector of the system equation will bet set to \(u_{\text{D}}\)

time : float

Will be forwarded to value generator.

userData : class

Will be forwarded to value generator.

assembleForceVector

pygimli.solver.assembleForceVector(mesh, f, userData=None)[source]

DEPRECATED use assembleLoadVector instead

assembleLoadVector

pygimli.solver.assembleLoadVector(mesh, f, userData=None)[source]

Create right hand side vector based on the given mesh and load or force values.

Create right hand side vector based on the given mesh and load or force values.

Parameters:
f: float, array, callable(cell, [rhs], [userData]), […]

Your callable need to return a load value for each cell with optional userData. f_cell = f(cell, [userData=None])

Force Values float -> ones(mesh.nodeCount()) * vals, for each node [0 .. mesh.nodeCount()] for each cell [0 .. mesh.cellCount()]

Returns:
rhs: pg.Vector(mesh.nodeCount())

Right hand side Load vector of

assembleNeumannBC

pygimli.solver.assembleNeumannBC(rhs, boundaryPairs, a=None, time=0.0, userData=None)[source]

Apply Neumann condition to the system matrix S.

Apply Neumann condition to the system matrix S. The right hand side values for g can be given for each boundary element individually by setting proper boundary pair arguments.

\[\frac{\partial u(\textbf{r}, t)}{\partial\textbf{n}} = \textbf{n}\nabla u(\textbf{r}, t) = g \quad\text{for}\quad\textbf{r}\quad\text{on}\quad\delta\Omega=\Gamma_{\text{Neumann}}\]
Parameters:
rhs : Vector

Right hand side vector of length node count.

boundaryPair : list()

List of pairs [ GIMLI::Boundary, g ]. The value \(g\) will assigned to the nodes of the boundaries. Later assignment overwrites prior.

\(g\) need to be a scalar value (float or int) or a value generator function that will be executed at run time.

See pygimli.solver.solver.parseArgToBoundaries and Modelling with Boundary Conditions for example syntax,

a : iterable

Per cell values to scale the Neumann part regarding weak formulation.

time : float

Will be forwarded to value generator.

userData : class

Will be forwarded to value generator.

assembleRobinBC

pygimli.solver.assembleRobinBC(S, boundaryPairs, rhs=None, time=0.0, userData=None)[source]

Apply Robin boundary condition.

Apply Robin boundary condition to the system matrix S and rhs vector (if needed for b != 1 and g != 0).

\[\begin{split}\alpha u(\textbf{r}, t) + \beta\frac{\partial u(\textbf{r}, t)}{\partial\textbf{n}} = \gamma \quad\text{for}\quad\textbf{r}\quad\text{on}\quad\delta\Omega=\Gamma_{\text{Robin}}\\ \quad\text{currently only with}\quad \beta = 1 \quad\text{and}\quad \gamma = 0\end{split}\]
TODO
  • b!=1 and g!=0 variable
  • check for b = 0 and move to dirichlet
Parameters:
S : GIMLI::SparseMatrix

System matrix of the system equation.

boundaryPair : list()

List of pairs [ GIMLI::Boundary, \(\alpha\) ]. The value \(\alpha\) will assigned to the nodes of the boundaries. Later assignment overwrites prior.

\(\alpha\) needs to be a scalar value (float or int) or a value generator function that will be executed at run time. See pygimli.solver.solver.parseArgToBoundaries and Modelling with Boundary Conditions for example syntax,

time : float

Will be forwarded to value generator.

userData : class

Will be forwarded to value generator.

boundaryNormals

pygimli.solver.boundaryNormals(mesh)[source]

Collect all boundary outer normal vectors.

cellDataToBoundaryData

pygimli.solver.cellDataToBoundaryData(mesh, v)[source]

TODO Documentme.

cellDataToBoundaryDataMatrix

pygimli.solver.cellDataToBoundaryDataMatrix(mesh)[source]

TODO Documentme.

cellDataToBoundaryGrad

pygimli.solver.cellDataToBoundaryGrad(mesh, v, vGrad)[source]

TODO Documentme.

cellDataToCellGrad

pygimli.solver.cellDataToCellGrad(mesh, v, CtB)[source]

TODO Documentme.

cellDataToCellGrad2

pygimli.solver.cellDataToCellGrad2(mesh, v)[source]

TODO Documentme.

cellToFace

pygimli.solver.cellToFace(boundary, vec, harmonic=False)[source]

DEPRECATED.

Interpolate cell to face values by weighted arithmetic/harmonic mean.

cellToFaceArithmetic

pygimli.solver.cellToFaceArithmetic(boundary, AMM)[source]

TODO Documentme.

checkCFL

pygimli.solver.checkCFL(times, mesh, vMax)[source]

Check Courant-Friedrichs-Lewy condition.

For advection and flow problems. CFL Number should be lower then 1 to ensure stability.

crankNicolson

pygimli.solver.crankNicolson(times, theta, S, I, f, u0=None, progress=None, debug=None)[source]

S = constant over time f = constant over time

createFVPostProzessMesh

pygimli.solver.createFVPostProzessMesh(mesh, u, uDirichlet)[source]

Create node based post processing of cell centered FV solutions.

Create a mesh suitable for node based post processing of cell centered Finite Volume solutions. This is something like cellDataToNodeData with extra Dirichlet points but without smoothing due to interpolation.

IMPROVE DOC!!

createLoadVector

pygimli.solver.createLoadVector(mesh, f, userData=None)[source]

createMassMatrix

pygimli.solver.createMassMatrix(mesh, b=None)[source]

Create the mass matrix.

Calculates the Mass matrix (Finite element identity matrix) the given mesh.

..math::
Parameters:
mesh : GIMLI::Mesh

Arbitrary mesh to calculate the mass element matrix. Type of base and shape functions depends on the cell types.

b : array

Per cell values. If None given default is 1.

Returns:
A : GIMLI::RSparseMatrix

Mass element matrix

createMesh

pygimli.solver.createMesh(poly, quality=30, area=0.0, smooth=None, switches=None, verbose=False, **kwargs)[source]

Create a mesh for a given geometry polygon.

The mesh is created by triangle or tetgen if the gimli support for these mesh generators are installed. The geometry needs to contain nodes and boundaries and should be valid in the sense that the boundaries are non intersecting.

If poly is a list of coordinates a simple Delaunay mesh of the convex hull will be created. TODO: Tetgen support need to be implemented

Parameters:
poly: :gimliapi:`GIMLI::Mesh` or list or ndarray
  • 2D or 3D gimli mesh that contains the PLC.
  • 2D mesh needs edges
  • 3D mesh needs a plc and tetgen as system component
  • List of x y pairs [[x0, y0], … ,[xN, yN]]
  • ndarray [x_i, y_i]
  • PLC or list of PLC
quality: float

2D triangle quality sets a minimum angle constraint. Be careful with values above 34 degrees.

area: float

2D maximum triangle size in m*m

smooth: tuple

[smoothing algorithm, number of iterations] 0, no smoothing 1, node center 2, weighted node center

switches: str

Force triangle to use the gives command switches.

Returns:
mesh: :gimliapi:`GIMLI::Mesh`

Examples

>>> # no need to import matplotlib. pygimli's show does
>>> import pygimli as pg
>>> import pygimli.meshtools as mt
>>> rect = mt.createRectangle(start=[0, 0], end=[4, 1])
>>> ax, _ = pg.show(mt.createMesh(rect, quality=10))
>>> ax, _ = pg.show(mt.createMesh(rect, quality=33))
>>> ax, _ = pg.show(mt.createMesh(rect, quality=33, area=0.01))
>>> pg.wait()
../../_images/pygimli-solver-1_00.png

(png, pdf)

../../_images/pygimli-solver-1_01.png

(png, pdf)

../../_images/pygimli-solver-1_02.png

(png, pdf)

createStiffnessMatrix

pygimli.solver.createStiffnessMatrix(mesh, a=None)[source]

Create the Stiffness matrix.

Calculates the Stiffness matrix \({\bf S}\) for the given mesh scaled with the per cell values a.

..math::
Parameters:
mesh : GIMLI::Mesh

Arbitrary mesh to calculate the stiffness for. Type of base and shape functions depends on the cell types.

a : array, either complex or real, callable

Per cell values., e.g., physical parameter. If None given default is 1.

Returns:
A : GIMLI::RSparseMatrix

Stiffness matrix

deepcopy

pygimli.solver.deepcopy(x, memo=None, _nil=[])[source]

Deep copy operation on arbitrary Python objects.

See the module’s __doc__ string for more info.

diffusionConvectionKernel

pygimli.solver.diffusionConvectionKernel(mesh, a=None, b=0.0, uB=None, duB=None, vel=0, fn=None, scheme='CDS', sparse=False, time=0.0, userData=None)[source]

Generate system matrix for diffusion and convection in a velocity field.

Particle concentration u inside a velocity field.

Peclet Number - ratio between convection/diffusion = F/D
F = velocity flow trough volume boundary, D = diffusion coefficient
Parameters:
mesh : GIMLI::Mesh

Mesh represents spatial discretization of the calculation domain

a : value | array | callable(cell, userData)

Diffusion coefficient per cell

b : value | array | callable(cell, userData)

TODO What is b

fn : iterable(cell)

TODO What is fn

vel : ndarray (N,dim) | RMatrix(N,dim)

velocity field [[v_i,]_j,] with i=[1..3] for the mesh dimension and j = [0 .. N-1] per Cell or per Node so N is either mesh.cellCount() or mesh.nodeCount()

scheme : str [CDS]

Finite volume scheme

  • CDS – Central Difference Scheme.
    maybe irregular for Peclet no. |F/D| > 2 Diffusion dominant. Error of order 2
  • UDS – Upwind Scheme.
    Convection dominant. Error of order 1
  • HS – Hybrid Scheme.
    Diffusion dominant for Peclet-number |(F/D)| < 2 Convection dominant else.
  • PS – Power Law Scheme.
    Identical to HS for Peclet-number |(F/D)| > 10 and near to ES else Convection dominant.
  • ES – Exponential scheme
    Only stationary one-dimensional but exact solution
Returns:
S : GIMLI::SparseMatrix | numpy.ndarray(nCells, nCells)

Kernel matrix, depends on vel, a, b, scheme, uB, duB .. if some of this has been changed you cannot cache these matrix

rhsBoundaryScales : ndarray(nCells)

RHS offset vector

div

pygimli.solver.div(mesh, v)[source]

Return the discrete interpolated divergence field.

Return the discrete interpolated divergence field. \(\mathbf{u}\) for each cell for a given vector field \(\mathbf{v}\). First order integration via boundary center.
\[\begin{split}d(cells) & = \nabla\cdot\vec{v} \\ d(c_i) & = \sum_{j=0}^{N_B}\vec{v}_{B_j} \cdot \vec{n}_{B_j}\end{split}\]
Parameters:
mesh : GIMLI::Mesh

Discretization base, interpolation will be performed via finite element base shape functions.

V : array(N,3) | R3Vector

Vector field at cell centers or boundary centers

Returns:
d : array(M)

Array of divergence values for each cell in the given mesh.

Examples

>>> import pygimli as pg
>>> import numpy as np
>>> v = lambda p: p
>>> mesh = pg.createGrid(x=np.linspace(0, 1, 4))
>>> print(pg.round(pg.solver.div(mesh, v(mesh.boundaryCenters())), 1e-5))
3 [1.0, 1.0, 1.0]
>>> print(pg.round(pg.solver.div(mesh, v(mesh.cellCenters())), 1e-5))
3 [0.5, 1.0, 0.5]
>>> mesh = pg.createGrid(x=np.linspace(0, 1, 4),
...                      y=np.linspace(0, 1, 4))
>>> print(pg.round(pg.solver.div(mesh, v(mesh.boundaryCenters())), 1e-5))
9 [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0]
>>> divCells = pg.solver.div(mesh, v(mesh.cellCenters()))
>>> # divergence from boundary values are exact where the divergence from
>>> # interpolated cell center values are wrong due to boundary interp
>>> print(sum(divCells))
12.0
>>> mesh = pg.createGrid(x=np.linspace(0, 1, 4),
...                      y=np.linspace(0, 1, 4),
...                      z=np.linspace(0, 1, 4))
>>> print(sum(pg.solver.div(mesh, v(mesh.boundaryCenters()))))
81.0
>>> divCells = pg.solver.div(mesh, v(mesh.cellCenters()))
>>> print(sum(divCells))
54.0

divergence

pygimli.solver.divergence(mesh, F=None, normMap=None, order=1)[source]

Divergence for callable function F((x,y,z)).

MOVE THIS to a better place

Divergence for callable function F((x,y,z)). Return sum div over boundary.

findDiffusion

pygimli.solver.findDiffusion(mesh, a, b, c, nc=None)[source]

TODO Documentme.

Parameters:
a :

Attribute for diffusion coefficient. Cell based or Boundary based

b : GIMLI::Boundary

Boundary

c : GIMLI::Cell

associated Cell in flow direction

nc : GIMLI::Cell

associated neighbor cell .. if one given from flow direction

findVelocity

pygimli.solver.findVelocity(mesh, v, b, c, nc=None)[source]

Find velocity for boundary b based on vector field v.

Parameters:
mesh : GIMLI::Mesh
v : array [(N,dim)]

velocity field [[v_i,]_j,] with i=[1..3] for the mesh dimension and j = [0 .. N-1] per Cell or per Node so N is either mesh.cellCount() or mesh.nodeCount()

b : GIMLI::Boundary

Boundary

c : GIMLI::Cell

Associated Cell in flow direction

nc : GIMLI::Cell

Associated neighbor cell .. if one given from flow direction

generateBoundaryValue

pygimli.solver.generateBoundaryValue(boundary, arg, time=0.0, userData=None)[source]

Generate a value for the given Boundary.

Parameters:
boundary : GIMLI::Boundary or list of ..

The related boundary.

arg : convertible | iterable | callable or list of ..
  • convertible into float
  • iterable of minimum length = boundary.id()
  • callable generator function

If arg is a callable it must fulfill:

:: arg(boundary=:gimliapi:GIMLI::Boundary, time=0.0, userData=None)

The callable function arg have to return appropriate values for all nodes of the boundary or one scalar value for all nodes.

Returns:
val : float or list of ..

Examples

>>> import pygimli as pg
>>>
>>> # def uBoundary(boundary):
>>> #    u = []
>>> #    for n in boundary.nodes():
>>> #        u.append(n[0] + n[1])
>>> #    return u

grad

pygimli.solver.grad(mesh, u, r=None)[source]

Return the discrete interpolated gradient \(\mathbf{v}\) for a given scalar field \(\mathbf{u}\).

\[\begin{split}\mathbf{v}(\mathbf{r}_{\mathcal{C}}) &= \nabla u(\mathbf{r}_{\mathcal{N}}) \\ (\mathbf{v_x}(\mathbf{r}_{\mathcal{C}}), \mathbf{v_y}(\mathbf{r}_{\mathcal{C}}), \mathbf{v_z}(\mathbf{r}_{\mathcal{C}}))^{\text{T}} &= \left(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}, \frac{\partial u}{\partial z}\right)^{\text{T}}\end{split}\]

With \(\mathcal{N}=\cup_{i=0}^{N} \text{Node}_i\), \(\mathcal{C}=\cup_{j=0}^{M} \text{Cell}_j\), \(\mathbf{u}=\{u(\mathbf{r}_{i})\}\in I\!R\) and \(\mathbf{r}_{i} = (x_i, y_i, z_i)^{\text{T}}\)

The discrete scalar field \(\mathbf{u}(\mathbf{r}_{\mathcal{N}})\) need to be defined for each node position \(\mathbf{r}_{\mathcal{N}}\). The resulting vector field \(\mathbf{v}(\mathbf{r}_{\mathcal{C}})\) is defined for each cell center position \(\mathbf{r}_{\mathcal{C}}\). If you need other positions than the cell center, provide an appropriate array of coordinates \(\mathbf{r}\).

Parameters:
mesh : GIMLI::Mesh

Discretization base, interpolation will be performed via finite element base shape functions.

u : array | callable

Scalar field per mesh node position or an appropriate callable([[x,y,z]])

r : ndarray((M, 3)) [mesh.cellCenter()]

Alternative target coordinates :math:`mathbf{r} for the resulting gradient field. i.e., the positions where the vector field is defined. Default are all cell centers.

Returns:
v : ndarray((M, 3))

Resulting vector field defined on \(\mathbf{v}(\mathbf{r}_{\mathcal{C}})\). M is number of cells or length of given alternative coordinates r.

See also

GIMLI
:Mesh::cellDataToBoundaryGradient
GIMLI
:Mesh::boundaryDataToCellGradient

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import pygimli as pg
>>> fig, ax = plt.subplots()
>>> mesh = pg.createGrid(x=np.linspace(0, 1, 20), y=np.linspace(0, 1, 20))
>>> u = lambda p: pg.x(p)**2 * pg.y(p)
>>> _ = pg.show(mesh, u(mesh.nodeCenters()), ax=ax)
>>> _ = pg.show(mesh, [2.*pg.y(mesh.cellCenters())*pg.x(mesh.cellCenters()),
...             pg.x(mesh.cellCenters())**2], ax=ax)
>>> _ = pg.show(mesh, pg.solver.grad(mesh, u), ax=ax, color='w',
...             linewidth=0.4)
>>> plt.show()

(png, pdf)

../../_images/pygimli-solver-2.png

greenDiffusion1D

pygimli.solver.greenDiffusion1D(x, t=0, a=1, dim=1)[source]

Greens function for diffusion operator.

Provides the elementary solution for:

\[g(x,t) = \partial t + a \Delta\]

To find a solution for:

\[u(x,t) \quad\text{for}\quad\frac{\partial u(x,t)}{\partial t} + a \Delta u(x,t) = f(x)\]
\[x = [-x, 0, x]\]
\[u(x,t) = g(x,t) * f(x)\]
Parameters:
x : array_like

Discrete spatial coordinates. Should better be symmetric [-x, 0, x].

t : float

Time

a : float, optional

Scale for Laplacian operator

dim : int, optional

Spatial dimension [1]

Returns:
g : array_like

Discrete Greens’function

Examples

>>> import numpy as np
>>> from pygimli.solver import greenDiffusion1D
>>> dx = 0.001
>>> x = np.arange(0, 0.1+dx, dx)
>>> g = greenDiffusion1D(np.hstack((-x[:0:-1], x)), t=1.0, a=1e-4)
>>> g *= dx
>>> f = np.zeros(len(x))
>>> f[int(len(x)/2)] = 1.
>>> u = np.convolve(g, f)[len(x)-1:2*len(x)-1]
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> _ = ax.plot(x, u, label="u(x)")
>>> _ = ax.plot(x, g[::2], label="g(x)")
>>> _ = ax.legend()
>>> fig.show()

(png, pdf)

../../_images/pygimli-solver-3.png

identity

pygimli.solver.identity(dom, start=0, end=-1, scale=1)[source]

Create identity matrix.

linSolve

pygimli.solver.linSolve(A, b, solver=None, verbose=False)[source]

Direct solution after \(\textbf{x}\) using core LinSolver.

\[\textbf{A}\textbf{x} = \textbf{b}\]

If \(\textbf{A}\) is symmetric, sparse and positive definite.

Parameters:
A : GIMLI::RSparseMatrix | :gimliapi:`GIMLI::RSparseMapMatrix`|

numpy.array System matrix. Need to be symmetric, sparse and positive definite.

b : iterable array

Right hand side of the equation.

solver : str [None]

Try to choose a solver, ‘pg’ for pygimli core cholmod or umpfack. ‘np’ for numpy linalg or scipy.sparse.linalg. Automatic choosing if solver is None depending on matrixtype.

verbose : bool [False]

Be verbose.

Returns:
x : GIMLI::RVector

Solution vector

linsolve

pygimli.solver.linsolve(A, b, verbose=False)[source]

DEPRECATED wrong name style

parseArgPairToBoundaryArray

pygimli.solver.parseArgPairToBoundaryArray(pair, mesh)[source]

Parse boundary related pair argument to [ GIMLI::Boundary, value|callable ] list.

Parameters:
pair : tuple
  • [marker, arg]
  • [[marker, …], arg]
  • [boundary, arg]
  • [[boundary,…], arg]
  • [node, arg]
  • [marker, callable, *kwargs]
  • [[marker, …], callable, *kwargs]

arg will be parsed by pygimli.solver.solver.generateBoundaryValue and distributed to each boundary. Callable functions will be executed at run time.

mesh : GIMLI::Mesh

Used to find boundaries by marker.

Returns:
boundaries : list()

[ GIMLI::Boundary, value|callable ]

parseArgToArray

pygimli.solver.parseArgToArray(arg, nDof, mesh=None, userData=None)[source]

Parse array related arguments to create a valid value array.

Parameters:
arg : float | int | iterable | callable

The target array value that will be converted to an array.

If arg is a callable with it must fulfill:

:: arg(cell|node|boundary, userData=None)

Where MeshEntity is one of GIMLI::Cell , GIMLI::Node or GIMLI::Boundary depeding on nDof, where nDof is mesh.cellCount(), mesh.nodeCount() or mesh.boundaryCount(), respectively.

nDof : int | [int]

Desired array size.

mesh : GIMLI::Mesh

Used if arg is callable

userData : class

Used if arg is callable

Returns:
ret : GIMLI::RVector

Array of desired length filled with the appropriate values.

parseArgToBoundaries

pygimli.solver.parseArgToBoundaries(args, mesh)[source]

Parse boundary related arguments to create a valid boundary value list: [ GIMLI::Boundary, value|callable ]

Parameters:
args : callable | pair | [pair, …]

If args is just a callable than every boundary will be evaluated at runtime with this function as args(boundary). Else see pygimli.solver.solver.parseArgPairToBoundaryArray

mesh : GIMLI::Mesh

Used to find boundaries by marker

Returns:
boundaries : list()

[ GIMLI::Boundary, value|callable ]

Examples

>>> # no need to import matplotlib. pygimli's show does
>>> import pygimli as pg
>>> mesh = pg.meshtools.createWorld([0, 0], [1, -1], worldMarker=0)
>>> ax, _ = pg.show(mesh, boundaryMarker=True)
>>> # all edges with marker 1 get value 1.0
>>> b = pg.solver.parseArgToBoundaries([1, 1.0], mesh)
>>> print(len(b))
1
>>> # same as above with marker 2 get value 2
>>> b = pg.solver.parseArgToBoundaries([[1, 1.0], [2, 2.0]], mesh)
>>> print(len(b))
2
>>> # same as above with marker 3 get value 3
>>> b = pg.solver.parseArgToBoundaries([[1, 1.], [2, 2.], [3, 3.]], mesh)
>>> print(len(b))
3
>>> # edges with marker 1 and 2 get value 1
>>> b = pg.solver.parseArgToBoundaries([[1, 2], 1.0], mesh)
>>> print(len(b))
2
>>> b = pg.solver.parseArgToBoundaries([[1, 2, 3], 1.0], mesh)
>>> print(len(b))
3
>>> b = pg.solver.parseArgToBoundaries([[[1, 2, 3], 1.0], [4, 4.0]], mesh)
>>> print(len(b))
4
>>> b = pg.solver.parseArgToBoundaries([mesh.node(0), 0.0], mesh)
>>> print(len(b))
1
>>> def bCall(boundary):
...     u = []
...     for i, n in enumerate(boundary.nodes()):
...         u.append(i)
...     return u
>>> b = pg.solver.parseArgToBoundaries([1, bCall], mesh)
>>> print(len(b),b[0][1](b[0][0]))
1 [0, 1]
>>> def bCall(boundary, a1, a2):
...     return a1 + a2
>>> b = pg.solver.parseArgToBoundaries([1, bCall, {'a1':2, 'a2':3}], mesh)
>>> print(len(b), b[0][1][0](b[0][0], **b[0][1][1]))
1 5
>>> b = pg.solver.parseArgToBoundaries([[1, bCall, {'a1':1, 'a2':2}],
...                                     [2, bCall, {'a1':3, 'a2':4}]], mesh)
>>> print(len(b), b[0][1][0](b[0][0], **b[0][1][1]))
2 3
>>> b = pg.solver.parseArgToBoundaries([[1,2], bCall, {'a1':4, 'a2':5}], mesh)
>>> print(len(b), b[1][1][0](b[1][0], **b[1][1][1]))
2 9
>>> pg.wait()

(png, pdf)

../../_images/pygimli-solver-4.png

parseMapToCellArray

pygimli.solver.parseMapToCellArray(attributeMap, mesh, default=0.0)[source]

Parse a value map to cell attributes.

A map should consist of pairs of marker and value. A marker is an integer and corresponds to the cell.marker().

Parameters:
mesh : GIMLI::Mesh

For each cell of mesh a value will be returned.

attributeMap : list | dict

List of pairs [marker, value] ] || [[marker, value]], or dictionary with marker keys

default : float [0.0]

Fill all unmapped atts to this default.

Returns:
atts : array

Array of length mesh.cellCount()

Examples using pygimli.solver.parseMapToCellArray

show

pygimli.solver.show(mesh=None, data=None, **kwargs)[source]

Mesh and model visualization.

Syntactic sugar to show a mesh with data. Forwards to pygimli.viewer.showMesh or pygimli.viewer.mayaview.showMesh3D to show most of the typical 2D and 3D content. See tutorials and examples for usage hints. An empty show call create an empty ax window.

Parameters:
mesh : GIMLI::Mesh or list of meshes

2D or 3D GIMLi mesh

**kwargs :
  • fitView : bool [True]
    Scale x and y limits to match the view.
  • ax : axe [None]
    Matplotlib axes object. Create a new if necessary.
  • Will be forwarded to the appropriate show functions.
Returns:
Return the results from the showMesh* functions.

See also

showMesh

showSparseMatrix

pygimli.solver.showSparseMatrix(A, full=False)[source]

Show the content of a sparse matrix.

Parameters:
A : GIMLI::SparseMatrix | GIMLI::SparseMapMatrix

Matrix to be shown.

full : bool [False]

Show as dense matrix.

solve

pygimli.solver.solve(mesh, **kwargs)[source]

Solve partial differential equation.

This is a syntactic sugar convenience function for solving partial differential equation on a given mesh. Using the best guess method for the given parameter. Currently only Finite Element calculation using pygimli.solver.solveFiniteElements

TODO pygimli.solver.solveFiniteVolume

solveFiniteElements

pygimli.solver.solveFiniteElements(mesh, a=1.0, b=0.0, f=0.0, bc=None, times=None, userData=None, verbose=False, **kwargs)[source]

Solve partial differential equation with Finite Elements.

This is a syntactic sugar convenience function for using the Finite Element functionality of the library core to solve partial differential equation (PDE) that match the following form:

\[\begin{split}\frac{\partial u}{\partial t} & = \nabla\cdot(a \nabla u) + b u + f(\mathbf{r},t)~~|~~\Omega_{\text{Mesh}}\\ u & = h~~|~~\Gamma_{\text{Dirichlet}}\\ \frac{\partial u}{\partial \mathbf{n}} & = g~~|~~\Gamma_{\text{Neumann}}\\ \alpha u + \beta\frac{\partial u}{\partial \mathbf{n}} & = \gamma~~|~~\Gamma_{\text{Robin}}\end{split}\]

for the scalar solution \(u(\mathbf{r}, t)\) at each node of a given mesh. The Domain \(\Omega\) and the Boundary \(\Gamma\) are defined through the mesh with appropriate boundary marker.

TODO:

  • unsteady ub and dub
  • ‘Infinity’ Boundary condition (u vanishes at infinity)
  • ‘Cauchy’ Boundary condition

(guaranties u and du on same boundary, will never work here because the problem becomes ill posed and would need some inverse strategy to solve.)

Parameters:
mesh : GIMLI::Mesh

Mesh represents spatial discretization of the calculation domain

a : value | array | callable(cell, userData)

Cell values

b : value | array | callable(cell, userData)

Cell values

f : value | array(cells) | array(nodes) | callable(args, kwargs)

force values

bc : dict()

Dictionary of boundary conditions. Current supported boundary conditions are by dictionary keys: ‘Dirichlet’, ‘Neumann’, ‘Robin’.

The dictionary can contain multiple “key: Arg” Arg will be parsed by pygimli.solver.solver.parseArgPairToBoundaryArray

If the dictionary key is ‘Node’ then fixed values for single node indices can by be given. e.g., bc={‘Node’: [nodeID, value]}. Note this is only a shortcut for bc={‘Dirichlet’: [mesh.node(nodeID), value]}.

times : array [None]

Solve as time dependent problem for the given times.

**kwargs:
u0 : value | array | callable(pos, userData)

Node values

uB : value | array | callable(pos, userData)

DEPRECATED use bc={‘Dirichlet’ | uB}

uN : list([node, value])

DEPRECATED use bc={‘Node’ | uN}

duB : value | array | callable(pos, userData)

DEPRECATED use bc={‘Neumann’ | duB}

theta : float [1]
  • \(theta = 0\) means explicit Euler, maybe stable for

\(\Delta t \quad\text{near}\quad h\) - \(theta = 0.5\), Crank-Nicolsen, maybe instable - \(theta = 1\), implicit Euler

If unsure choose \(\theta = 0.5 + \epsilon\), which is probably stable.

stats : bool

Give some statistics.

progress : bool

Give some calculation progress.

assembleOnly : bool

Stops after matrix asssemblation. Returns the system matrix A and the rhs vector.

ws : dict

The WorkSpace is a dictionary that will get some temporary data during the calculation. Any keyvalue ‘u’ in the dictionary will be used for the resulting array.

Returns:
u : array

Returns the solution u either 1,n array for stationary problems or for m,n array for m time steps

Examples

>>> import pygimli as pg
>>> from pygimli.meshtools import polytools as plc
>>> from pygimli.mplviewer import drawField, drawMesh
>>> import matplotlib.pyplot as plt
>>> world = plc.createWorld(start=[-10, 0], end=[10, -10],
...                         marker=1, worldMarker=False)
>>> c1 = plc.createCircle(pos=[0.0, -5.0], radius=3.0, area=.1, marker=2)
>>> mesh = pg.meshtools.createMesh([world, c1], quality=34.3)
>>> u = pg.solver.solveFiniteElements(mesh, a=[[1, 100], [2, 1]],
...                                   bc={'Dirichlet':[[4, 1.0], [2, 0.0]]})
>>> fig, ax = plt.subplots()
>>> pc = drawField(ax, mesh, u)
>>> drawMesh(ax, mesh)
>>> plt.show()

(png, pdf)

../../_images/pygimli-solver-5.png

Examples using pygimli.solver.solveFiniteElements

solveFiniteVolume

pygimli.solver.solveFiniteVolume(mesh, a=1.0, b=0.0, f=0.0, fn=0.0, vel=None, u0=0.0, times=None, uL=None, relax=1.0, ws=None, scheme='CDS', **kwargs)[source]

Solve partial differential equation with Finite Volumes.

This function is a syntactic sugar proxy for using the Finite Volume functionality of the library core to solve elliptic and parabolic partial differential of the following type:

\[\begin{split}\frac{\partial u}{\partial t} + \mathbf{v}\cdot\nabla u & = \nabla\cdot(a \nabla u) + b u + f(\mathbf{r},t)\\ u(\mathbf{r}, t) & = u_B \quad\mathbf{r}\in\Gamma_{\text{Dirichlet}}\\ \frac{\partial u(\mathbf{r}, t)}{\partial \mathbf{n}} & = u_{\partial \text{B}} \quad\mathbf{r}\in\Gamma_{\text{Neumann}}\\ u(\mathbf{r}, t=0) & = u_0 \quad\text{with} \quad\mathbf{r}\in\Omega\end{split}\]

The Domain \(\Omega\) and the Boundary \(\Gamma\) are defined through the given mesh with appropriate boundary marker.

The solution \(u(\mathbf{r}, t)\) is given for each cell in the mesh.

TODO:

  • Refactor with solver class and Runga-Kutte solver
  • non steady boundary conditions
Parameters:
mesh : GIMLI::Mesh

Mesh represents spatial discretization of the calculation domain

a : value | array | callable(cell, userData)

Stiffness weighting per cell values.

b : value | array | callable(cell, userData)

Scale for mass values b

f : iterable(cell)

Load vector

fn : iterable(cell)

TODO What is fn

vel : ndarray (N,dim) | RMatrix(N,dim)

Velocity field \(\mathbf{v}(\mathbf{r}, t=\text{const}) = \{[v_i]_j,\}\) with \(i=[1\ldots 3]\) for the mesh dimension and \(j = [0\ldots N-1]\) with N either the amount of cells, nodes, or boundaries. Velocities per boundary are preferred and will be interpolated on demand.

u0 : value | array | callable(cell, userData)

Starting field

times : iterable

Time steps to calculate for.

ws : Workspace

This can be an empty class that will used as an Workspace to store and cache data.

If ws is given: The system matrix is taken from ws or calculated once and stored in ws for further usage.

The system matrix is cached in this Workspace as ws.S The LinearSolver with the factorized matrix is cached in this Workspace as ws.solver The rhs vector is only stored in this Workspace as ws.rhs

scheme : str [CDS]

Finite volume scheme: pygimli.solver.diffusionConvectionKernel

**kwargs:
  • uB : Dirichlet boundary conditions
  • duB : Neumann boundary conditions
Returns:
u : ndarray(nTimes, nCells)

solution field for all time steps

triDiagToeplitz

pygimli.solver.triDiagToeplitz(dom, a, l, r, start=0, end=-1)[source]

Create tri-diagonal Toeplitz matrix.

unique

pygimli.solver.unique(a)[source]

Return list of unique elements ever seen with preserving order.

See also

unique_everseen, unique_rows

Examples

>>> from pygimli.utils import unique
>>> unique((1,1,2,2,3,1))
[1, 2, 3]

Classes

RungeKutta

class pygimli.solver.RungeKutta(solver, verbose=False)[source]

TODO DOCUMENT ME

Methods

run(u0, dt[, tMax]) TODO DOCUMENT_ME
start(u0, dt[, tMax]) TODO DOCUMENT_ME
step() TODO DOCUMENT ME
__init__(solver, verbose=False)[source]

TODO DOCUMENT_ME

rk4a = [0.0, -0.41789047449985195, -1.192151694642677, -1.6977846924715279, -1.5141834442571558]
rk4b = [0.14965902199922912, 0.37921031299962726, 0.8229550293869817, 0.6994504559491221, 0.15305724796815198]
rk4c = [0.0, 0.14965902199922912, 0.37040095736420475, 0.6222557631344432, 0.9582821306746903]
run(u0, dt, tMax=1)[source]

TODO DOCUMENT_ME

start(u0, dt, tMax=1)[source]

TODO DOCUMENT_ME

step()[source]

TODO DOCUMENT ME

WorkSpace

class pygimli.solver.WorkSpace[source]
__init__($self, /, *args, **kwargs)

Initialize self. See help(type(self)) for accurate signature.



2018 - GIMLi Development Team
Created using Bootstrap, Sphinx and pyGIMLi 1.0.9+87.g68804698 on Dec 16, 2018.