Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
LagrangePolynomials = "fb83b288-ce1a-45f4-939e-c34c34051e49"

[sources]
LagrangePolynomials = {rev = "main", url = "https://github.com/moment-kinetics/LagrangePolynomials"}
LagrangePolynomials = {rev = "main", url = "https://github.com/moment-kinetics/LagrangePolynomials.jl.git"}

[compat]
LagrangePolynomials = "0.1.0"
Expand Down
26 changes: 13 additions & 13 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,13 +10,13 @@ This package computes matrices of the following form.
A_{ij} = \int^{x_u}_{x_l} P_i(x) Q_j(x) x^p d x,
```
where $`P_i`$ and $`Q_i`$ are from the set of Lagrange polynomials basis
functions (and their first derivatives) used for
functions (and their first derivatives) used for
$`C^0`$ continuous Galerkin
finite element models, $`x_u`$ and $`x_l`$
are the upper and lower limits of the element in the
coordinate $`x`$, and $`p`$ is an integer.

The basis functions we use here are
The basis functions we use here are
```math
\Phi(x) = \Theta(x-x_l)\Theta(x_u - x)l_j(z(x)),
```
Expand All @@ -33,7 +33,7 @@ is required for the problem of interest.
To compute the matrix $`A_{ij}`$, we use the reference coordinate
$`z`$ to write
```math
A_{ij} = \int^{1}_{-1} P_i(z) Q_j(z) (s z + c)^p d z.
A_{ij} = \int^{1}_{-1} P_i(z) Q_j(z) (s z + c)^p d z.
```
We can also compute nonlinear matrices which have the form
```math
Expand All @@ -45,8 +45,8 @@ N_{ijk} = \int^{x_u}_{x_l} P_i(x) Q_j(x) R_k(x) x^p d x = \int^{1}_{-1} P_i(z) Q
To create these matrices, the user must supply the key grid
information for the element of interest, i.e., the set
$`\{z_j\}`$, and the factors $`s`$ and $`c`$.
This is done in the following function call, defining a
type `element_coordinates`. We use [FastGaussQuadrature](https://juliaapproximation.github.io/FastGaussQuadrature.jl/stable/)
This is done in the following function call, defining a
type `ElementCoordinates`. We use [FastGaussQuadrature](https://juliaapproximation.github.io/FastGaussQuadrature.jl/stable/)
to make a set of reference nodes, but any set of nodes with good
interpolation properties would be adequate.
```
Expand All @@ -57,23 +57,23 @@ x, w = gausslobatto(ngrid)
scale = 1.0 # change these to match the size and position of elements in the problem
shift = 0.0
# the coordinate information
coordinate = element_coordinates(x,scale,shift)
coordinate = ElementCoordinates(x,scale,shift)
```
Once an instance of `element_coordinates` is available,
Once an instance of `ElementCoordinates` is available,
we can compute the higher-order finite element matrices
as follows. We use the enums
`lagrange_x`, and `d_lagrange_dx` and the integer `p` to choose the form of the desired matrix. For a mass matrix, with
Jacobian $`1`$ for example,
we can use the following commands.
```
using FiniteElementMatrices
coordinate = element_coordinates(x,scale,shift)
M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
coordinate = ElementCoordinates(x,scale,shift)
M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
```
For a stiffness matrix describing a second derivative, we could
use the following.
```
K = -finite_element_matrix(d_lagrange_dx,d_lagrange_dx,0,coordinate)
K = -finite_element_matrix(d_lagrange_dx,d_lagrange_dx,0,coordinate)
```
The second derivative of a function $`f`$ can be found by solving the system
```math
Expand All @@ -86,8 +86,8 @@ in $`x`$.
For a system with Jacobian $`x`$ (e.g., cylindrical polar coordinates), we could form the appropriate matrices to evaluate
the 1D, radial Laplacian via
```
coordinate = element_coordinates(x,scale,shift)
M = finite_element_matrix(lagrange_x,lagrange_x,1,coordinate)
coordinate = ElementCoordinates(x,scale,shift)
M = finite_element_matrix(lagrange_x,lagrange_x,1,coordinate)
K = -finite_element_matrix(d_lagrange_dx,d_lagrange_dx,1,coordinate)
```
i.e.,
Expand All @@ -96,7 +96,7 @@ i.e.,
```
Finally, we can form a nonlinear stiffness matrix by inserting an extra enum argument, as follows.
```
N = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,1,coordinate)
N = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,1,coordinate)
```

# Examples
Expand Down
32 changes: 16 additions & 16 deletions src/FiniteElementMatrices.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,24 +4,24 @@ module FiniteElementMatrices

using LagrangePolynomials: lagrange_poly,
lagrange_poly_derivative,
lagrange_poly_data
LagrangePolyData
using FastGaussQuadrature: gausslegendre

export lagrange_x,
d_lagrange_dx,
finite_element_matrix,
element_coordinates
ElementCoordinates

@enum lagrange_function_type begin
@enum LagrangeFunctionType begin
lagrange_x
d_lagrange_dx
end

struct element_coordinates
struct ElementCoordinates
# precomputed data for calculating
# the Lagrange polynomials, including
# the reference x nodes
lpoly_data::lagrange_poly_data
lpoly_data::LagrangePolyData
# scale factor s for transform
# v = s x + c
# from reference grid x on [-1,1] (or (-1,1], [-1,1), (-1,1))
Expand All @@ -37,19 +37,19 @@ struct element_coordinates
`shift` : normalisation factor, such that the physical coordinate
`v = `scale * x_nodes + shift` in this element.
"""
function element_coordinates(x_nodes::AbstractArray{Float64,1},
function ElementCoordinates(x_nodes::AbstractArray{Float64,1},
scale::Float64,
shift::Float64)
# initialise and save Lagrange polynomial
# data for this 1D element.
lpoly_data = lagrange_poly_data(x_nodes)
lpoly_data = LagrangePolyData(x_nodes)
return new(lpoly_data,
scale,
shift)
end
end

function select_lagrange_function(fn_type::lagrange_function_type)
function select_lagrange_function(fn_type::LagrangeFunctionType)
if fn_type == lagrange_x
func = lagrange_poly
elseif fn_type == d_lagrange_dx
Expand All @@ -58,7 +58,7 @@ function select_lagrange_function(fn_type::lagrange_function_type)
return func
end

function select_lagrange_prefactor(fn_type::lagrange_function_type,
function select_lagrange_prefactor(fn_type::LagrangeFunctionType,
scale::Float64)
if fn_type == lagrange_x
prefactor = 1.0
Expand All @@ -82,10 +82,10 @@ function get_polyz(power::Int64,
end

function finite_element_matrix(
fn1_type::lagrange_function_type,
fn2_type::lagrange_function_type,
fn1_type::LagrangeFunctionType,
fn2_type::LagrangeFunctionType,
power::Int64,
coordinate::element_coordinates
coordinate::ElementCoordinates
)
lpoly_data = coordinate.lpoly_data
ngrid = length(coordinate.lpoly_data.x_nodes)
Expand Down Expand Up @@ -124,11 +124,11 @@ function finite_element_matrix(
end

function finite_element_matrix(
fn1_type::lagrange_function_type,
fn2_type::lagrange_function_type,
fn3_type::lagrange_function_type,
fn1_type::LagrangeFunctionType,
fn2_type::LagrangeFunctionType,
fn3_type::LagrangeFunctionType,
power::Int64,
coordinate::element_coordinates
coordinate::ElementCoordinates
)
lpoly_data = coordinate.lpoly_data
ngrid = length(coordinate.lpoly_data.x_nodes)
Expand Down
10 changes: 5 additions & 5 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ using Test: @testset, @test
using FiniteElementMatrices: lagrange_x,
d_lagrange_dx,
finite_element_matrix,
element_coordinates
ElementCoordinates
using FastGaussQuadrature: gausslobatto, gaussradau, gausslegendre
using LinearAlgebra: lu, ldiv!, mul!

Expand Down Expand Up @@ -45,7 +45,7 @@ function test_first_derivative(;nodes::node_type=GLL,
x = reference_nodes(nodes,ngrid)
scale = 0.5*(y1-y0)
shift = 0.5*(y0+y1)
coordinate = element_coordinates(x,scale,shift)
coordinate = ElementCoordinates(x,scale,shift)
M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
P = finite_element_matrix(lagrange_x,d_lagrange_dx,pp,coordinate)
S = finite_element_matrix(d_lagrange_dx,lagrange_x,pp,coordinate)
Expand Down Expand Up @@ -86,7 +86,7 @@ function test_first_derivative_nonlinear_operators(
x = reference_nodes(nodes,ngrid)
scale = 0.5*(y1-y0)
shift = 0.5*(y0+y1)
coordinate = element_coordinates(x,scale,shift)
coordinate = ElementCoordinates(x,scale,shift)
M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
Y100 = finite_element_matrix(d_lagrange_dx,lagrange_x,lagrange_x,pp,coordinate)
Y010 = finite_element_matrix(lagrange_x,d_lagrange_dx,lagrange_x,pp,coordinate)
Expand Down Expand Up @@ -131,7 +131,7 @@ function test_second_derivative(;nodes::node_type=GLL,
x = reference_nodes(nodes,ngrid)
scale = 0.5*(y1-y0)
shift = 0.5*(y0+y1)
coordinate = element_coordinates(x,scale,shift)
coordinate = ElementCoordinates(x,scale,shift)
M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
S = -finite_element_matrix(d_lagrange_dx,lagrange_x,0,coordinate)
K = -finite_element_matrix(d_lagrange_dx,d_lagrange_dx,0,coordinate)
Expand Down Expand Up @@ -189,7 +189,7 @@ function test_nonlinear_operators(;nodes::node_type=GLL,
x = reference_nodes(nodes,ngrid)
scale = 0.5*(y1-y0)
shift = 0.5*(y0+y1)
coordinate = element_coordinates(x,scale,shift)
coordinate = ElementCoordinates(x,scale,shift)
M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
Y000 = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,0,coordinate)
Y100 = finite_element_matrix(d_lagrange_dx,lagrange_x,lagrange_x,0,coordinate)
Expand Down