pygimli.solver¶
General physics independent solver interface.
Overview¶
Functions
|
Create anisotropy matrix with desired properties. |
|
TODO Documentme. |
|
Shortcut to apply all boundary conditions. |
|
Apply Dirichlet boundary condition. |
|
Assemble the load vector. |
|
Apply Neumann condition to the system matrix S. |
|
Apply Robin boundary condition. |
|
Collect all boundary outer normal vectors. |
|
TODO Documentme. |
TODO Documentme. |
|
|
TODO Documentme. |
|
TODO Documentme. |
|
TODO Documentme. |
|
DEPRECATED. |
|
TODO Documentme. |
|
Get a value for each cell. |
|
Check Courant-Friedrichs-Lewy condition. |
|
Create constitutive matrix for 2 or 3D isotropic media. |
|
S = constant over time f = constant over time |
|
Create node based post processing of cell centered FV solutions. |
|
Create a right hand side vector for vector valued solutions. |
|
Create right hand side vector based on the given mesh and load values (scalar solution) or force vectors (vector value solution). |
|
Create the mass matrix. |
|
Create a mesh for a given PLC or point list. |
|
Create the Stiffness matrix. |
|
Deep copy operation on arbitrary Python objects. |
|
Generate system matrix for diffusion and convection in a velocity field. |
|
Return the discrete interpolated divergence field. |
|
Divergence for callable function func((x,y,z)). |
|
TODO Documentme. |
|
Find velocity for boundary b based on vector field v. |
|
Generate a value for the given Boundary. |
|
Return the discrete interpolated gradient \(\mathbf{v}\) for a given scalar field \(\mathbf{u}\). |
|
Greens function for diffusion operator. |
|
Create identity matrix. |
|
Return integral over nodal solution \(u\). |
|
Direct solution after \(\textbf{x}\) using core LinSolver. |
|
Create (H1) norm for the finite element space. |
|
Create Lebesgue (L2) norm for finite element space. |
|
Parse boundary related pair argument to create a list of [ GIMLI::Boundary, value|callable ]. |
|
Parse array related arguments to create a valid value array. |
|
Parse boundary related arguments to create a valid boundary value list: [ GIMLI::Boundary, value|callable ] |
|
Parse dictionary key of type str to marker list. |
|
Parse a value map to cell attributes. |
|
Mesh and model visualization. |
|
Show the content of a sparse matrix. |
|
Solve partial differential equation. |
|
Solve partial differential equation with Finite Elements. |
|
Solve partial differential equation with Finite Volumes. |
|
Create tri-diagonal Toeplitz matrix. |
Classes
|
TODO DOCUMENT ME |
Functions¶
anisotropyMatrix¶
-
pygimli.solver.anisotropyMatrix(lon, trans, theta)[source]¶ Create anisotropy matrix with desired properties.
Anistropy tensor from longitudinal value lon, transverse value trans and the angle theta of the symmetry axis relative to the vertical after cite:WieseGreZho2015 https://www.researchgate.net/publication/249866312_Explicit_expressions_for_the_Frechet_derivatives_in_3D_anisotropic_resistivity_inversion
Examples using pygimli.solver.anisotropyMatrix
assembleBC¶
assembleDirichletBC¶
-
pygimli.solver.assembleDirichletBC(mat, boundaryPairs, rhs=None, time=0.0, userData={}, nodePairs=None, dofOffset=0)[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
mat (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.parseArgToBoundariesand Modelling with Boundary Conditions for example syntax,nodePairs (list() | callable) – 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.
dofOffset (int[0]) – Offset for matrix index.
Examples using pygimli.solver.assembleDirichletBC
assembleLoadVector¶
assembleNeumannBC¶
-
pygimli.solver.assembleNeumannBC(rhs, boundaryPairs, nDim=1, time=0.0, userData={}, dofOffset=0)[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.parseArgToBoundariesand Modelling with Boundary Conditions for example syntax,nDim (int [1]) – Number of dimensions for vector valued problems. The rhs array need top have the correct size, i.e., number of Nodes * mesh.dimension()
time (float) – Will be forwarded to value generator.
userData (class) – Will be forwarded to value generator.
dofOffset (int[0]) – Offset for matrix index.
assembleRobinBC¶
-
pygimli.solver.assembleRobinBC(mat, boundaryPairs, rhs=None, time=0.0, userData={}, dofOffset=0)[source]¶ Apply Robin boundary condition.
Apply Robin boundary condition to the system matrix and the rhs vector
\[\begin{split}\frac{\partial u(\textbf{r}, t)}{\partial\textbf{n}} & = \alpha(u_0-u) \quad\text{or} \\ \beta\frac{\partial u(\textbf{r}, t)}{\partial\textbf{n}} + \alpha u & = \gamma \\ & \quad\text{for}\quad\textbf{r}\quad\text{on}\quad\delta\Omega=\Gamma_{\text{Robin}}\\\end{split}\]- Parameters
mat (GIMLI::SparseMatrix) – System matrix of the system equation.
boundaryPair (list()) –
List of pairs [GIMLI::Boundary, \(a, u_0\) | \(\alpha, \beta, \gamma\)]. The values will assigned to the nodes of the boundaries. Later assignment overwrites prior.
Values can be a single value for \(\alpha\) or \(a\), two values will be interpreted as \(a, u_0\), and three values will be \(\alpha, \beta, \gamma\). Also generator (callable) is possible which will be executed at run time. See
pygimli.solver.solver.parseArgToBoundariesModelling with Boundary Conditions or testing/test_FEM.py for example syntax.time (float) – Will be forwarded to value generator.
userData (dict) – Will be forwarded to value generator.
dofOffset (int[0]) – Offset for matrix index.
boundaryNormals¶
cellDataToBoundaryDataMatrix¶
cellDataToBoundaryGrad¶
cellToFace¶
cellValues¶
-
pygimli.solver.cellValues(mesh, arg, **kwargs)[source]¶ Get a value for each cell.
Returns a array or vector of length mesh.cellCount() based on arg. The preferable arg is a dictionary for the cell marker and the appropriate cell value. The designated value can be calculated using a callable(cell, **kwargs), which is called on demand
- Parameters
mesh (GIMLI::Mesh) – Used if arg is callable
arg (float | int | complex | ndarray | iterable | callable | dict) –
Argument to be parsed as cell data. If arg is a dictionary, the dict key will be interpreted as cell marker:
Dictionary is key: value. Value can be float, int, complex or ndarray. The last for anistropic or elastic tensors.
Key can be integer for cell marker or str, which will be interpreted as splice or list. See examples or parseDictKey.
Iterable of length mesh.nodeCount() will be interpolated to cellCenters.
userData (class) – Used if arg is callable
- Returns
ret – Array of desired length filled with the appropriate values.
- Return type
GIMLI::RVector | ndarray(mesh.cellCount(), xx )
Examples
>>> import pygimli as pg >>> >>> mesh = pg.createGrid(x=range(5)) >>> mesh.setCellMarkers([1, 1, 2, 2]) >>> print(mesh.cellCount()) 4 >>> print(pg.solver.cellValues(mesh, [1, 2, 3, 4])) [1, 2, 3, 4] >>> print(pg.solver.cellValues(mesh, {1:1.0, 2:10})) [1.0, 1.0, 10, 10] >>> print(pg.solver.cellValues(mesh, {':':2.0})) [2.0, 2.0, 2.0, 2.0] >>> print(pg.solver.cellValues(mesh, {'0:2':3.0})) [3.0, 3.0, None, None] >>> print(pg.solver.cellValues(mesh, np.ones(mesh.nodeCount()))) 4 [1.0, 1.0, 1.0, 1.0] >>> print(np.array(pg.solver.cellValues(mesh, {'1:3' : np.diag([1.0, 2.0])}))) [[[1. 0.] [0. 2.]] [[1. 0.] [0. 2.]] [[1. 0.] [0. 2.]] [[1. 0.] [0. 2.]]] >>> print(np.array(pg.solver.cellValues(mesh, {':' : pg.core.CMatrix(2, 2)}))) [[[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]] [[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]] [[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]] [[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]]] >>> print(pg.solver.cellValues(mesh, {'1,2':1 + 1j*2.0})) [(1+2j), (1+2j), (1+2j), (1+2j)] >>> def cellVal(c, b=1): ... return c.center()[0]*b >>> t = pg.solver.cellValues(mesh, {':' : cellVal}) >>> print([t[c.id()](c) for c in mesh.cells()]) [0.5, 1.5, 2.5, 3.5]
checkCFL¶
constitutiveMatrix¶
crankNicolson¶
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!!
createForceVector¶
-
pygimli.solver.createForceVector(mesh, f, userData={})[source]¶ Create a right hand side vector for vector valued solutions.
- Parameters
f ([ convertable ]) – List of rhs side options. Must be convertable to createLoadVector. See
createLoadVectorrhs (np.array()) – Squeezed vector of length mesh.nodeCount() * mesh.dimensions()
createLoadVector¶
-
pygimli.solver.createLoadVector(mesh, f=1.0, userData={})[source]¶ Create right hand side vector based on the given mesh and load values (scalar solution) or force vectors (vector value solution).
Create right hand side based on the given mesh and load or force values.
- Parameters
f (float[1.0], array, callable(cell, [userData]), [f_x, f_y, f_z]) –
float will be assumed as constant for all cells
like rhs = rhs(np.ones(mesh.cellCount() * f), - array of length mesh.cellCount() will be processed as load value for each cell: rhs = rhs(f), - array of length mesh.nodeCount() will be assumed to be allready processed correct: rhs = f - callable is evaluated on once for each cell and need to return a load value for each cell and can have optional a userData dictionary: f_cell = f(cell, [userData={}]) rhs = rhs(f(c, userData) for c in mesh.cells()) - list with length of mesh.dimension() of float or array entries will create a squeezed rhs for vector valued problems rhs = squeeze([rhs(f[0]), rhs(f[1]), rhs(f[2])])
- Returns
rhs – Right hand side load vector for scalar values or squeezed vector values.
- Return type
pg.Vector(mesh.nodeCount())
Examples using pygimli.solver.createLoadVector
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 – Mass element matrix
- Return type
GIMLI::RSparseMatrix
createMesh¶
-
pygimli.solver.createMesh(poly, quality=32, area=0.0, smooth=None, switches=None, verbose=False, **kwargs)[source]¶ Create a mesh for a given PLC or point list.
The mesh is created by triangle or tetgen if the pgGIMLi support for these mesh generators is installed. A PLC 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 with a convex hull will be created. Quality and area arguments are ignored for this case to create a mesh with one node for each coordinate position.
- Parameters
poly (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 PLCs
quality (float) – 2D triangle quality sets a minimum angle constraint. Be careful with values above 34 degrees. 3D tetgen quality. Be careful with values below 1.12.
area (float) – Maximum element size (global). 2D maximum triangle size in m*², 3D maximum tetrahedral size in m³.
smooth (tuple) – [smoothing algorithm, number of iterations] 0: no smoothing 1: node center 2: weighted node center
switches (str) – Set additional triangle command switches. https://www.cs.cmu.edu/~quake/triangle.switch.html
- Returns
mesh
- Return type
>>> 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()
createStiffnessMatrix¶
-
pygimli.solver.createStiffnessMatrix(mesh, a=None, isVector=False)[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 (iterable of type float, int, complex, RMatrix, CMatrix) – Per cell values., e.g., physical parameter. Length of a need to be mesh.cellCount(). If None given default is 1.
isVector (bool [False]) – We want to solve for vector valued problems. Resulting SparseMatrix is a SparseMapMatrix and have the dimension (nNodes * nDims, nNodes * nDims) with nNodes = mesh.nodeCount() and nDims = mesh.dimension().
- Returns
A – Stiffness matrix, with real or complex values.
- Return type
GIMLI::[C]SparseMatrix | [C]SparseMapMatrix
Examples using pygimli.solver.createStiffnessMatrix
deepcopy¶
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 of divergence values for each cell in the given mesh.
- Return type
array(M)
Examples
>>> import pygimli as pg >>> import numpy as np >>> v = lambda p: p >>> mesh = pg.createGrid(x=np.linspace(0, 1, 4)) >>> print(pg.math.round(pg.solver.div(mesh, v(mesh.boundaryCenters())), 1e-5)) 3 [1.0, 1.0, 1.0] >>> print(pg.math.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.math.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 interpolation to the boundary >>> 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¶
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={}, expectList=False)[source]¶ Generate a value for the given Boundary.
- Parameters
boundary (GIMLI::Boundary or list of ..) – The related boundary.
expectList (bool[False]) – Allow list values for Robin BC.
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={})
The callable function arg have to return appropriate values for all nodes of the boundary or one value for all nodes (scalar field only). Value can be scalar or vector field value, e.g., return force values for all nodes at a boundary need to return a ndarray((nodes, dims)) like ‘lambda _b: np.array([[force_x, force_y, force_z] for n in _b.nodes()]).T’
- Returns
val – Value for all nodes of the boundary.
- Return type
[float]
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 – Resulting vector field defined on \(\mathbf{v}(\mathbf{r}_{\mathcal{C}})\). M is number of cells or length of given alternative coordinates r.
- Return type
ndarray((M, 3))
>>> 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.positions()), 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()
Examples using pygimli.solver.grad
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
- Returns
g – Discrete Greens’function
- Return type
array_like
>>> 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()
intDomain¶
linSolve¶
-
pygimli.solver.linSolve(mat, 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
mat (GIMLI::RSparseMatrix | 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 umfpack. ‘np’ for numpy linalg or scipy.sparse.linalg. Automatic choosing if solver is None depending on matrixtype.
verbose (bool [False]) – Be verbose.
- Returns
x – Solution vector
- Return type
GIMLI::RVector
Examples using pygimli.solver.linSolve
normH1¶
-
pygimli.solver.normH1(u, mat=None, mesh=None)[source]¶ Create (H1) norm for the finite element space.
- Parameters
u (iterable) – Node based value to compute the H1 norm for.
mat (Matrix) – Stiffness matrix.
mesh (GIMLI::Mesh) – Mesh with the FE space to generate S if necessary.
- Returns
ret – \(H1(u)\) norm.
- Return type
normL2¶
-
pygimli.solver.normL2(u, mat=None, mesh=None)[source]¶ Create Lebesgue (L2) norm for finite element space.
Find the L2 Norm for a solution for 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 \mathrm{d}\:x)^{1/2} \\ & \approx h (\sum |f(x)|^2 )^{1/2} \\ L2(u) = || u ||_{L^2} & = (\int |u|^2 \mathrm{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.
mat (Matrix) – Mass element matrix.
mesh (GIMLI::Mesh) – Mesh with the FE space to generate M if necessary.
- Returns
ret – \(L2(u)\) norm.
- Return type
Examples using pygimli.solver.normL2
parseArgPairToBoundaryArray¶
-
pygimli.solver.parseArgPairToBoundaryArray(pair, mesh)[source]¶ Parse boundary related pair argument to create a list of [ GIMLI::Boundary, value|callable ].
- Parameters
pair (tuple) –
[marker, arg]
[marker, [callable, *kwargs]]
[marker, [arg_x, arg_y, arg_z]]
[boundary, arg]
[‘*’, arg]
[node, arg]
[[marker, …], arg] (REMOVE ME because of bad design)
[[boundary,…], arg] (REMOVE ME because of bad design)
[marker, callable, *kwargs] (REMOVE ME because of bad design)
[[marker, …], callable, *kwargs] (REMOVE ME because of bad design)
arg will be parsed by
pygimli.solver.solver.generateBoundaryValueand distributed to each boundary. Callable functions will be executed at run time. ‘*’ will be interpreted as all boundary elements with one neighboring cellmesh (GIMLI::Mesh) – Used to find boundaries by marker.
- Returns
bc – [GIMLI::Boundary, value|callable]
- Return type
list()
parseArgToArray¶
-
pygimli.solver.parseArgToArray(arg, nDof, mesh=None, userData={})[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={})
Where MeshEntity is one of GIMLI::Cell , GIMLI::Node or GIMLI::Boundary depending 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 – Array of desired length filled with the appropriate values.
- Return type
GIMLI::RVector
Examples using pygimli.solver.parseArgToArray
parseArgToBoundaries¶
-
pygimli.solver.parseArgToBoundaries(args, mesh)[source]¶ Parse boundary related arguments to create a valid boundary value list: [ GIMLI::Boundary, value|callable ]
- Parameters
args (dict, float, callable) –
Dictionary is preferred (key=value|callable). If args is just a callable or float every outer boundary will be processed with args.
List pairs will be removed or not correct parsed for vector valued problems. Callable will be evaluated at runtime. See examples. Else see
pygimli.solver.solver.parseArgPairToBoundaryArraymesh (GIMLI::Mesh) – Used to find boundaries by marker
- Returns
boundaries – [ GIMLI::Boundary, value|callable ]
- Return type
list()
>>> # no need to import matplotlib. pygimli show does >>> import pygimli as pg >>> import pygimli.meshtools as mt >>> plc = mt.createWorld([0, 0], [1, -1], worldMarker=0) >>> ax, _ = pg.show(plc, boundaryMarker=True) >>> mesh = mt.createMesh(plc) >>> # all four outer boundaries get value = 1.0 >>> b = pg.solver.parseArgToBoundaries(1.0, mesh) >>> print(len(b)) 4 >>> # 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 >>> # Boundary values for vector valued problem >>> b = pg.solver.parseArgToBoundaries({1:[1.0, 1.0]}, mesh) >>> print(len(b), b[0][1]) 1 [1.0, 1.0] >>> # 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:4':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()
Examples using pygimli.solver.parseArgToBoundaries
parseDictKey¶
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 attributes to this default.
- Returns
att – Array of length mesh.cellCount()
- Return type
array
Examples using pygimli.solver.parseMapToCellArray
show¶
-
pygimli.solver.show(obj=None, data=None, **kwargs)[source]¶ Mesh and model visualization.
Syntactic sugar to show a obj with data. Forwards to a known visualization for obj. Typical is
pygimli.viewer.showMeshorpygimli.viewer.mayaview.showMesh3Dto 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
obj (obj) – obj can be so far. * GIMLI::Mesh or list of meshes * DataContainer * pg.core.Sparse[Map]Matrix
data (iterable) – Optionally data to visualize. See appropriate show function.
- Keyword Arguments
**kwargs –
Additional kwargs forward to appropriate show functions.
- axaxe [None]
Matplotlib axes object. Create a new if necessary.
- fitViewbool [True]
Scale x and y limits to match the view.
- Returns
Return the results from the showMesh functions. Usually the axe object*
and a colorbar.
See also
showMesh()
showSparseMatrix¶
-
pygimli.solver.showSparseMatrix(mat, full=False)[source]¶ Show the content of a sparse matrix.
- Parameters
mat (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
solveFiniteElements¶
-
pygimli.solver.solveFiniteElements(mesh, a=1.0, b=None, f=0.0, bc=None, times=None, userData={}, 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}}\\ \frac{\partial u}{\partial \mathbf{n}} & = \alpha(u_0-u)~~|~~\Gamma_{\text{Robin}}\end{split}\]for the scalar \(u(\mathbf{r}, t)\) or vector \(\mathbf(u)(\mathbf{r}, t)\) solution at each node of a given mesh. The Domain \(\Omega\) and the Boundary \(\Gamma\) are defined through the mesh with appropriate boundary marker.
Note, to ensure vector solution either set vector forces or at least on vector component boundary condition.
- Parameters
mesh (GIMLI::Mesh) – Mesh represents spatial discretization of the calculation domain
a (value | array | callable(cell, userData)) – Cell values of type float or complex can be scalar, anisotropic matrix or elastic tensor.
b (value | array | callable(cell, userData) [None]) – Cell values. None means the term is unused.
f (value | array(cells) | array(nodes) | callable(args, kwargs)) – force values, for vector fields use (n x dim) values.
bc (dict()) –
Dictionary of boundary conditions. Current supported boundary conditions by dictionary keys: ‘Dirichlet’, ‘Neumann’, ‘Robin’, ‘Node’.
The dictionary can contain multiple “key: Arg” Arg will be parsed by
pygimli.solver.solver.parseArgPairToBoundaryArrayIf 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.
- Keyword Arguments
**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-Nicolson, maybe instable - \(theta = 1\), implicit Euler
If unsure choose \(\theta = 0.5 + \epsilon\), which is probably stable.
- dynamic: bool [False]
Boundary conditions for time depending problems will be considered dynamic for each time step.
- 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.
- pureNeumann: bool [auto]
If set or detected automatic, we add the additional condition: \(\int_domain u dv = 0\) which makes elliptic problems well posed again.
- rhs: iterable
Pre assembled rhs. Will preferred on any f settings.
- 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 – Returns the solution u either 1,n array for stationary problems or for m,n array for m time steps
- Return type
array
>>> import pygimli as pg >>> from pygimli.meshtools import polytools as plc >>> from pygimli.viewer.mpl 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.0, 2: 1.0}, ... 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
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.
- 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.
Workspace (ws) –
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 –
bc : Boundary Conditions dictionary, see pg.solver
uB : Dirichlet boundary conditions DEPRECATED
duB : Neumann boundary conditions DEPRECATED
- Returns
u – Solution field for all time steps.
- Return type
ndarray(nTimes, nCells)
Examples using pygimli.solver.solveFiniteVolume
Classes¶
RungeKutta
-
class
pygimli.solver.RungeKutta(solver, verbose=False)[source]¶ Bases:
objectTODO 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]¶
-
WorkSpace
