Discretization

AbstractDiscretization

Discrete function spaces for the potential u and the conductivity σ over one mesh. u and σ live in different finite element spaces on the same grid (e.g. P1/P0, but any pair is allowed). Carries the boundary description and everything needed to assemble mass, stiffness and weighted stiffness matrices. After adaptive mesh refinement a new discretization is built.

Interface: ndofs_u, ndofs_σ.

FerriteDiscretization(grid; ip_u, ip_σ, boundary, qr_order)

Finite element spaces for the potential u (interpolation ip_u, default linear Lagrange) and the conductivity σ (interpolation ip_σ, default piecewise constant DiscontinuousLagrange{…,0}) on the same Ferrite grid.

• boundary: boundary facets (a collection of FacetIndex or the name of a facet set). Default: all facets that belong to exactly one cell.
• qr_order: quadrature order. Default: exact for the weighted stiffness matrix ∫ σ ∇φᵢ⋅∇φⱼ and the mass matrices on affine cells.

On a non-conforming grid (hanging nodes, see AdaptiveMesh) the u space is the conforming subspace u_full = C u; only linear Lagrange ip_u and discontinuous ip_σ are supported there.

Fields: grid, dh_u, dh_σ, ip_u, ip_σ, cv_u, cv_σ (cell values on one common quadrature rule), fv_u (facet values), boundary_facets, boundary_dofs (u dofs on the boundary), interior_facets ((cell₁, facet₁, cell₂) for every facet shared by two cells, including the pieces of coarse–fine interfaces), C_u (conformity constraints ndofs(dh_u) × ndofs_u, nothing on conforming grids).

ndofs_u(disc)

Number of degrees of freedom of the potential (the free ones on non-conforming meshes). Back end contract.

ndofs_σ(disc)

Number of degrees of freedom of the conductivity. Back end contract.

FEMatrices(disc)

Assembled matrices of a discretization: M_u, K_u (mass and stiffness of the u space), M_Γ (boundary mass of the u space on the boundary), M_σ, K_σ (mass and stiffness of the σ space; K_σ = 0 for piecewise constants) and M_σ_fac (Cholesky factorisation of M_σ, used for L² projections and L² gradients). The constructor for a discretization is part of the back end contract.

assemble_mass!(M, dh, cv)
assemble_mass(dh, cv)

Mass matrix ∫ φᵢ φⱼ dΩ of the (single-field) DofHandler dh with cell values cv.

assemble_stiffness!(K, dh, cv)
assemble_stiffness(dh, cv)

Stiffness matrix ∫ ∇φᵢ⋅∇φⱼ dΩ.

assemble_boundary_mass!(M, dh, fv, facets)
assemble_boundary_mass(dh, fv, facets)

Boundary mass matrix ∫_Γ φᵢ φⱼ ds over the facets facets, with facet values fv.

assemble_boundary_load!(f, dh, fv, facets)

Boundary load vector fᵢ = ∫_Γ φᵢ ds over facets (added to f).

ConductivityTensor(disc; pattern = <u-space pattern>, to_device = identity)

Sparse tensor T (nnz(pattern) × n_σ) with nzval(L(σ)) = T σ for the weighted stiffness matrix L(σ) = ∫ σ ∇φᵢ⋅∇φⱼ in the storage order of pattern. pattern may be larger than the u–u block (e.g. the augmented matrix of the complete electrode model), as long as the u dofs come first; entries outside the u–u block get zero rows. The constructor for a discretization is part of the back end contract.

Fields: pattern (host SparseMatrixCSC), T and Tt (T and Tᵀ, on the device if to_device is given), rows, cols (row/column of each stored entry) and a buffer w.

See assemble_weighted_stiffness! and tensor_gradient!.

assemble_weighted_stiffness!(L, disc, σ)
assemble_weighted_stiffness(disc, σ)

