pygimli.solver¶

Overview¶

Functions

Classes

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: 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()

(png, pdf)¶

(png, pdf)¶

(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()

()

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()

()

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()

()

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

Gravimetry in 2D - Part I¶

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 creates an empty ax window.

Parameters

mesh : GIMLI::Mesh or list of meshes

2D or 3D GIMLi mesh

**kwargs :

• fitViewbool [True]

  Scale x and y limits to match the view.

• axaxe [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=None, 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: :gimliapi:`GIMLI::Mesh`

Mesh represents spatial discretization of the calculation domain

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

Cell values

b: value | array | callable(cell, userData) [None]

Cell values. None means the termin is unused.

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.

Returns

u : array

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

Other Parameters

**kwargs

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

Node values

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.

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()

()

Examples using pygimli.solver.solveFiniteElements

Heat equation in 2D¶

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.