Finite Element Basis
LFAToolkit.AbstractBasis
— TypeFinite element basis for function spaces and test spaces
LFAToolkit.transformquadrature
— Functiontransformquadrature(points, weights = nothing, mapping = nothing)
Transformed quadrature by applying a smooth mapping = (g, gprime) from the original domain.
Arguments:
points::AbstractArray{Float64}
: array of quadrature pointsweights::Union{AbstractArray{Float64},Nothing} = nothing
: optional array of weights to transformmapping::Union{Tuple{Function,Function},Nothing} = nothing
: choice of conformal map
Returns:
- transformed quadrature points and weights
Example:
# generate transformed quadrature points, weights with choice of conformal map
quadraturepoints = [
-√(5 + 2 * √(10 / 7)) / 3,
-√(5 - 2 * √(10 / 7)) / 3,
0.0,
√(5 - 2 * √(10 / 7)) / 3,
√(5 + 2 * √(10 / 7)) / 3,
];
quadratureweights = [
(322 - 13 * √(70)) / 900,
(322 + 13 * √(70)) / 900,
128 / 225,
(322 + 13 * √(70)) / 900,
(322 - 13 * √(70)) / 900,
];
mapping = sausage_transformation(9);
mappedpoints, mappedweights =
LFAToolkit.transformquadrature(quadraturepoints, quadratureweights, mapping);
# verify:
weightsum = sum(mappedweights);
@assert weightsum ≈ 2.0
# output
LFAToolkit.sausage_transformation
— Functionsausage_transformation(d)
Compute the conformal mapping of Gauss ellipses to sausage_transformations using a truncated Taylor expansion of arcsin(s). See Figure 4.1 of Hale and Trefethen (2008).
Arguments:
d::Int
: polynomial degree of truncated Taylor series expansion of arcsin(s).
Returns:
- conformal mapping
Example:
# sausage_transformation conformal map
g, gprime = LFAToolkit.sausage_transformation(9);
# verify
truegonrange = [
-0.9999999999999999,
-0.39765215163451934,
0.0,
0.39765215163451934,
0.9999999999999999,
];
@assert g.(LinRange(-1, 1, 5)) ≈ truegonrange
# output
LFAToolkit.kosloff_talezer_transformation
— Functionkosloff_talezer_transformation(α)
Compute the Kosloff and Tal-Ezer conformal map derived from the inverse sine function
Arguments:
α::Float64
: polynomial degree of truncated Taylor series expansion of arcsin(s).
Returns:
- conformal mapping
Example:
# kosloff_talezer_transformation conformal map
g, gprime = LFAToolkit.kosloff_talezer_transformation(0.95);
# verify
truegonrange = [-1.0, -0.39494881426787537, 0.0, 0.39494881426787537, 1.0];
@assert g.(LinRange(-1, 1, 5)) ≈ truegonrange
# output
LFAToolkit.hale_trefethen_strip_transformation
— Functionhale_trefethen_strip_transformation(ρ)
Compute the Hale and Trefethen strip transformation
Arguments:
ρ::Float64
: sum of the semiminor and semimajor axis
Returns:
- conformal mapping
Example:
# hale_trefethen_strip_transformation conformal map
g, gprime = hale_trefethen_strip_transformation(1.4);
# verify
truegonrange = [-1.0, -0.36812132798370184, 0.0, 0.36812132798370184, 1.0];
@assert g.(LinRange(-1, 1, 5)) ≈ truegonrange
# output
LFAToolkit.TensorMacroElementBasisFrom1D
— FunctionTensorMacroElementBasisFrom1D(
numbernodes1d,
numberquadraturepoints1d,
numbercomponents,
dimension,
numberelements1d,
basis1dmicro;
overlapquadraturepoints = false,
)
Tensor product macro-element basis from 1d single element tensor product basis
Arguments:
numbernodes1d::Int
: number of basis nodes in 1 dimensionnumberquadraturepoints1d::Int
: number of quadrature points in 1 dimensionnumbercomponents::Int
: number of componentsdimension::Int
: dimension of basisnumberelements1d::Int
: number of elements in macro-elementbasis1dmicro::TensorBasis
: 1d micro element basis to replicate
Keyword Arguments:
overlapquadraturepoints::Bool = false
: overlap quadrature points between elements, for prolongation
Returns:
- tensor product macro-element basis object
LFAToolkit.TensorHProlongationBasis
— FunctionTensorHProlongationBasis(
coarsenodes1d,
finenodes1d,
numbercomponents,
dimension,
numberfineelements1d,
)
Tensor product h-prolongation basis
Arguments:
coarsenodes1d::AbstractArray{Float64,1}
: coarse grid node coordinates in 1 dimensionfinenodes1d::AbstractArray{Float64,1}
: fine grid node coordinates in 1 dimensionnumbercomponents::Int
: number of componentsdimension::Int
: dimension of basisnumberfineelements1d::Int
: number of fine grid elements
Returns:
- H1 Lagrange tensor product h-prolongation basis object
LFAToolkit.buildinterpolationandgradient
— Functionbuildinterpolationandgradient(nodes1d, quadraturepoints1d)
Build one dimensional interpolation and gradient matrices, from Fornberg 1998
Arguments:
nodes1d::AbstractArray{Float64}
: 1d basis nodesquadraturepoints1d::AbstractArray{Float64}
: 1d basis quadrature points
Returns:
- one dimensional interpolation and gradient matrices
Example:
# nodes, quadrature points, and weights
numberquadraturepoints = 4;
nodes = [-1.0, 0.0, 1.0];
quadraturepoints = [
-√(3 / 7 + 2 / 7 * √(6 / 5)),
-√(3 / 7 - 2 / 7 * √(6 / 5)),
√(3 / 7 - 2 / 7 * √(6 / 5)),
√(3 / 7 + 2 / 7 * √(6 / 5)),
];
# build interpolation, gradient matrices
interpolation, gradient = LFAToolkit.buildinterpolationandgradient(nodes, quadraturepoints);
# verify
for i = 1:numberquadraturepoints
total = sum(interpolation[i, :])
@assert total ≈ 1.0
total = sum(gradient[i, :])
@assert abs(total) < 1e-14
end
# output
LFAToolkit.getnumbernodes
— Functiongetnumbernodes(basis)
Get the number of nodes for the basis
Returns:
- integer number of basis nodes
Example:
# setup basis
basis = TensorH1LagrangeBasis(4, 3, 2, 2);
# verify number of nodes
@assert basis.numbernodes == 4^2
# output
LFAToolkit.getnodes
— Functiongetnodes(basis)
Get nodes for basis
Returns:
- basis nodes array
Example:
# test for all supported dimensions
for dimension = 1:3
# setup basis
basis = TensorH1LagrangeBasis(4, 3, 2, dimension)
# get true nodes
truenodes1d = [-1, -√(1 / 5), √(1 / 5), 1]
truenodes = []
if dimension == 1
truenodes = truenodes1d
elseif dimension == 2
truenodes = transpose(hcat([[[x, y] for x in truenodes1d, y in truenodes1d]...]...))
elseif dimension == 3
truenodes = transpose(
hcat(
[
[[x, y, z] for x in truenodes1d, y in truenodes1d, z in truenodes1d]...,
]...,
),
)
end
# verify basis nodes
@assert basis.nodes ≈ truenodes
end
# output
LFAToolkit.getnumberquadraturepoints
— Functiongetnumberquadraturepoints(basis)
Get the number of quadrature points for the basis
Returns:
- integer number of basis quadrature points
Example:
# setup basis
basis = TensorH1LagrangeBasis(4, 3, 2, 2);
# verify number of quadrature points
@assert basis.numberquadraturepoints == 3^2
# output
LFAToolkit.getquadraturepoints
— Functiongetquadraturepoints(basis)
Get quadrature points for basis
Returns:
- basis quadrature points array
Example:
# test for all supported dimensions
for dimension = 1:3
# setup basis
basis = TensorH1LagrangeBasis(4, 3, 2, dimension)
# get true quadrature points
truequadraturepoints1d = [-√(3 / 5), 0, √(3 / 5)]
truequadraturepoints = []
if dimension == 1
truequadraturepoints = truequadraturepoints1d
elseif dimension == 2
truequadraturepoints = transpose(
hcat(
[
[
[x, y] for x in truequadraturepoints1d, y in truequadraturepoints1d
]...,
]...,
),
)
elseif dimension == 3
truequadraturepoints = transpose(
hcat(
[
[
[x, y, z] for x in truequadraturepoints1d,
y in truequadraturepoints1d, z in truequadraturepoints1d
]...,
]...,
),
)
end
# verify quadrature points
@assert basis.quadraturepoints ≈ truequadraturepoints
end
# output
LFAToolkit.getquadratureweights
— Functiongetquadratureweights(basis)
Get full quadrature weights vector for basis
Returns:
- basis quadrature weights vector
Example:
# test for all supported dimensions
for dimension = 1:3
# setup basis
basis = TensorH1LagrangeBasis(4, 3, 2, dimension)
# get true quadratrue weights
trueweights1d = [5 / 9, 8 / 9, 5 / 9]
trueweights = []
if dimension == 1
trueweights = trueweights1d
elseif dimension == 2
trueweights = kron(trueweights1d, trueweights1d)
elseif dimension == 3
trueweights = kron(trueweights1d, trueweights1d, trueweights1d)
end
# verify
@assert basis.quadratureweights ≈ trueweights
end
# output
LFAToolkit.getinterpolation
— Functiongetinterpolation(basis)
Get full interpolation matrix for basis
Returns:
- basis interpolation matrix
Example:
# test for all supported dimensions
for dimension = 1:3
# setup basis
basis = TensorH1LagrangeBasis(4, 3, 2, dimension)
# get basis interpolation matrix
interpolation = basis.interpolation
# verify
for i = 1:3^dimension
total = sum(interpolation[i, :])
@assert total ≈ 1.0
end
end
# output
LFAToolkit.getgradient
— Functiongetgradient(basis)
Get full gradient matrix for basis
Returns:
- basis gradient matrix
Example:
# test for all supported dimensions
for dimension = 1:3
# setup basis
basis = TensorH1LagrangeBasis(4, 3, 2, dimension)
# get basis gradient matrix
gradient = basis.gradient
# verify
for i = 1:dimension*3^dimension
total = sum(gradient[i, :])
@assert abs(total) < 1e-14
end
end
# output
LFAToolkit.getnumbermodes
— Functiongetnumbermodes(basis)
Get number of modes for basis
Returns:
- number of modes for basis
Example:
# test for all supported dimensions
for dimension = 1:3
# setup basis
basis = TensorH1LagrangeBasis(4, 3, 2, dimension)
# verify number of basis modes
@assert basis.numbermodes == 2 * 3^dimension
end
# output
LFAToolkit.getmodemap
— Functiongetmodemap(basis)
Get mode mapping vector for basis
Returns:
- basis mode map vector
Example:
# test for all supported dimensions
for dimension = 1:3
# setup basis
basis = TensorH1LagrangeBasis(4, 3, 1, dimension)
# get true mode map
truemodemap1d = [1, 2, 3, 1]
truemodemap = []
if dimension == 1
truemodemap = truemodemap1d
elseif dimension == 2
truemodemap = [[i + (j - 1) * 3 for i in truemodemap1d, j in truemodemap1d]...]
elseif dimension == 3
truemodemap = [
[
i + (j - 1) * 3 + (k - 1) * 3^2 for i in truemodemap1d,
j in truemodemap1d, k in truemodemap1d
]...,
]
end
# verify
@assert basis.modemap == truemodemap
end
# output
LFAToolkit.getnumberelements
— Functiongetnumberelements(basis)
Get the number of elements for the basis
Returns:
- integer number of basis micro-elements
Example:
# get number of nodes for basis
basis = TensorH1LagrangeMacroBasis(4, 4, 1, 2, 2);
# verify number of elements
@assert basis.numberelements == 2^2
# output
LFAToolkit.getdiagonal
— Functiongetdiagonal(operator)
Compute or retrieve the symbol matrix diagonal for an operator
Returns:
- symbol matrix diagonal for the operator
Example:
# setup
mesh = Mesh1D(1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);
# verify diagonal
@assert diffusion.diagonal ≈ [7/6 0; 0 4/3]
# output
LFAToolkit.getdXdxgradient
— FunctiongetdXdxgradient(basis, mesh)
Get gradient adjusted for mesh stretching
Arguments:
basis::TensorBasis
: basis to compute gradientmesh::Mesh
: mesh to compute gradient
Returns:
- gradient matrix multiplied by change of coordinates adjoint
Example:
for dimension = 1:3
# mesh
mesh = []
if dimension == 1
mesh = Mesh1D(2.0)
elseif dimension == 2
mesh = Mesh2D(2.0, 3.0)
elseif dimension == 3
mesh = Mesh3D(2.0, 3.0, 4.0)
end
# basis
basis = TensorH1LagrangeBasis(4, 3, 1, dimension)
# get gradient on mesh
gradient = LFAToolkit.getdXdxgradient(basis, mesh)
# verify
nodes = basis.nodes
linearfunction = []
truegradient = []
if dimension == 1
linearfunction = nodes
truegradient = [1 / 2 * ones(basis.numberquadraturepoints)...]
elseif dimension == 2
linearfunction = (nodes[:, 1] + nodes[:, 2])
truegradient = [
1 / 2 * ones(basis.numberquadraturepoints)...
1 / 3 * ones(basis.numberquadraturepoints)...
]
elseif dimension == 3
linearfunction = (nodes[:, 1] + nodes[:, 2] + nodes[:, 3])
truegradient = [
1 / 2 * ones(basis.numberquadraturepoints)...
1 / 3 * ones(basis.numberquadraturepoints)...
1 / 4 * ones(basis.numberquadraturepoints)...
]
end
@assert gradient * linearfunction ≈ truegradient
end
# output
LFAToolkit.getdxdXquadratureweights
— FunctiongetdxdXquadratureweights(basis, mesh)
Get quadrature weights adjusted for mesh stretching
Arguments:
basis::TensorBasis
: basis to compute quadratureweightsmesh::Mesh
: mesh to compute quadratureweights
Returns:
- quadrature weights multiplied by change of coordinates adjoint
Example:
mesh = Mesh1D(2.0)
# basis
basis = TensorH1LagrangeBasis(4, 3, 1, 1);
# get gradient on mesh
weights = LFAToolkit.getdxdXquadratureweights(basis, mesh);
# verify
@assert basis.quadratureweights * mesh.volume / basis.volume ≈ weights
# output
LFAToolkit.getprimalnodes
— Methodgetprimalnodes(basis)
Get primal nodes for basis
Returns:
- basis primal nodes vector
Example:
# test for all supported dimensions
for dimension = 1:3
# setup basis
p = 4
basis = TensorH1LagrangeBasis(p, p, 1, dimension)
# get true primal modes
trueprimalnodes = []
if dimension == 1
trueprimalnodes = [1, p]
elseif dimension == 2
trueprimalnodes = [1, p, p^2 - p + 1, p^2]
elseif dimension == 3
trueprimalnodes =
[1, p, p^2 - p + 1, p^2, p^3 - p^2 + 1, p^3 - p^2 + p, p^3 - p + 1, p^3]
end
# verify
@assert basis.primalnodes == trueprimalnodes
end
# output
LFAToolkit.getinterfacenodes
— Methodgetinterfacenodes(basis)
Get interface nodes for basis
Returns:
- basis primal nodes vector
Example:
# test for all supported dimensions
for dimension = 1:3
# setup basis
p = 4
basis = TensorH1LagrangeBasis(p, p, 1, dimension)
# get true interface nodes
trueinterfacenodes = []
if dimension == 1
trueinterfacenodes = [1, p]
elseif dimension == 2
trueinterfacenodes = [1:p..., p^2-p+1:p^2...]
for i = 1:p-2
push!(trueinterfacenodes, i * p + 1)
push!(trueinterfacenodes, (i + 1) * p)
end
elseif dimension == 3
trueinterfacenodes = [1:p^2..., p^3-p^2+1:p^3...]
for i = 1:p-2
push!(trueinterfacenodes, i*p^2+1:i*p^2+p...)
push!(trueinterfacenodes, (i+1)*p^2-p+1:(i+1)*p^2...)
push!(trueinterfacenodes, i*p^2+p+1:p:(i+1)*p^2-2*p+1...)
push!(trueinterfacenodes, i*p^2+2*p:p:(i+1)*p^2-p...)
end
end
trueinterfacenodes = sort(trueinterfacenodes)
# verify
@assert basis.interfacenodes == trueinterfacenodes
end
# output