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

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

()

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.