Weighted stiffness matrix ∫ σ ∇φᵢ⋅∇φⱼ dΩ by a classical element loop, with σ given by its coefficients in the σ space of disc. The ! version fills a matrix with the pattern of allocate_matrix(disc.dh_u) (all dofs, before conformity constraints); the allocating version returns the matrix of the u space of disc (condensed on non-conforming grids). For repeated assembly on a fixed mesh, the ConductivityTensor method is a single sparse matrix-vector product.

assemble_weighted_stiffness!(L, ct::ConductivityTensor, σ; A₀ = nothing)

nzval(L) ← T σ (+ A₀): the weighted stiffness matrix for the conductivity coefficients σ by one sparse matrix-vector product. L must have the sparsity pattern ct.pattern (e.g. copy(ct.pattern)). A₀ is an optional constant part of the stored values (e.g. the contact impedance terms of the complete electrode model).

weighted_stiffness_values!(nzval, ct::ConductivityTensor, σ; A₀ = nothing)

Stored values nzval ← T σ (+ A₀) of the weighted stiffness matrix, on the device of ct.

pair_products!(w, ct::ConductivityTensor, Λ, U)

wₖ = Σₛ Λ[rowₖ, s] U[colₖ, s] for every stored entry k of the pattern (Λ, U: vectors or n × s blocks with n ≥ the size of the pattern's u–u block). KernelAbstractions kernel: runs multithreaded on the CPU and on every GPU backend.

tensor_gradient!(g, ct::ConductivityTensor, Λ, U; α = 1, β = 0)

g ← α Σₛ ∂/∂σ (λₛᵀ L(σ) uₛ) + β g, i.e. gₐ = α Σₛ λₛᵀ Lₐ uₛ + β gₐ = α Σₛ ∫ ψₐ ∇uₛ⋅∇λₛ + β gₐ for the columns λₛ, uₛ of Λ, U. This is the conductivity gradient of the adjoint state method (with α = -1) and of the Kohn–Vogelius functional (with Λ = U). Two parallel passes: pair_products! and one sparse Tᵀ w product.

AbstractRieszMap

Representation of the gradient: the coefficient gradient ∂J/∂σₐ of the discrete functional (discretize-then-optimize, CoefficientGradient) or its Riesz representative in a function space, e.g. the L² gradient M_σ⁻¹ ∂J/∂σ (L2Gradient).

CoefficientGradient()

The gradient as the vector of partial derivatives ∂J/∂σₐ (a dual vector). This is the exact gradient of the discrete objective (discretize-then-optimize); it depends on the mesh.

L2Gradient(mats::FEMatrices)
L2Gradient(M_σ_factorization)

The L² Riesz representative M_σ⁻¹ ∂J/∂σ: the L² projection of ∫ ψₐ ∇u⋅∇λ onto the σ space (optimize-then-discretize with projection). It approximates the continuous gradient -∇u⋅∇λ independently of the mesh. For piecewise constant σ, M_σ is diagonal with the cell areas.

riesz_map!(g, R::AbstractRieszMap)
riesz_map(R, g)

Apply the Riesz map R to the coefficient gradient g (in place / out of place).

interpolate_function(disc, f; field = :σ)

Coefficients of the interpolant of the function f(x) in the σ space (field = :σ) or the u space (field = :u). Back end contract.

l2_project(disc, f; field = :σ, kwargs...)

Coefficients of the L² projection of the function f(x) onto the σ or u space. Back end contract.

fe_inner(disc, a, b; field = :σ, kind = :L2, mats = nothing)

Inner product of two coefficient vectors as functions: :L2, :H1semi or :H1. Back end contract.

fe_norm(disc, a; kwargs...)

Norm √fe_inner(disc, a, a; kwargs...). Back end contract.

total_variation(disc, σ; ε = 0)
total_variation!(g, disc, σ; ε = 0)

Total variation ∫ √(|∇σ|² + ε²) of a conductivity (for piecewise constants: the facet jumps), and its coefficient gradient (in place). Back end contract.