Finite element basis for function spaces and test spaces

transformquadrature(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 points
• weights::Union{AbstractArray{Float64},Nothing} = nothing: optional array of weights to transform
• mapping::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

sausage_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

kosloff_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

hale_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

TensorMacroElementBasisFrom1D(
    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 dimension
• numberquadraturepoints1d::Int: number of quadrature points in 1 dimension
• numbercomponents::Int: number of components
• dimension::Int: dimension of basis
• numberelements1d::Int: number of elements in macro-element
• basis1dmicro::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

TensorHProlongationBasis(
    coarsenodes1d,
    finenodes1d,
    numbercomponents,
    dimension,
    numberfineelements1d,
)

Tensor product h-prolongation basis

Arguments:

• coarsenodes1d::AbstractArray{Float64,1}: coarse grid node coordinates in 1 dimension
• finenodes1d::AbstractArray{Float64,1}: fine grid node coordinates in 1 dimension
• numbercomponents::Int: number of components
• dimension::Int: dimension of basis
• numberfineelements1d::Int: number of fine grid elements

Returns:

• H1 Lagrange tensor product h-prolongation basis object

buildinterpolationandgradient(nodes1d, quadraturepoints1d)

Build one dimensional interpolation and gradient matrices, from Fornberg 1998

Arguments:

• nodes1d::AbstractArray{Float64}: 1d basis nodes
• quadraturepoints1d::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

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

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

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

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

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

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

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

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

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

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

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

getdXdxgradient(basis, mesh)

Get gradient adjusted for mesh stretching

Arguments:

• basis::TensorBasis: basis to compute gradient
• mesh::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

getdxdXquadratureweights(basis, mesh)

Get quadrature weights adjusted for mesh stretching

Arguments:

• basis::TensorBasis: basis to compute quadratureweights
• mesh::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

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

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