From 54fa1133b235653228eb4e904f14aa1d43f21dad Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Sun, 28 Jun 2026 14:35:39 -0400 Subject: [PATCH] reworking a lot of internals to not require dof coordinates. --- Project.toml | 2 +- src/{AbstractTypes.jl => AbstractElements.jl} | 135 +---- src/PolynomialTypes.jl | 42 ++ src/QuadratureRules.jl | 568 ++++++++++++++++++ src/ReferenceFE.jl | 218 ++----- src/ReferenceFiniteElements.jl | 8 +- src/elements/Edge.jl | 86 +-- src/elements/Hex.jl | 341 +++-------- src/elements/Quad.jl | 202 +++---- src/elements/Tet.jl | 201 +------ src/elements/Tri.jl | 219 +------ src/elements/Vertex.jl | 6 +- test/runtests.jl | 198 +++--- 13 files changed, 992 insertions(+), 1234 deletions(-) rename src/{AbstractTypes.jl => AbstractElements.jl} (68%) create mode 100644 src/PolynomialTypes.jl create mode 100644 src/QuadratureRules.jl diff --git a/Project.toml b/Project.toml index 2f9214f..e65ff42 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ReferenceFiniteElements" uuid = "6dc62d09-f8eb-43fd-9672-074e490a997f" -version = "0.14.5" +version = "0.14.6" authors = ["Craig M. Hamel and contributors"] [deps] diff --git a/src/AbstractTypes.jl b/src/AbstractElements.jl similarity index 68% rename from src/AbstractTypes.jl rename to src/AbstractElements.jl index 820433a..6ed00ae 100644 --- a/src/AbstractTypes.jl +++ b/src/AbstractElements.jl @@ -1,102 +1,3 @@ -""" -$(TYPEDEF) -""" -abstract type AbstractPolynomialType end - -# TODO move somewhere else eventually -""" -$(TYPEDEF) -""" -struct Hermite <: AbstractPolynomialType -end -""" -$(TYPEDEF) -""" -struct Lagrange <: AbstractPolynomialType -end -# special type for Vertex -""" -$(TYPEDEF) -""" -struct NoInterpolation <: AbstractPolynomialType -end -""" -$(TYPEDEF) -""" -struct RaviartThomas <: AbstractPolynomialType -end -""" -$(TYPEDEF) -""" -struct Serendipity <: AbstractPolynomialType -end - -""" -$(TYPEDEF) -""" -abstract type AbstractQuadratureType{CD, SD} end - -""" -$(TYPEDSIGNATURES) -""" -cell_quadrature_degree(::AbstractQuadratureType{CD, SD}) where {CD, SD} = CD -""" -$(TYPEDSIGNATURES) -""" -surface_quadrature_degree(::AbstractQuadratureType{CD, SD}) where {CD, SD} = SD - -""" -$(TYPEDEF) -""" -struct GaussLobattoLegendre{CD, SD} <: AbstractQuadratureType{CD, SD} - function GaussLobattoLegendre(degree::Int) - return GaussLobattoLegendre{degree, degree}() - end - - function GaussLobattoLegendre(cell_degree::Int, surf_degree::Int) - return GaussLobattoLegendre{cell_degree, surf_degree}() - end - - function GaussLobattoLegendre{CD, SD}() where {CD, SD} - @assert isa(CD, Integer) - @assert isa(SD, Integer) - @assert CD > 0 "Cell quadrature degree must be greater than zero" - @assert SD > 0 "Surface quadrature degree must be greater than zero" - new{CD, SD}() - end -end - -""" -$(TYPEDEF) -""" -struct GaussLegendre{CD, SD} <: AbstractQuadratureType{CD, SD} - function GaussLegendre(degree::Int) - return GaussLegendre{degree, degree}() - end - - function GaussLegendre(cell_degree::Int, surf_degree::Int) - return GaussLegendre{cell_degree, surf_degree}() - end - - function GaussLegendre{CD, SD}() where {CD, SD} - @assert isa(CD, Integer) - @assert isa(SD, Integer) - @assert CD > 0 "Cell quadrature degree must be greater than zero" - @assert SD > 0 "Surface quadrature degree must be greater than zero" - new{CD, SD}() - end -end - -# methods to define for quadrature types -""" -$(TYPEDSIGNATURES) -""" -function cell_quadrature_points_and_weights end -""" -$(TYPEDSIGNATURES) -""" -function surface_quadrature_points_and_weights end - """ $(TYPEDEF) Base type for the package that all element types @@ -138,10 +39,6 @@ function boundary_element end $(TYPEDSIGNATURES) """ function boundary_normals end # always returns 3 x num boundaries -# """ -# $(TYPEDSIGNATURES) -# """ -# function dimension end """ $(TYPEDSIGNATURES) """ @@ -157,11 +54,11 @@ function num_boundaries end """ $(TYPEDSIGNATURES) """ -function num_edges end +function num_edges_per_cell end """ $(TYPEDSIGNATURES) """ -function num_faces end +function num_faces_per_cell end """ $(TYPEDSIGNATURES) """ @@ -200,14 +97,14 @@ boundary_normals(::AbstractVertex) = Matrix{Float64}(undef, 3, 0) edge_vertices(::AbstractVertex) = Matrix{Int}(undef, 0, 0) face_vertices(::AbstractVertex) = Matrix{Int}(undef, 0, 0) num_boundaries(::AbstractVertex) = 0 -num_edges(::AbstractVertex) = 0 -num_faces(::AbstractVertex) = 0 +num_edges_per_cell(::AbstractVertex) = 0 +num_faces_per_cell(::AbstractVertex) = 0 num_vertices_per_cell(::AbstractVertex) = 1 vertex_coordinates(::AbstractVertex) = [0. 0. 0.]' |> collect # this one is an oddity boundary_dofs(::AbstractVertex) = Matrix{Int}(undef, 0, 0) -dof_coordinates(::AbstractVertex) = [0. 0. 0.]' |> collect +# dof_coordinates(::AbstractVertex) = [0. 0. 0.]' |> collect interior_dofs(::AbstractVertex) = Vector{Int}(undef, 0) num_cell_dofs(::AbstractVertex) = 1 num_dofs_on_boundary(::AbstractVertex, ::Int) = 0 @@ -223,8 +120,8 @@ boundary_normals(::AbstractEdge) = [-1. 0. 0.; 1. 0. 0.]' |> collect edge_vertices(::AbstractEdge) = [1 2]' |> collect face_vertices(::AbstractEdge) = Matrix{Int}(undef, 0, 0) num_boundaries(::AbstractEdge) = 2 -num_edges(::AbstractEdge) = 1 -num_faces(::AbstractEdge) = 0 +num_edges_per_cell(::AbstractEdge) = 1 +num_faces_per_cell(::AbstractEdge) = 0 num_vertices_per_cell(::AbstractEdge) = 2 function vertex_coordinates(e::AbstractEdge) # if e.shifted @@ -249,8 +146,8 @@ $(TYPEDEF) """ abstract type AbstractFace{PT, PD} <: AbstractElementType{2, PT, PD} end face_vertices(e::AbstractFace) = 1:num_vertices_per_cell(e) |> collect -num_boundaries(e::AbstractFace) = num_edges(e) -num_faces(::AbstractFace) = 1 +num_boundaries(e::AbstractFace) = num_edges_per_cell(e) +num_faces_per_cell(::AbstractFace) = 1 """ $(TYPEDEF) @@ -266,7 +163,7 @@ edge_vertices(::AbstractQuad) = [ 1 2 3 4; 2 3 4 1 ] -num_edges(::AbstractQuad) = 4 +num_edges_per_cell(::AbstractQuad) = 4 num_vertices_per_cell(::AbstractQuad) = 4 vertex_coordinates(::AbstractQuad) = [ -1. 1. 1. -1.; @@ -288,7 +185,7 @@ edge_vertices(::AbstractTri) = [ 1 2 3; 2 3 1 ] -num_edges(::AbstractTri) = 3 +num_edges_per_cell(::AbstractTri) = 3 num_vertices_per_cell(::AbstractTri) = 3 vertex_coordinates(::AbstractTri) = [ 0. 1. 0.; @@ -302,7 +199,7 @@ vertex_coordinates(::AbstractTri) = [ $(TYPEDEF) """ abstract type AbstractVolume{PT, PD} <: AbstractElementType{3, PT, PD} end -num_boundaries(e::AbstractVolume) = num_faces(e) +num_boundaries(e::AbstractVolume) = num_faces_per_cell(e) """ $(TYPEDEF) @@ -325,8 +222,8 @@ face_vertices(::AbstractHex) = [ 5 6 7 4 2 8 ] num_dofs_on_boundary(::AbstractHex{Lagrange, PD}, ::Int) where PD = PD == 0 ? 4 : (PD + 1) * (PD + 1) -num_edges(::AbstractHex) = 12 -num_faces(::AbstractHex) = 6 +num_edges_per_cell(::AbstractHex) = 12 +num_faces_per_cell(::AbstractHex) = 6 num_vertices_per_cell(::AbstractHex) = 8 vertex_coordinates(::AbstractHex) = [ -1. 1. 1. -1. -1. 1. 1. -1. @@ -360,8 +257,8 @@ face_vertices(::AbstractTet) = [ 4 4 3 2 ] num_dofs_on_boundary(::AbstractTet{Lagrange, PD}, ::Int) where PD = PD == 0 ? 3 : (PD + 1) * (PD + 2) ÷ 2 -num_edges(::AbstractTet) = 6 -num_faces(::AbstractTet) = 4 +num_edges_per_cell(::AbstractTet) = 6 +num_faces_per_cell(::AbstractTet) = 4 num_vertices_per_cell(::AbstractTet) = 4 vertex_coordinates(::AbstractTet) = [ 0. 1. 0. 0. diff --git a/src/PolynomialTypes.jl b/src/PolynomialTypes.jl new file mode 100644 index 0000000..8cabb94 --- /dev/null +++ b/src/PolynomialTypes.jl @@ -0,0 +1,42 @@ +""" +$(TYPEDEF) +""" +abstract type AbstractPolynomialType end +""" +$(TYPEDEF) +""" +struct Hermite <: AbstractPolynomialType +end +""" +$(TYPEDEF) +""" +struct Lagrange <: AbstractPolynomialType +end +""" +$(TYPEDEF) +Planned +""" +struct NedelecFirstKind <: AbstractPolynomialType +end +""" +$(TYPEDEF) +Planned +""" +struct NedelecSecondKind <: AbstractPolynomialType +end +""" +$(TYPEDEF) +special type for ``Vertex`` +""" +struct NoInterpolation <: AbstractPolynomialType +end +""" +$(TYPEDEF) +""" +struct RaviartThomas <: AbstractPolynomialType +end +""" +$(TYPEDEF) +""" +struct Serendipity <: AbstractPolynomialType +end diff --git a/src/QuadratureRules.jl b/src/QuadratureRules.jl new file mode 100644 index 0000000..c54db41 --- /dev/null +++ b/src/QuadratureRules.jl @@ -0,0 +1,568 @@ +""" +$(TYPEDEF) +""" +abstract type AbstractQuadratureType{CD, SD} end +""" +$(TYPEDSIGNATURES) +""" +function cell_quadrature_points_and_weights end +""" +$(TYPEDSIGNATURES) +""" +function surface_quadrature_points_and_weights end +""" +$(TYPEDSIGNATURES) +""" +cell_quadrature_degree(::AbstractQuadratureType{CD, SD}) where {CD, SD} = CD +""" +$(TYPEDSIGNATURES) +""" +surface_quadrature_degree(::AbstractQuadratureType{CD, SD}) where {CD, SD} = SD + +""" +$(TYPEDEF) +""" +struct GaussLegendre{CD, SD} <: AbstractQuadratureType{CD, SD} + function GaussLegendre(degree::Int) + return GaussLegendre{degree, degree}() + end + + function GaussLegendre(cell_degree::Int, surf_degree::Int) + return GaussLegendre{cell_degree, surf_degree}() + end + + function GaussLegendre{CD, SD}() where {CD, SD} + @assert isa(CD, Integer) + @assert isa(SD, Integer) + @assert CD > 0 "Cell quadrature degree must be greater than zero" + @assert SD > 0 "Surface quadrature degree must be greater than zero" + new{CD, SD}() + end +end + +######################################################################## +# Edge implementation +######################################################################## +function cell_quadrature_points_and_weights(e::AbstractEdge, q_rule::GaussLegendre) + ξs, ws = gausslegendre(cell_quadrature_degree(q_rule)) + + # if e.shifted + if _is_shifted(e) + ξs .= (ξs .+ 1.) ./ 2. + ws .= ws ./ 2. + end + + return reshape(ξs, 1, length(ξs)), ws +end + +num_cell_quadrature_points(::AbstractEdge, ::Type{GaussLegendre{CD, SD}}) where {CD, SD} = CD + +function surface_quadrature_points_and_weights(e::AbstractEdge, ::GaussLegendre) + if _is_shifted(e) + x_min = 0. + else + x_min = -1. + end + + ξs = zeros(1, 1, 2) + ξs[1, 1, 1] = x_min + ξs[1, 1, 2] = 1. + ws = ones(1, 2) + return ξs, ws +end + +######################################################################## +# Hex implementation +######################################################################## +function cell_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(boundary_element(e, 0), 0), q_rule) + n = length(ws) + ξ_return = Matrix{eltype(ξs)}(undef, 3, n * n * n) + w_return = Vector{eltype(ξs)}(undef, n * n * n) + for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs, ξs)) + ξ_return[1, q] = ξ[1] + ξ_return[2, q] = ξ[2] + ξ_return[3, q] = ξ[3] + end + for (q, w) in enumerate(Base.Iterators.product(ws, ws, ws)) + w_return[q] = w[1] * w[2] * w[3] + end + return ξ_return, w_return +end + +num_cell_quadrature_points(::AbstractHex, ::Type{GaussLegendre{CD, SD}}) where {CD, SD} = CD * CD * CD + +function surface_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) + + ξ_return = zeros(3, length(ws), 6) + w_return = zeros(length(ws), 6) + + ξ_return[1:2, :, 1] .= ξs + ξ_return[3, :, 1] .= -1. + ξ_return[1, :, 2] .= 1. + ξ_return[2:3, :, 2] .= ξs + ξ_return[1:2, :, 3] .= ξs + ξ_return[3, :, 3] .= 1. + ξ_return[1, :, 4] .= -1. + ξ_return[2:3, :, 4] .= ξs + ξ_return[1, :, 5] .= ξs[1, :] + ξ_return[2, :, 5] .= -1. + ξ_return[3, :, 5] .= ξs[2, :] + ξ_return[1, :, 6] .= ξs[1, :] + ξ_return[2, :, 6] .= 1. + ξ_return[3, :, 6] .= ξs[2, :] + + for n in 1:6 + w_return[:, n] .= ws + end + return ξ_return, w_return +end + +######################################################################## +# Quad implementation +######################################################################## +function cell_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) + ξ_return = Matrix{eltype(ξs)}(undef, 2, length(ws) * length(ws)) + w_return = Vector{eltype(ξs)}(undef, length(ws) * length(ws)) + for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs)) + ξ_return[1, q] = ξ[1] + ξ_return[2, q] = ξ[2] + end + for (q, w) in enumerate(Base.Iterators.product(ws, ws)) + w_return[q] = w[1] * w[2] + end + return ξ_return, w_return +end + +num_cell_quadrature_points(::AbstractQuad, ::Type{GaussLegendre{CD, SD}}) where {CD, SD} = CD * CD + +function surface_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) + + ξ_return = zeros(2, length(ws), 4) + w_return = zeros(length(ws), 4) + + ξ_return[1, :, 1] .= ξs[1, :] + ξ_return[2, :, 1] .= -1. + ξ_return[1, :, 2] .= 1. + ξ_return[2, :, 2] .= ξs[1, :] + ξ_return[1, :, 3] .= ξs[1, :] + ξ_return[2, :, 3] .= 1. + ξ_return[1, :, 4] .= -1. + ξ_return[2, :, 4] .= ξs[1, :] + + for n in 1:4 + w_return[:, n] .= ws + end + return ξ_return, w_return +end + +######################################################################## +# Tet implementation +######################################################################## +function cell_quadrature_points_and_weights(::AbstractTet, q_rule::GaussLegendre) + deg = cell_quadrature_degree(q_rule) + if deg == 1 + # 1-point centroid rule (degree 1) + ξs = Matrix{Float64}(undef, 3, 1) + ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.] + ws = [1. / 6.] + elseif deg == 2 + # 4-point symmetric rule (degree 2) + s = sqrt(5.0) + a = (5. + 3. * s) / 20. + b = (5. - s) / 20. + ξs = Matrix{Float64}(undef, 3, 4) + ξs[:, 1] = [b, b, b] + ξs[:, 2] = [a, b, b] + ξs[:, 3] = [b, a, b] + ξs[:, 4] = [b, b, a] + ws = [1. / 24., 1. / 24., 1. / 24., 1. / 24.] + elseif deg == 3 + # 5-point rule (degree 3) + ξs = Matrix{Float64}(undef, 3, 5) + ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.] + ξs[:, 2] = [1. / 6., 1. / 6., 1. / 6.] + ξs[:, 3] = [1. / 6., 1. / 6., 1. / 2.] + ξs[:, 4] = [1. / 6., 1. / 2., 1. / 6.] + ξs[:, 5] = [1. / 2., 1. / 6., 1. / 6.] + ws = [-2. / 15., 3. / 40., 3. / 40., 3. / 40., 3. / 40.] + else + @assert false "GaussLegendre degree 1 through 3 supported for Tet." + end + return ξs, ws +end + +num_cell_quadrature_points(::AbstractTet, ::Type{GaussLegendre{1, SD}}) where SD = 1 +num_cell_quadrature_points(::AbstractTet, ::Type{GaussLegendre{2, SD}}) where SD = 4 +num_cell_quadrature_points(::AbstractTet, ::Type{GaussLegendre{3, SD}}) where SD = 5 + +function surface_quadrature_points_and_weights(e::AbstractTet, q_rule::GaussLegendre) + return surface_quadrature_points_and_weights(e, GaussLobattoLegendre(cell_quadrature_degree(q_rule), surface_quadrature_degree(q_rule))) +end + +######################################################################## +# Tri implementation +######################################################################## +function cell_quadrature_points_and_weights(::AbstractTri, q_rule::GaussLegendre) + deg = cell_quadrature_degree(q_rule) + if deg == 1 + # 1-point centroid rule (degree 1) + ξs = Matrix{Float64}(undef, 2, 1) + ξs[:, 1] = [1. / 3., 1. / 3.] + ws = [0.5] + elseif deg == 2 + # 3-point rule (degree 2) + ξs = Matrix{Float64}(undef, 2, 3) + ξs[:, 1] = [1. / 6., 1. / 6.] + ξs[:, 2] = [4. / 6., 1. / 6.] + ξs[:, 3] = [1. / 6., 4. / 6.] + ws = [1. / 6., 1. / 6., 1. / 6.] + elseif deg == 3 + # 4-point rule (degree 3): centroid + 3 edge midpoints + ξs = Matrix{Float64}(undef, 2, 4) + ξs[:, 1] = [1. / 3., 1. / 3.] + ξs[:, 2] = [1. / 5., 3. / 5.] + ξs[:, 3] = [3. / 5., 1. / 5.] + ξs[:, 4] = [1. / 5., 1. / 5.] + ws = [-27. / 96., 25. / 96., 25. / 96., 25. / 96.] + else + @assert false "GaussLegendre degree 1 through 3 supported for Tri." + end + return ξs, ws +end + +num_cell_quadrature_points(::AbstractTri, ::Type{GaussLegendre{1, SD}}) where SD = 1 +num_cell_quadrature_points(::AbstractTri, ::Type{GaussLegendre{2, SD}}) where SD = 3 +num_cell_quadrature_points(::AbstractTri, ::Type{GaussLegendre{3, SD}}) where SD = 4 + +function surface_quadrature_points_and_weights(e::AbstractTri, q_rule::GaussLegendre) + return surface_quadrature_points_and_weights(e, GaussLobattoLegendre(cell_quadrature_degree(q_rule), surface_quadrature_degree(q_rule))) +end + +""" +$(TYPEDEF) +""" +struct GaussLobattoLegendre{CD, SD} <: AbstractQuadratureType{CD, SD} + function GaussLobattoLegendre(degree::Int) + return GaussLobattoLegendre{degree, degree}() + end + + function GaussLobattoLegendre(cell_degree::Int, surf_degree::Int) + return GaussLobattoLegendre{cell_degree, surf_degree}() + end + + function GaussLobattoLegendre{CD, SD}() where {CD, SD} + @assert isa(CD, Integer) + @assert isa(SD, Integer) + @assert CD > 0 "Cell quadrature degree must be greater than zero" + @assert SD > 0 "Surface quadrature degree must be greater than zero" + new{CD, SD}() + end +end + +######################################################################## +# Edge implementation +######################################################################## +function cell_quadrature_points_and_weights(e::AbstractEdge, q_rule::GaussLobattoLegendre) + ξs, ws = gausslegendre(cell_quadrature_degree(q_rule)) + + if _is_shifted(e) + ξs .= (ξs .+ 1.) ./ 2. + ws .= ws ./ 2. + end + + return reshape(ξs, 1, length(ξs)), ws +end + +num_cell_quadrature_points(::AbstractEdge, ::Type{GaussLobattoLegendre{CD, SD}}) where {CD, SD} = CD + +function surface_quadrature_points_and_weights(e::AbstractEdge, ::GaussLobattoLegendre) + if _is_shifted(e) + x_min = 0. + else + x_min = -1. + end + + ξs = zeros(1, 1, 2) + ξs[1, 1, 1] = x_min + ξs[1, 1, 2] = 1. + ws = ones(1, 2) + return ξs, ws +end + +######################################################################## +# Hex implementation +######################################################################## +function cell_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLobattoLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(boundary_element(e, 0), 0), q_rule) + ξ_return = Matrix{eltype(ξs)}(undef, 3, length(ξs) * length(ξs) * length(ξs) * length(ξs)) + w_return = Vector{eltype(ξs)}(undef, length(ξs) * length(ξs) * length(ξs) * length(ξs)) + for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs, ξs)) + ξ_return[1, q] = ξ[1] + ξ_return[2, q] = ξ[2] + ξ_return[3, q] = ξ[3] + end + for (q, w) in enumerate(Base.Iterators.product(ws, ws, ws)) + w_return[q] = w[1] * w[2] * w[3] + end + return ξ_return, w_return +end + +num_cell_quadrature_points(::AbstractHex, ::Type{GaussLobattoLegendre{CD, SD}}) where {CD, SD} = CD * CD * CD + +function surface_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLobattoLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) + + ξ_return = zeros(3, length(ws), 6) + w_return = zeros(length(ws), 6) + + ξ_return[1:2, :, 1] .= ξs + ξ_return[3, :, 1] .= -1. + # + ξ_return[1, :, 2] .= 1. + ξ_return[2:3, :, 2] .= ξs + # + ξ_return[1:2, :, 3] .= ξs + ξ_return[3, :, 3] .= 1. + # + ξ_return[1, :, 4] .= -1. + ξ_return[2:3, :, 4] .= ξs + # + ξ_return[1, :, 5] .= ξs[1, :] + ξ_return[2, :, 5] .= -1. + ξ_return[3, :, 5] .= ξs[2, :] + # + ξ_return[1, :, 5] .= ξs[1, :] + ξ_return[2, :, 5] .= 1. + ξ_return[3, :, 5] .= ξs[2, :] + # + # + # ξ_return[1, :, 2] .= 1. + # ξ_return[2, :, 2] .= ξs[1, :] + # ξ_return[1, :, 3] .= ξs[1, :] + # ξ_return[2, :, 3] .= 1. + # ξ_return[1, :, 4] .= -1. + # ξ_return[2, :, 4] .= ξs[1, :] + + for n in 1:6 + w_return[:, n] .= ws + end + return ξ_return, w_return +end + +######################################################################## +# Quad implementation +######################################################################## +function cell_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLobattoLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) + ξ_return = Matrix{eltype(ξs)}(undef, 2, length(ws) * length(ws)) + w_return = Vector{eltype(ξs)}(undef, length(ws) * length(ws)) + for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs)) + ξ_return[1, q] = ξ[1] + ξ_return[2, q] = ξ[2] + end + for (q, w) in enumerate(Base.Iterators.product(ws, ws)) + w_return[q] = w[1] * w[2] + end + return ξ_return, w_return +end + +num_cell_quadrature_points(::AbstractQuad, ::Type{GaussLobattoLegendre{CD, SD}}) where {CD, SD} = CD * CD + +function surface_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLobattoLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) + + ξ_return = zeros(2, length(ws), 4) + w_return = zeros(length(ws), 4) + + ξ_return[1, :, 1] .= ξs[1, :] + ξ_return[2, :, 1] .= -1. + ξ_return[1, :, 2] .= 1. + ξ_return[2, :, 2] .= ξs[1, :] + ξ_return[1, :, 3] .= ξs[1, :] + ξ_return[2, :, 3] .= 1. + ξ_return[1, :, 4] .= -1. + ξ_return[2, :, 4] .= ξs[1, :] + + for n in 1:4 + w_return[:, n] .= ws + end + return ξ_return, w_return +end + +######################################################################## +# Tet implementation +######################################################################## +function cell_quadrature_points_and_weights(::AbstractTet, q_rule::GaussLobattoLegendre) + if cell_quadrature_degree(q_rule) == 1 + ξs = Matrix{Float64}(undef, 3, 1) + ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.] + ws = [1. / 6.] + elseif cell_quadrature_degree(q_rule) == 2 + ξs = Matrix{Float64}(undef, 3, 5) + ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.] + ξs[:, 2] = [1. / 6., 1. / 6., 1. / 6.] + ξs[:, 3] = [1. / 6., 1. / 6., 1. / 2.] + ξs[:, 4] = [1. / 6., 1. / 2., 1. / 6.] + ξs[:, 5] = [1. / 2., 1. / 6., 1. / 6.] + + # + ws = [ + -2. / 15. + 3. / 40. + 3. / 40. + 3. / 40. + 3. / 40. + ] + else + @assert false "Quadrature 1 through 2 currently supported." + end + return ξs, ws +end + +num_cell_quadrature_points(::AbstractTet, ::Type{GaussLobattoLegendre{1, SD}}) where SD = 1 +num_cell_quadrature_points(::AbstractTet, ::Type{GaussLobattoLegendre{2, SD}}) where SD = 5 + +function surface_quadrature_points_and_weights(e::AbstractTet, q_rule::GaussLobattoLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) + + ξ_return = zeros(3, length(ws), 4) + w_return = zeros(length(ws), 4) + + ξ_return[1, :, 1] .= ξs[1, :] + ξ_return[2, :, 1] .= 0. + ξ_return[3, :, 1] .= ξs[2, :] + # + ξ_return[1, :, 2] .= ξs[1, :] + ξ_return[2, :, 2] .= ξs[2, :] + ξ_return[3, :, 2] .= 1. .- ξs[1, :] .- ξs[2, :] + # + ξ_return[1, :, 3] .= 0. + ξ_return[2, :, 3] .= ξs[1, :] + ξ_return[3, :, 3] .= ξs[2, :] + # + ξ_return[1, :, 4] .= ξs[1, :] + ξ_return[2, :, 4] .= ξs[2, :] + ξ_return[3, :, 4] .= 0. + + for n in 1:4 + w_return[:, n] .= ws + end + + return ξ_return, w_return +end + +######################################################################## +# Tri implementation +######################################################################## +function cell_quadrature_points_and_weights(::AbstractTri, q_rule::GaussLobattoLegendre) + if cell_quadrature_degree(q_rule) == 1 + ξs = Matrix{Float64}(undef, 2, 1) + ξs[:, 1] = [1. / 3., 1. / 3.] + ws = [0.5] + elseif cell_quadrature_degree(q_rule) == 2 + ξs = Matrix{Float64}(undef, 2, 3) + ξs[:, 1] = [2. / 3., 1. / 6.] + ξs[:, 2] = [1. / 6., 2. / 3.] + ξs[:, 3] = [1. / 6., 1. / 6.] + ws = [1. / 6., 1. / 6., 1. / 6.] + elseif cell_quadrature_degree(q_rule) <= 4 + ξs = Matrix{Float64}(undef, 2, 6) + ξs[:, 1] = [1.081030181680700E-01, 4.459484909159650E-01] + ξs[:, 2] = [4.459484909159650E-01, 1.081030181680700E-01] + ξs[:, 3] = [4.459484909159650E-01, 4.459484909159650E-01] + ξs[:, 4] = [8.168475729804590E-01, 9.157621350977100E-02] + ξs[:, 5] = [9.157621350977100E-02, 8.168475729804590E-01] + ξs[:, 6] = [9.157621350977100E-02, 9.157621350977100E-02] + + ws = [ + 1.116907948390055E-01, + 1.116907948390055E-01, + 1.116907948390055E-01, + 5.497587182766100E-02, + 5.497587182766100E-02, + 5.497587182766100E-02 + ] + elseif cell_quadrature_degree(q_rule) <= 5 + ξs = Matrix{Float64}(undef, 2, 7) + ξs[:, 1] = [3.33333333333333E-01, 3.33333333333333E-01] + ξs[:, 2] = [5.97158717897700E-02, 4.70142064105115E-01] + ξs[:, 3] = [4.70142064105115E-01, 5.97158717897700E-02] + ξs[:, 4] = [4.70142064105115E-01, 4.70142064105115E-01] + ξs[:, 5] = [7.97426985353087E-01, 1.01286507323456E-01] + ξs[:, 6] = [1.01286507323456E-01, 7.97426985353087E-01] + ξs[:, 7] = [1.01286507323456E-01, 1.01286507323456E-01] + + ws = [ + 1.12500000000000E-01, + 6.61970763942530E-02, + 6.61970763942530E-02, + 6.61970763942530E-02, + 6.29695902724135E-02, + 6.29695902724135E-02, + 6.29695902724135E-02 + ] + elseif cell_quadrature_degree(q_rule) <= 6 + ξs = Matrix{Float64}(undef, 2, 12) + ξs[:, 1] = [5.01426509658179E-01, 2.49286745170910E-01] + ξs[:, 2] = [2.49286745170910E-01, 5.01426509658179E-01] + ξs[:, 3] = [2.49286745170910E-01, 2.49286745170910E-01] + ξs[:, 4] = [8.73821971016996E-01, 6.30890144915020E-02] + ξs[:, 5] = [6.30890144915020E-02, 8.73821971016996E-01] + ξs[:, 6] = [6.30890144915020E-02, 6.30890144915020E-02] + ξs[:, 7] = [5.31450498448170E-02, 3.10352451033784E-01] + ξs[:, 8] = [6.36502499121399E-01, 5.31450498448170E-02] + ξs[:, 9] = [3.10352451033784E-01, 6.36502499121399E-01] + ξs[:, 10] = [5.31450498448170E-02, 6.36502499121399E-01] + ξs[:, 11] = [6.36502499121399E-01, 3.10352451033784E-01] + ξs[:, 12] = [3.10352451033784E-01, 5.31450498448170E-02] + + ws = [ + 5.83931378631895E-02, + 5.83931378631895E-02, + 5.83931378631895E-02, + 2.54224531851035E-02, + 2.54224531851035E-02, + 2.54224531851035E-02, + 4.14255378091870E-02, + 4.14255378091870E-02, + 4.14255378091870E-02, + 4.14255378091870E-02, + 4.14255378091870E-02, + 4.14255378091870E-02 + ] + else + @assert false "Quadrature degree 1 through 6 currently supported." + end + + return ξs, ws +end + +num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{1, SD}}) where SD = 1 +num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{2, SD}}) where SD = 3 +num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{3, SD}}) where SD = 6 +num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{4, SD}}) where SD = 6 +num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{5, SD}}) where SD = 7 +num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{6, SD}}) where SD = 12 + +function surface_quadrature_points_and_weights(e::AbstractTri, q_rule::GaussLobattoLegendre) + ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) + + ξ_return = zeros(2, length(ws), 3) + w_return = zeros(length(ws), 3) + + ξ_return[1, :, 1] .= ξs[1, :] + ξ_return[2, :, 1] .= -1. + ξ_return[1, :, 2] .= ξs[1, :] + ξ_return[2, :, 2] .= 1. .- ξs[1, :] + ξ_return[1, :, 3] .= -1. + ξ_return[2, :, 3] .= ξs[1, :] + + for n in 1:3 + w_return[:, n] .= ws + end + return ξ_return, w_return +end diff --git a/src/ReferenceFE.jl b/src/ReferenceFE.jl index 61e0bd8..52b6b6b 100644 --- a/src/ReferenceFE.jl +++ b/src/ReferenceFE.jl @@ -14,35 +14,6 @@ quadrature_weight(interps::AbstractMatrix{T}, q::Int, f::Int) where T <: Abstrac shape_function_value(interps::AbstractVector{T}, q::Int) where T <: AbstractStaticInterpolants = interps[q].N shape_function_value(interps::AbstractMatrix{T}, q::Int, f::Int) where T <: AbstractStaticInterpolants = interps[q, f].N -# abstract type AbstractH1OrL2Interpolants <: AbstractInterpolants end -# abstract type AbstractHdivInterpolants <: AbstractInterpolants end - -# struct H1OrL2Interpolants{ -# RT <: Number, -# R1 <: Union{<:AbstractArray{RT, 1}, <:AbstractArray{RT, 2}}, -# R2 <: Union{<:AbstractArray{RT, 2}, <:AbstractArray{RT, 3}}, -# R3 <: Union{<:AbstractArray{RT, 3}, <:AbstractArray{RT, 4}}, -# R4 <: Union{<:AbstractArray{RT, 4}, <:AbstractArray{RT, 5}} -# } <: AbstractH1OrL2Interpolants -# ws::R1 -# ξs::R2 -# Ns::R2 -# ∇N_ξs::R3 -# ∇∇N_ξs::R4 -# end - -# num_quadrature_points(interps::H1OrL2Interpolants) = size(interps.ws, 1) -# quadrature_point(interps::H1OrL2Interpolants, q::Int) = view(interps.ξs, :, q) -# quadrature_point(interps::H1OrL2Interpolants, q::Int, f::Int) = view(interps.ξs, :, q, f) -# quadrature_weight(interps::H1OrL2Interpolants, q::Int) = interps.ws[q] -# quadrature_weight(interps::H1OrL2Interpolants, q::Int, f::Int) = interps.ws[q, f] -# shape_function_gradient(interps::H1OrL2Interpolants, q::Int) = view(interps.∇N_ξs, :, :, q) -# shape_function_gradient(interps::H1OrL2Interpolants, q::Int, f::Int) = view(interps.∇N_ξs, :, :, q, f) -# shape_function_hessian(interps::H1OrL2Interpolants, q::Int) = view(interps.∇∇N_ξs, :, :, :, q) -# shape_function_hessian(interps::H1OrL2Interpolants, q::Int, f::Int) = view(interps.∇∇N_ξs, :, :, :, q, f) -# shape_function_value(interps::H1OrL2Interpolants, q::Int) = view(interps.Ns, :, q) -# shape_function_value(interps::H1OrL2Interpolants, q::Int, f::Int) = view(interps.Ns, :, q, f) - abstract type AbstractStaticH1OrL2Interpolants <: AbstractStaticInterpolants end shape_function_gradient(interps::AbstractVector{T}, q::Int) where T <: AbstractStaticH1OrL2Interpolants = interps[q].∇N_ξ shape_function_gradient(interps::AbstractMatrix{T}, q::Int, f::Int) where T <: AbstractStaticH1OrL2Interpolants = interps[q, f].∇N_ξ @@ -57,13 +28,13 @@ struct StaticH1OrL2Interpolants{ ∇N_ξ::SMatrix{NN, ND, RT, NDNN} end -function StaticH1OrL2Interpolants(el_type::AbstractElementType, X, ξ, w) +function StaticH1OrL2Interpolants(el_type::AbstractElementType, ξ, w) ND = dimension(el_type) NN = num_cell_dofs(el_type) return StaticH1OrL2Interpolants( w, SVector{ND, Float64}(ξ), - SVector{NN, Float64}(shape_function_value(el_type, X, ξ)), - SMatrix{NN, ND, Float64, ND * NN}(shape_function_gradient(el_type, X, ξ)) + SVector{NN, Float64}(shape_function_value(el_type, ξ)), + SMatrix{NN, ND, Float64, ND * NN}(shape_function_gradient(el_type, ξ)) ) end @@ -78,14 +49,14 @@ struct StaticH1OrL2InterpolantsWithHessians{ ∇∇N_ξ::SArray{Tuple{NN, ND, ND}, RT, 3, NDNDNN} end -function StaticH1OrL2InterpolantsWithHessians(el_type::AbstractElementType, X, ξ, w) +function StaticH1OrL2InterpolantsWithHessians(el_type::AbstractElementType, ξ, w) ND = dimension(el_type) NN = num_cell_dofs(el_type) return StaticH1OrL2InterpolantsWithHessians( w, SVector{ND, Float64}(ξ), - SVector{NN, Float64}(shape_function_value(el_type, X, ξ)), - SMatrix{NN, ND, Float64, ND * NN}(shape_function_gradient(el_type, X, ξ)), - SArray{Tuple{NN, ND, ND}, Float64, 3, ND * ND * NN}(shape_function_hessian(el_type, X, ξ)) + SVector{NN, Float64}(shape_function_value(el_type, ξ)), + SMatrix{NN, ND, Float64, ND * NN}(shape_function_gradient(el_type, ξ)), + SArray{Tuple{NN, ND, ND}, Float64, 3, ND * ND * NN}(shape_function_hessian(el_type, ξ)) ) end @@ -94,7 +65,6 @@ shape_function_hessian(interps::AbstractMatrix{T}, q::Int, f::Int) where T <: St abstract type AbstractStaticHdivInterpolants <: AbstractStaticInterpolants end - struct StaticHdivInterpolants{ RT <: Number, ND, NN, NDNN @@ -107,84 +77,26 @@ struct StaticHdivInterpolants{ geometry_∇N_ξ::SMatrix{NN, ND, RT, NDNN} end -function StaticHdivInterpolants(el_type::AbstractElementType, X, ξ, w) +function StaticHdivInterpolants(el_type::AbstractElementType, ξ, w) ND = dimension(el_type) NN = num_cell_dofs(el_type) return StaticHdivInterpolants( w, SVector{ND, eltype(ξ)}(ξ), - SMatrix{NN, ND, eltype(ξ), ND * NN}(shape_function_value(el_type, X, ξ)), - SVector{NN, eltype(ξ)}(shape_function_divergence(el_type, X, ξ)), - SVector{NN, eltype(ξ)}(geometry_shape_function_value(el_type, X, ξ)), - SMatrix{NN, ND, eltype(ξ), ND * NN}(geometry_shape_function_gradient(el_type, X, ξ)) + SMatrix{NN, ND, eltype(ξ), ND * NN}(shape_function_value(el_type, ξ)), + SVector{NN, eltype(ξ)}(shape_function_divergence(el_type, ξ)), + SVector{NN, eltype(ξ)}(geometry_shape_function_value(el_type, ξ)), + SMatrix{NN, ND, eltype(ξ), ND * NN}(geometry_shape_function_gradient(el_type, ξ)) ) end shape_function_divergence(interps::AbstractVector{T}, q::Int) where T <: StaticHdivInterpolants = interps[q].divN_ξ shape_function_divergence(interps::AbstractVector{T}, q::Int, f::Int) where T <: StaticHdivInterpolants = interps[q, f].divN_ξ -# function _setup_cell_interpolants( -# el_type::AbstractElementType{Lagrange, PD}, -# q_rule::AbstractQuadratureType, -# ::Type{<:H1OrL2Interpolants} -# ) where PD -# Xs = dof_coordinates(el_type) -# ξs, ws = cell_quadrature_points_and_weights(el_type, q_rule) - -# ND = dimension(el_type) -# NE = num_cell_dofs(el_type) -# NQ = length(ws) - -# Ns = zeros(NE, NQ) -# ∇N_ξs = zeros(ND, NE, NQ) -# ∇∇N_ξs = zeros(ND, ND, NE, NQ) - -# for (q, ξ) in enumerate(eachcol(ξs)) -# if ND == 1 -# Ns[:, q] = shape_function_value(el_type, Xs, ξ[1]) -# else -# Ns[:, q] = shape_function_value(el_type, Xs, ξ) -# end -# ∇N_ξs[:, :, q] = shape_function_gradient(el_type, Xs, ξ) -# ∇∇N_ξs[:, :, :, q] = shape_function_hessian(el_type, Xs, ξ) -# end - -# return H1OrL2Interpolants(ws, ξs, Ns, ∇N_ξs, ∇∇N_ξs) -# end - -# function _setup_cell_interpolants_old( -# el_type::AbstractElementType, -# q_rule::AbstractQuadratureType, -# type::Type{<:AbstractStaticInterpolants} -# ) -# Xs = dof_coordinates(el_type) -# ξs, ws = cell_quadrature_points_and_weights(el_type, q_rule) -# NN = num_cell_dofs(el_type) -# ND = dimension(el_type) -# if type <: StaticH1OrL2Interpolants -# interps = Vector{StaticH1OrL2Interpolants{Float64, ND, NN, ND * NN}}(undef, length(ws)) -# elseif type <: StaticH1OrL2InterpolantsWithHessians -# interps = Vector{StaticH1OrL2InterpolantsWithHessians{Float64, ND, NN, ND * NN, ND * ND * NN}}(undef, length(ws)) -# elseif type <: StaticHdivInterpolants -# interps = Vector{StaticHdivInterpolants{Float64, ND, NN, ND * NN}}(undef, length(ws)) -# else -# @assert false "Unsupported type $type" -# end - -# for (q, (w, ξ)) in enumerate(zip(ws, eachcol(ξs))) -# if typeof(el_type) <: Edge -# ξ = ξ[1] -# end -# interps[q] = type(el_type, Xs, ξ, w) -# end -# return interps -# end - function _setup_cell_interpolants( el_type::AbstractElementType, q_rule::AbstractQuadratureType, ::Type{<:StaticH1OrL2Interpolants} ) - Xs = dof_coordinates(el_type) ξs, ws = cell_quadrature_points_and_weights(el_type, q_rule) NN = num_cell_dofs(el_type) ND = dimension(el_type) @@ -194,7 +106,7 @@ function _setup_cell_interpolants( if typeof(el_type) <: Edge ξ = ξ[1] end - interps[q] = StaticH1OrL2Interpolants(el_type, Xs, ξ, w) + interps[q] = StaticH1OrL2Interpolants(el_type, ξ, w) end return interps end @@ -204,7 +116,6 @@ function _setup_cell_interpolants( q_rule::AbstractQuadratureType, ::Type{<:StaticH1OrL2InterpolantsWithHessians} ) - Xs = dof_coordinates(el_type) ξs, ws = cell_quadrature_points_and_weights(el_type, q_rule) NN = num_cell_dofs(el_type) ND = dimension(el_type) @@ -214,61 +125,36 @@ function _setup_cell_interpolants( if typeof(el_type) <: Edge ξ = ξ[1] end - interps[q] = StaticH1OrL2InterpolantsWithHessians(el_type, Xs, ξ, w) + interps[q] = StaticH1OrL2InterpolantsWithHessians(el_type, ξ, w) end return interps end -# function _setup_surface_interpolants( -# el_type::AbstractElementType{Lagrange, PD}, -# q_rule::AbstractQuadratureType, -# ::Type{<:H1OrL2Interpolants} -# ) where PD -# if dimension(el_type) == 1 -# Xs = zeros(1, 1, 2) -# Xs[1, 1, 1] = -1. -# Xs[1, 1, 2] = 1. -# else -# Xs = dof_coordinates(el_type)[:, boundary_dofs(el_type)] -# end -# ξs, ws = surface_quadrature_points_and_weights(el_type, q_rule) - -# ND = dimension(el_type) -# NE = num_cell_dofs(el_type) -# NQ = size(ξs, 2) -# NF = size(ξs, 3) - -# Ns = zeros(NE, NQ, NF) -# ∇N_ξs = zeros(ND, NE, NQ, NF) -# ∇∇N_ξs = zeros(ND, ND, NE, NQ, NF) - -# if dimension(el_type) == 1 -# Ns[1, :, 1] .= 1. -# Ns[2, :, 2] .= 1. -# else -# for f in axes(ξs, 3) -# for q in axes(ξs, 2) -# Ns[:, q, f] .= shape_function_value(el_type, Xs[:, :, f], ξs[:, q, f]) -# ∇N_ξs[:, :, q, f] .= shape_function_gradient(el_type, Xs[:, :, f], ξs[:, q, f]) -# ∇∇N_ξs[:, :, :, q, f] .= shape_function_hessian(el_type, Xs[:, :, f], ξs[:, q, f]) -# end -# end -# end -# return H1OrL2Interpolants(ws, ξs, Ns, ∇N_ξs, ∇∇N_ξs) -# end +function _setup_cell_interpolants( + el_type::AbstractElementType, + q_rule::AbstractQuadratureType, + ::Type{<:StaticHdivInterpolants} +) + # Xs = dof_coordinates(el_type) + ξs, ws = cell_quadrature_points_and_weights(el_type, q_rule) + NN = num_cell_dofs(el_type) + ND = dimension(el_type) + interps = Vector{StaticHdivInterpolants{Float64, ND, NN, ND * NN}}(undef, length(ws)) + + for (q, (w, ξ)) in enumerate(zip(ws, eachcol(ξs))) + if typeof(el_type) <: Edge + ξ = ξ[1] + end + interps[q] = StaticHdivInterpolants(el_type, ξ, w) + end + return interps +end function _setup_surface_interpolants( el_type::AbstractElementType, q_rule::AbstractQuadratureType, type::Type{<:StaticH1OrL2Interpolants} ) - if dimension(el_type) == 1 - Xs = zeros(1, 1, 2) - Xs[1, 1, 1] = -1. - Xs[1, 1, 2] = 1. - else - Xs = dof_coordinates(el_type)[:, boundary_dofs(el_type)] - end ξs, ws = surface_quadrature_points_and_weights(el_type, q_rule) NN = num_cell_dofs(el_type) @@ -296,7 +182,7 @@ function _setup_surface_interpolants( else ξ = ξs[:, q, f] end - interps[q, f] = type(el_type, Xs[:, :, f], ξ, ws[q, f]) + interps[q, f] = type(el_type, ξ, ws[q, f]) end end end @@ -308,13 +194,6 @@ function _setup_surface_interpolants( q_rule::AbstractQuadratureType, type::Type{<:StaticH1OrL2InterpolantsWithHessians} ) - if dimension(el_type) == 1 - Xs = zeros(1, 1, 2) - Xs[1, 1, 1] = -1. - Xs[1, 1, 2] = 1. - else - Xs = dof_coordinates(el_type)[:, boundary_dofs(el_type)] - end ξs, ws = surface_quadrature_points_and_weights(el_type, q_rule) NN = num_cell_dofs(el_type) @@ -344,14 +223,41 @@ function _setup_surface_interpolants( else ξ = ξs[:, q, f] end - interps[q, f] = type(el_type, Xs[:, :, f], ξ, ws[q, f]) + interps[q, f] = type(el_type, ξ, ws[q, f]) + end + end + end + return interps +end + +function _setup_surface_interpolants( + el_type::AbstractElementType, + q_rule::AbstractQuadratureType, + type::Type{<:StaticHdivInterpolants} +) + Xs = dof_coordinates(el_type)[:, boundary_dofs(el_type)] + ξs, ws = surface_quadrature_points_and_weights(el_type, q_rule) + + NN = num_cell_dofs(el_type) + ND = dimension(el_type) + NB = num_boundaries(el_type) + interps = Matrix{StaticHdivInterpolants{Float64, ND, NN, ND * NN}}(undef, size(ws, 1), NB) + + for f in axes(ξs, 3) + for q in axes(ξs, 2) + if typeof(el_type) <: Edge + ξ = ξs[1, q, f] + else + ξ = ξs[:, q, f] end + interps[q, f] = type(el_type, Xs[:, :, f], ξ, ws[q, f]) end end return interps end _default_interpolants_type(::AbstractElementType{D, Lagrange, PD}) where {D, PD} = StaticH1OrL2Interpolants +_default_interpolants_type(::AbstractElementType{D, RaviartThomas, PD}) where {D, PD} = StaticHdivInterpolants struct ReferenceFE{E, Q, BDofs, BNorms, CellInterps, SurfInterps, NEPE} element::E diff --git a/src/ReferenceFiniteElements.jl b/src/ReferenceFiniteElements.jl index 5dbf4d9..2750a53 100644 --- a/src/ReferenceFiniteElements.jl +++ b/src/ReferenceFiniteElements.jl @@ -9,7 +9,9 @@ using SpecialPolynomials using StaticArrays # includes -include("AbstractTypes.jl") +include("PolynomialTypes.jl") +include("AbstractElements.jl") +include("QuadratureRules.jl") include("ReferenceFE.jl") # 0-d elements @@ -90,8 +92,8 @@ export edge_vertices export element export face_vertices export num_boundaries -export num_edges -export num_faces +export num_edges_per_cell +export num_faces_per_cell export num_vertices_per_cell export vertex_coordinates diff --git a/src/elements/Edge.jl b/src/elements/Edge.jl index f437635..04271e2 100644 --- a/src/elements/Edge.jl +++ b/src/elements/Edge.jl @@ -147,103 +147,48 @@ num_cell_dofs(::Edge{Lagrange, PD}) where PD = PD + 1 num_dofs_on_boundary(::Edge{Lagrange, PD}, ::Int) where PD = 1 num_interior_dofs(::Edge{Lagrange, PD}) where PD = PD < 2 ? 0 : PD - 1 -function cell_quadrature_points_and_weights(e::AbstractEdge, q_rule::GaussLegendre) - ξs, ws = gausslegendre(cell_quadrature_degree(q_rule)) - - # if e.shifted - if _is_shifted(e) - ξs .= (ξs .+ 1.) ./ 2. - ws .= ws ./ 2. - end - - return reshape(ξs, 1, length(ξs)), ws -end - -num_cell_quadrature_points(::AbstractEdge, ::Type{GaussLegendre{CD, SD}}) where {CD, SD} = CD - -function surface_quadrature_points_and_weights(e::AbstractEdge, ::GaussLegendre) - if _is_shifted(e) - x_min = 0. - else - x_min = -1. - end - - ξs = zeros(1, 1, 2) - ξs[1, 1, 1] = x_min - ξs[1, 1, 2] = 1. - ws = ones(1, 2) - return ξs, ws -end - -function cell_quadrature_points_and_weights(e::AbstractEdge, q_rule::GaussLobattoLegendre) - ξs, ws = gausslegendre(cell_quadrature_degree(q_rule)) - - if _is_shifted(e) - ξs .= (ξs .+ 1.) ./ 2. - ws .= ws ./ 2. - end - - return reshape(ξs, 1, length(ξs)), ws -end - -num_cell_quadrature_points(::AbstractEdge, ::Type{GaussLobattoLegendre{CD, SD}}) where {CD, SD} = CD - -function surface_quadrature_points_and_weights(e::AbstractEdge, ::GaussLobattoLegendre) - if _is_shifted(e) - x_min = 0. - else - x_min = -1. - end - - ξs = zeros(1, 1, 2) - ξs[1, 1, 1] = x_min - ξs[1, 1, 2] = 1. - ws = ones(1, 2) - return ξs, ws -end - # 0th order Lagrange implementation -function shape_function_value(::Edge{Lagrange, 0, Shifted}, _, ::Number) where Shifted +function shape_function_value(::Edge{Lagrange, 0, Shifted}, ::Number) where Shifted return ones(1) end -function shape_function_gradient(::Edge{Lagrange, 0, Shifted}, _, ::Number) where Shifted +function shape_function_gradient(::Edge{Lagrange, 0, Shifted}, ::Number) where Shifted return zeros(1, 1) end -function shape_function_hessian(::Edge{Lagrange, 0, Shifted}, _, ::Number) where Shifted +function shape_function_hessian(::Edge{Lagrange, 0, Shifted}, ::Number) where Shifted return zeros(1, 1, 1) end -function shape_function_value(::Edge{Lagrange, 1, false}, _, ξ::Number) +function shape_function_value(::Edge{Lagrange, 1, false}, ξ::Number) return [ 0.5 * (1 - ξ), 0.5 * (1 + ξ) ] end -function shape_function_value(::Edge{Lagrange, 1, true}, _, ξ::Number) +function shape_function_value(::Edge{Lagrange, 1, true}, ξ::Number) return [ 1 - ξ, ξ ] end -function shape_function_gradient(::Edge{Lagrange, 1, false}, _, ξ::Number) +function shape_function_gradient(::Edge{Lagrange, 1, false}, ξ::Number) return [ -0.5, 0.5 ] end -function shape_function_gradient(::Edge{Lagrange, 1, true}, _, ξ::Number) +function shape_function_gradient(::Edge{Lagrange, 1, true}, ξ::Number) return [ -1.0, 1.0 ] end -function shape_function_value(::Edge{Lagrange, 2, false}, _, ξ::Number) +function shape_function_value(::Edge{Lagrange, 2, false}, ξ::Number) return [ 0.5 * ξ * (ξ - 1.0), 0.5 * ξ * (ξ + 1.0), @@ -251,7 +196,7 @@ function shape_function_value(::Edge{Lagrange, 2, false}, _, ξ::Number) ] end -function shape_function_value(::Edge{Lagrange, 2, true}, _, ξ::Number) +function shape_function_value(::Edge{Lagrange, 2, true}, ξ::Number) return [ (1.0 - ξ) * (1.0 - 2.0 * ξ), 4.0 * ξ * (1.0 - ξ), @@ -259,7 +204,7 @@ function shape_function_value(::Edge{Lagrange, 2, true}, _, ξ::Number) ] end -function shape_function_gradient(::Edge{Lagrange, 2, false}, _, ξ::Number) +function shape_function_gradient(::Edge{Lagrange, 2, false}, ξ::Number) return [ 0.5 * (2.0 * ξ - 1.0), 0.5 * (2.0 * ξ + 1.0), @@ -268,7 +213,7 @@ function shape_function_gradient(::Edge{Lagrange, 2, false}, _, ξ::Number) end # Second order shifted -function shape_function_gradient(::Edge{Lagrange, 2, true}, _, ξ::Number) +function shape_function_gradient(::Edge{Lagrange, 2, true}, ξ::Number) return [ 4.0 * ξ - 3.0, 4.0 - 8.0 * ξ, @@ -276,7 +221,8 @@ function shape_function_gradient(::Edge{Lagrange, 2, true}, _, ξ::Number) ] end -function shape_function_value(e::Edge{Lagrange, PD, Shifted}, Xs, ξ::Number) where {PD, Shifted} +function shape_function_value(e::Edge{Lagrange, PD, Shifted}, ξ::Number) where {PD, Shifted} + Xs = dof_coordinates(e) type = _interp_type(e) n_nodes = polynomial_degree(e) + 1 A = zeros(length(Xs), n_nodes) @@ -293,7 +239,8 @@ function shape_function_value(e::Edge{Lagrange, PD, Shifted}, Xs, ξ::Number) wh return N[:, 1] end -function shape_function_gradient(e::Edge{Lagrange, PD, Shifted}, Xs, ξ) where {PD, Shifted} +function shape_function_gradient(e::Edge{Lagrange, PD, Shifted}, ξ) where {PD, Shifted} + Xs = dof_coordinates(e) type = _interp_type(e) n_nodes = polynomial_degree(e) + 1 A = zeros(length(Xs), n_nodes) @@ -310,7 +257,8 @@ function shape_function_gradient(e::Edge{Lagrange, PD, Shifted}, Xs, ξ) where { return reshape(∇N_ξ, n_nodes, 1) end -function shape_function_hessian(e::Edge{Lagrange, PD, Shifted}, Xs, ξ) where {PD, Shifted} +function shape_function_hessian(e::Edge{Lagrange, PD, Shifted}, ξ) where {PD, Shifted} + Xs = dof_coordinates(e) type = _interp_type(e) n_nodes = polynomial_degree(e) + 1 A = zeros(length(Xs), n_nodes) diff --git a/src/elements/Hex.jl b/src/elements/Hex.jl index 63108a8..84f7061 100644 --- a/src/elements/Hex.jl +++ b/src/elements/Hex.jl @@ -4,24 +4,86 @@ $(TYPEDEF) struct Hex{PT, PD} <: AbstractHex{PT, PD} end +function _edge_dof_indices(v1::Int, v2::Int, n::Int, PD::Int) + # Number of vertices + nv = 8 + + # Check valid edge node + @assert 1 ≤ n ≤ PD-1 + + # Edges are ordered as in the standard Hex: + # edge 1: 1→2 + # edge 2: 2→3 + # edge 3: 3→4 + # edge 4: 4→1 + # edge 5: 5→6 + # edge 6: 6→7 + # edge 7: 7→8 + # edge 8: 8→5 + # edge 9: 1→5 + # edge10:2→6 + # edge11:3→7 + # edge12:4→8 + + # List of edges in vertex ordering + edges = [ + (1,2), (2,3), (3,4), (4,1), + (5,6), (6,7), (7,8), (8,5), + (1,5), (2,6), (3,7), (4,8) + ] + + # find edge number + edge_num = findfirst(e -> e == (v1, v2) || e == (v2, v1), edges) + @assert edge_num !== nothing "Edge not found" + + # Each edge has PD-1 interior DOFs + edge_offset = nv + (edge_num-1)*(PD-1) + + # global DOF index + return edge_offset + n +end + +function _face_dof_indices(face::Int, i::Int, j::Int, PD::Int) + # i,j in 2:PD (interior nodes) + # returns (xi_idx, eta_idx, zeta_idx) + if face == 1 # bottom, ζ = 1 + return (i, j, 1) + elseif face == 2 # top, ζ = PD+1 + return (i, j, PD+1) + elseif face == 3 # left, ξ = 1 + return (1, i, j) + elseif face == 4 # right, ξ = PD+1 + return (PD+1, i, j) + elseif face == 5 # front, η = 1 + return (i, 1, j) + elseif face == 6 # back, η = PD+1 + return (i, PD+1, j) + else + error("Invalid face index $face") + end +end + +######################################################################## +# Lagrange implementation +######################################################################## function boundary_dofs(e::Hex{Lagrange, PD}) where PD # 6 faces nv = num_vertices_per_cell(e) # 8 - ne = num_edges(e) # 12 + ne = num_edges_per_cell(e) # 12 faces = face_vertices(e) # 4 × 6 edges = edge_vertices(e) # 2 × 12 # Prepare container # max DOFs per face: 4 vertices + 4 * (PD - 1) edge DOFs + (PD - 1)^2 face interior max_rows = 4 + 4 * (PD - 1) + (PD - 1)^2 - F = zeros(Int, max_rows, num_faces(e)) + F = zeros(Int, max_rows, num_faces_per_cell(e)) # Offset counters edge_offset = nv + 1 # first edge DOF index face_offset = nv + ne * (PD - 1) + 1 # first interior face DOF index if PD > 1 - for f = 1:num_faces(e) + for f = 1:num_faces_per_cell(e) # corner vertices F[1:4, f] .= faces[:, f] @@ -64,67 +126,7 @@ function boundary_dofs(e::Hex{Lagrange, PD}) where PD return faces end end -function dof_coordinates(e::Hex{Lagrange, PD}) where PD - if PD == 0 - return zeros(3, 1) - end - - Xv = vertex_coordinates(e) # 3 × 8 - coords = copy(Xv) # start with vertices - - if PD > 1 - # -------------------------- - # edge DOFs - # -------------------------- - edge_pts_1d = range(-1.0, 1.0, PD + 1)[2:end - 1] # interior points along [-1,1] - - for (v1, v2) in eachcol(edge_vertices(e)) - for ξ in edge_pts_1d - pt = (1.0 - ξ) / 2 * Xv[:, v1] + - (1.0 + ξ) / 2 * Xv[:, v2] - coords = hcat(coords, pt) - end - end - - # -------------------------- - # face interior DOFs - # -------------------------- - face_pts_1d = edge_pts_1d # same 1D points along face - for f in 1:num_faces(e) - v = face_vertices(e)[:, f] - v1, v2, v3, v4 = v - for i in face_pts_1d, j in face_pts_1d - # bilinear interpolation in the face - pt = (1.0 - i) * (1.0 - j) / 4 * Xv[:, v1] + - (1.0 + i) * (1.0 - j) / 4 * Xv[:, v2] + - (1.0 + i) * (1.0 + j) / 4 * Xv[:, v3] + - (1.0 - i) * (1.0 + j) / 4 * Xv[:, v4] - coords = hcat(coords, pt) - end - end - - # -------------------------- - # cell interior DOFs - # -------------------------- - pts_1d = edge_pts_1d - for i in pts_1d, j in pts_1d, k in pts_1d - # trilinear interpolation inside the hex - pt = zeros(eltype(coords), 3) - for (α, β, γ, vtx) in ( - (-1, -1, -1, 1), (1, -1, -1, 2), - (1, 1, -1, 3), (-1, 1, -1, 4), - (-1, -1, 1, 5), (1, -1, 1, 6), - (1, 1, 1, 7), (-1, 1, 1, 8) - ) - w = (1 + α * i) / 2 * (1 + β * j) / 2 * (1 + γ * k) / 2 - pt += w * Xv[:, vtx] - end - coords = hcat(coords, pt) - end - end - return coords -end function interior_dofs(e::Hex{Lagrange, PD}) where PD if PD < 2 return Int[] @@ -133,7 +135,9 @@ function interior_dofs(e::Hex{Lagrange, PD}) where PD return offset:offset + num_interior_dofs(e) end end + num_cell_dofs(::Hex{Lagrange, PD}) where PD = (PD + 1)^3 + function num_interior_dofs(::Hex{Lagrange, PD}) where PD if PD == 0 return 1 @@ -144,180 +148,19 @@ function num_interior_dofs(::Hex{Lagrange, PD}) where PD end end -function cell_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLegendre) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(boundary_element(e, 0), 0), q_rule) - n = length(ws) - ξ_return = Matrix{eltype(ξs)}(undef, 3, n * n * n) - w_return = Vector{eltype(ξs)}(undef, n * n * n) - for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs, ξs)) - ξ_return[1, q] = ξ[1] - ξ_return[2, q] = ξ[2] - ξ_return[3, q] = ξ[3] - end - for (q, w) in enumerate(Base.Iterators.product(ws, ws, ws)) - w_return[q] = w[1] * w[2] * w[3] - end - return ξ_return, w_return -end - -num_cell_quadrature_points(::AbstractHex, ::Type{GaussLegendre{CD, SD}}) where {CD, SD} = CD * CD * CD - -function surface_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLegendre) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) - - ξ_return = zeros(3, length(ws), 6) - w_return = zeros(length(ws), 6) - - ξ_return[1:2, :, 1] .= ξs - ξ_return[3, :, 1] .= -1. - ξ_return[1, :, 2] .= 1. - ξ_return[2:3, :, 2] .= ξs - ξ_return[1:2, :, 3] .= ξs - ξ_return[3, :, 3] .= 1. - ξ_return[1, :, 4] .= -1. - ξ_return[2:3, :, 4] .= ξs - ξ_return[1, :, 5] .= ξs[1, :] - ξ_return[2, :, 5] .= -1. - ξ_return[3, :, 5] .= ξs[2, :] - ξ_return[1, :, 6] .= ξs[1, :] - ξ_return[2, :, 6] .= 1. - ξ_return[3, :, 6] .= ξs[2, :] - - for n in 1:6 - w_return[:, n] .= ws - end - return ξ_return, w_return -end - -function cell_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLobattoLegendre) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(boundary_element(e, 0), 0), q_rule) - ξ_return = Matrix{eltype(ξs)}(undef, 3, length(ξs) * length(ξs) * length(ξs) * length(ξs)) - w_return = Vector{eltype(ξs)}(undef, length(ξs) * length(ξs) * length(ξs) * length(ξs)) - for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs, ξs)) - ξ_return[1, q] = ξ[1] - ξ_return[2, q] = ξ[2] - ξ_return[3, q] = ξ[3] - end - for (q, w) in enumerate(Base.Iterators.product(ws, ws, ws)) - w_return[q] = w[1] * w[2] * w[3] - end - return ξ_return, w_return -end - -num_cell_quadrature_points(::AbstractHex, ::Type{GaussLobattoLegendre{CD, SD}}) where {CD, SD} = CD * CD * CD - -function surface_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLobattoLegendre) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) - - ξ_return = zeros(3, length(ws), 6) - w_return = zeros(length(ws), 6) - - ξ_return[1:2, :, 1] .= ξs - ξ_return[3, :, 1] .= -1. - # - ξ_return[1, :, 2] .= 1. - ξ_return[2:3, :, 2] .= ξs - # - ξ_return[1:2, :, 3] .= ξs - ξ_return[3, :, 3] .= 1. - # - ξ_return[1, :, 4] .= -1. - ξ_return[2:3, :, 4] .= ξs - # - ξ_return[1, :, 5] .= ξs[1, :] - ξ_return[2, :, 5] .= -1. - ξ_return[3, :, 5] .= ξs[2, :] - # - ξ_return[1, :, 5] .= ξs[1, :] - ξ_return[2, :, 5] .= 1. - ξ_return[3, :, 5] .= ξs[2, :] - # - # - # ξ_return[1, :, 2] .= 1. - # ξ_return[2, :, 2] .= ξs[1, :] - # ξ_return[1, :, 3] .= ξs[1, :] - # ξ_return[2, :, 3] .= 1. - # ξ_return[1, :, 4] .= -1. - # ξ_return[2, :, 4] .= ξs[1, :] - - for n in 1:6 - w_return[:, n] .= ws - end - return ξ_return, w_return -end - -function _edge_dof_indices(v1::Int, v2::Int, n::Int, PD::Int) - # Number of vertices - nv = 8 - - # Check valid edge node - @assert 1 ≤ n ≤ PD-1 - - # Edges are ordered as in the standard Hex: - # edge 1: 1→2 - # edge 2: 2→3 - # edge 3: 3→4 - # edge 4: 4→1 - # edge 5: 5→6 - # edge 6: 6→7 - # edge 7: 7→8 - # edge 8: 8→5 - # edge 9: 1→5 - # edge10:2→6 - # edge11:3→7 - # edge12:4→8 - - # List of edges in vertex ordering - edges = [ - (1,2), (2,3), (3,4), (4,1), - (5,6), (6,7), (7,8), (8,5), - (1,5), (2,6), (3,7), (4,8) - ] - - # find edge number - edge_num = findfirst(e -> e == (v1, v2) || e == (v2, v1), edges) - @assert edge_num !== nothing "Edge not found" - - # Each edge has PD-1 interior DOFs - edge_offset = nv + (edge_num-1)*(PD-1) - - # global DOF index - return edge_offset + n -end - -function _face_dof_indices(face::Int, i::Int, j::Int, PD::Int) - # i,j in 2:PD (interior nodes) - # returns (xi_idx, eta_idx, zeta_idx) - if face == 1 # bottom, ζ = 1 - return (i, j, 1) - elseif face == 2 # top, ζ = PD+1 - return (i, j, PD+1) - elseif face == 3 # left, ξ = 1 - return (1, i, j) - elseif face == 4 # right, ξ = PD+1 - return (PD+1, i, j) - elseif face == 5 # front, η = 1 - return (i, 1, j) - elseif face == 6 # back, η = PD+1 - return (i, PD+1, j) - else - error("Invalid face index $face") - end -end - -function shape_function_value(::Hex{Lagrange, 0}, _, _) +function shape_function_value(::Hex{Lagrange, 0}, _) return ones(1) end -function shape_function_gradient(::Hex{Lagrange, 0}, _, _) +function shape_function_gradient(::Hex{Lagrange, 0}, _) return zeros(1, 3) end -function shape_function_hessian(::Hex{Lagrange, 0}, _, _) +function shape_function_hessian(::Hex{Lagrange, 0}, _) return zeros(1, 3, 3) end -function shape_function_value(::Hex{Lagrange, 1}, _, ξ) +function shape_function_value(::Hex{Lagrange, 1}, ξ) Ns = [ 0.125 * (1 - ξ[1]) * (1 - ξ[2]) * (1 - ξ[3]), 0.125 * (1 + ξ[1]) * (1 - ξ[2]) * (1 - ξ[3]), @@ -331,36 +174,9 @@ function shape_function_value(::Hex{Lagrange, 1}, _, ξ) return Ns end -function shape_function_gradient(::Hex{Lagrange, 1}, _, ξ) +function shape_function_gradient(::Hex{Lagrange, 1}, ξ) Ns = zeros(8, 3) - # Ns = reshape([ - # -0.125 * (1 - ξ[2]) * (1 - ξ[3]), - # 0.125 * (1 - ξ[2]) * (1 - ξ[3]), - # 0.125 * (1 + ξ[2]) * (1 - ξ[3]), - # -0.125 * (1 + ξ[2]) * (1 - ξ[3]), - # -0.125 * (1 - ξ[2]) * (1 + ξ[3]), - # 0.125 * (1 - ξ[2]) * (1 + ξ[3]), - # 0.125 * (1 + ξ[2]) * (1 + ξ[3]), - # -0.125 * (1 + ξ[2]) * (1 + ξ[3]), - # # - # -0.125 * (1 - ξ[1]) * (1 - ξ[3]), - # -0.125 * (1 + ξ[1]) * (1 - ξ[3]), - # 0.125 * (1 + ξ[1]) * (1 - ξ[3]), - # 0.125 * (1 - ξ[1]) * (1 - ξ[3]), - # -0.125 * (1 - ξ[1]) * (1 + ξ[3]), - # -0.125 * (1 + ξ[1]) * (1 + ξ[3]), - # 0.125 * (1 + ξ[1]) * (1 + ξ[3]), - # 0.125 * (1 - ξ[1]) * (1 + ξ[3]), - # # - # -0.125 * (1 - ξ[1]) * (1 - ξ[2]), - # -0.125 * (1 + ξ[1]) * (1 - ξ[2]), - # -0.125 * (1 + ξ[1]) * (1 + ξ[2]), - # -0.125 * (1 - ξ[1]) * (1 + ξ[2]), - # 0.125 * (1 - ξ[1]) * (1 - ξ[2]), - # 0.125 * (1 + ξ[1]) * (1 - ξ[2]), - # 0.125 * (1 + ξ[1]) * (1 + ξ[2]), - # 0.125 * (1 - ξ[1]) * (1 + ξ[2]) - # ], 8, 3)' |> collect + # Ns[1, 1] = -0.125 * (1 - ξ[2]) * (1 - ξ[3]) Ns[2, 1] = 0.125 * (1 - ξ[2]) * (1 - ξ[3]) Ns[3, 1] = 0.125 * (1 + ξ[2]) * (1 - ξ[3]) @@ -390,11 +206,11 @@ function shape_function_gradient(::Hex{Lagrange, 1}, _, ξ) return Ns end -function shape_function_hessian(::Hex{Lagrange, 1}, _, _) +function shape_function_hessian(::Hex{Lagrange, 1}, _) return zeros(3, 3, 8) end -function shape_function_value(e::Hex{Lagrange, PD}, _, ξ) where PD +function shape_function_value(e::Hex{Lagrange, PD}, ξ) where PD # ξ = [ξ, η, ζ] in reference coordinates le = boundary_element(boundary_element(e)) @@ -494,7 +310,7 @@ function shape_function_value(e::Hex{Lagrange, PD}, _, ξ) where PD return N end -function shape_function_gradient(e::Hex{Lagrange, PD}, _, ξ) where PD +function shape_function_gradient(e::Hex{Lagrange, PD}, ξ) where PD le = boundary_element(boundary_element(e)) coords_1d = dof_coordinates(le) @@ -645,3 +461,8 @@ end # return H # end + +######################################################################## +# Raviart-Thomas implementation +######################################################################## +# TODO \ No newline at end of file diff --git a/src/elements/Quad.jl b/src/elements/Quad.jl index bfd9af0..006db2d 100644 --- a/src/elements/Quad.jl +++ b/src/elements/Quad.jl @@ -19,43 +19,7 @@ function boundary_dofs(e::Quad{Lagrange, PD}) where PD return edges end end -function dof_coordinates(e::Quad{Lagrange, PD}) where PD - # if PD == 0 - # return zeros(2, 1) - # end - coords = vertex_coordinates(e)[1:2, :] - - if PD > 1 - # do edge midpoints - edge_coords = dof_coordinates(boundary_element(e, 0)) - - # face 1 - for n in 1:PD - 1 - coords = hcat(coords, [edge_coords[1, n + 2], -1.]) - end - # face 2 - for n in 1:PD - 1 - coords = hcat(coords, [1., edge_coords[1, n + 2]]) - end - # face 3 - for n in 1:PD - 1 - coords = hcat(coords, [edge_coords[1, n + 2], 1.]) - end - # face 4 - for n in 1:PD - 1 - coords = hcat(coords, [-1., edge_coords[1, n + 2]]) - end - - # now for interiors - for n in 1:PD - 1 - for m in 1:PD - 1 - coords = hcat(coords, [edge_coords[1, m + 2], edge_coords[1, n + 2]]) - end - end - end - return coords -end function interior_dofs(e::Quad{Lagrange, PD}) where PD if PD < 2 return Int[] @@ -64,9 +28,11 @@ function interior_dofs(e::Quad{Lagrange, PD}) where PD return offset:offset + num_interior_dofs(e) - 1 |> collect end end + num_cell_dofs(::Quad{Lagrange, PD}) where PD = (PD + 1) * (PD + 1) # provides bdofs...? num_dofs_on_boundary(::Quad{Lagrange, PD}, ::Int) where PD = PD == 0 ? 2 : PD + 1 + function num_interior_dofs(::Quad{Lagrange, PD}) where PD if PD == 0 return 1 @@ -77,98 +43,21 @@ function num_interior_dofs(::Quad{Lagrange, PD}) where PD end end -function cell_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLegendre) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) - ξ_return = Matrix{eltype(ξs)}(undef, 2, length(ws) * length(ws)) - w_return = Vector{eltype(ξs)}(undef, length(ws) * length(ws)) - for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs)) - ξ_return[1, q] = ξ[1] - ξ_return[2, q] = ξ[2] - end - for (q, w) in enumerate(Base.Iterators.product(ws, ws)) - w_return[q] = w[1] * w[2] - end - return ξ_return, w_return -end - -num_cell_quadrature_points(::AbstractQuad, ::Type{GaussLegendre{CD, SD}}) where {CD, SD} = CD * CD - -function surface_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLegendre) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) - - ξ_return = zeros(2, length(ws), 4) - w_return = zeros(length(ws), 4) - - ξ_return[1, :, 1] .= ξs[1, :] - ξ_return[2, :, 1] .= -1. - ξ_return[1, :, 2] .= 1. - ξ_return[2, :, 2] .= ξs[1, :] - ξ_return[1, :, 3] .= ξs[1, :] - ξ_return[2, :, 3] .= 1. - ξ_return[1, :, 4] .= -1. - ξ_return[2, :, 4] .= ξs[1, :] - - for n in 1:4 - w_return[:, n] .= ws - end - return ξ_return, w_return -end - -function cell_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLobattoLegendre) - # ξs, ws = gausslegendre(cell_quadrature_degree(e)) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) - ξ_return = Matrix{eltype(ξs)}(undef, 2, length(ws) * length(ws)) - w_return = Vector{eltype(ξs)}(undef, length(ws) * length(ws)) - for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs)) - ξ_return[1, q] = ξ[1] - ξ_return[2, q] = ξ[2] - end - for (q, w) in enumerate(Base.Iterators.product(ws, ws)) - w_return[q] = w[1] * w[2] - end - return ξ_return, w_return -end - -num_cell_quadrature_points(::AbstractQuad, ::Type{GaussLobattoLegendre{CD, SD}}) where {CD, SD} = CD * CD - -function surface_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLobattoLegendre) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) - - ξ_return = zeros(2, length(ws), 4) - w_return = zeros(length(ws), 4) - - ξ_return[1, :, 1] .= ξs[1, :] - ξ_return[2, :, 1] .= -1. - ξ_return[1, :, 2] .= 1. - ξ_return[2, :, 2] .= ξs[1, :] - ξ_return[1, :, 3] .= ξs[1, :] - ξ_return[2, :, 3] .= 1. - ξ_return[1, :, 4] .= -1. - ξ_return[2, :, 4] .= ξs[1, :] - - for n in 1:4 - w_return[:, n] .= ws - end - return ξ_return, w_return -end - -function shape_function_value(::Quad{Lagrange, 0}, _, _) +function shape_function_value(::Quad{Lagrange, 0}, _) return ones(1) end -function shape_function_gradient(::Quad{Lagrange, 0}, _, _) +function shape_function_gradient(::Quad{Lagrange, 0}, _) return zeros(1, 2) end -function shape_function_hessian(::Quad{Lagrange, 0}, _, _) +function shape_function_hessian(::Quad{Lagrange, 0}, _) return zeros(1, 2, 2) end -function shape_function_value(e::Quad{Lagrange, PD}, _, ξ) where PD - coords_x = dof_coordinates(boundary_element(e, 0)) - coords_y = dof_coordinates(boundary_element(e, 0)) - N_x = shape_function_value(boundary_element(e, 0), coords_x, ξ[1]) - N_y = shape_function_value(boundary_element(e, 0), coords_y, ξ[2]) +function shape_function_value(e::Quad{Lagrange, PD}, ξ) where PD + N_x = shape_function_value(boundary_element(e, 0), ξ[1]) + N_y = shape_function_value(boundary_element(e, 0), ξ[2]) N = Vector{eltype(ξ)}(undef, num_cell_dofs(e)) @@ -211,13 +100,11 @@ function shape_function_value(e::Quad{Lagrange, PD}, _, ξ) where PD return N end -function shape_function_gradient(e::Quad{Lagrange, PD}, X, ξ) where PD - coords_x = dof_coordinates(boundary_element(e, 0)) - coords_y = dof_coordinates(boundary_element(e, 0)) - N_x = shape_function_value(boundary_element(e, 0), coords_x, ξ[1]) - N_y = shape_function_value(boundary_element(e, 0), coords_y, ξ[2]) - ∇N_x = shape_function_gradient(boundary_element(e, 0), coords_x, ξ[1]) - ∇N_y = shape_function_gradient(boundary_element(e, 0), coords_y, ξ[2]) +function shape_function_gradient(e::Quad{Lagrange, PD}, ξ) where PD + N_x = shape_function_value(boundary_element(e, 0), ξ[1]) + N_y = shape_function_value(boundary_element(e, 0), ξ[2]) + ∇N_x = shape_function_gradient(boundary_element(e, 0), ξ[1]) + ∇N_y = shape_function_gradient(boundary_element(e, 0), ξ[2]) # return N_x * N_y @@ -275,15 +162,13 @@ function shape_function_gradient(e::Quad{Lagrange, PD}, X, ξ) where PD return ∇N end -function shape_function_hessian(e::Quad{Lagrange, PD}, X, ξ) where PD - coords_x = dof_coordinates(boundary_element(e, 0)) - coords_y = dof_coordinates(boundary_element(e, 0)) - N_x = shape_function_value(boundary_element(e, 0), coords_x, ξ[1]) - N_y = shape_function_value(boundary_element(e, 0), coords_y, ξ[2]) - ∇N_x = shape_function_gradient(boundary_element(e, 0), coords_x, ξ[1]) - ∇N_y = shape_function_gradient(boundary_element(e, 0), coords_y, ξ[2]) - ∇∇N_x = shape_function_hessian(boundary_element(e, 0), coords_x, ξ[1]) - ∇∇N_y = shape_function_hessian(boundary_element(e, 0), coords_y, ξ[2]) +function shape_function_hessian(e::Quad{Lagrange, PD}, ξ) where PD + N_x = shape_function_value(boundary_element(e, 0), ξ[1]) + N_y = shape_function_value(boundary_element(e, 0), ξ[2]) + ∇N_x = shape_function_gradient(boundary_element(e, 0), ξ[1]) + ∇N_y = shape_function_gradient(boundary_element(e, 0), ξ[2]) + ∇∇N_x = shape_function_hessian(boundary_element(e, 0), ξ[1]) + ∇∇N_y = shape_function_hessian(boundary_element(e, 0), ξ[2]) ∇∇N = Array{eltype(ξ), 3}(undef, num_cell_dofs(e), 2, 2) @@ -362,3 +247,50 @@ function shape_function_hessian(e::Quad{Lagrange, PD}, X, ξ) where PD return ∇∇N end + +######################################################################## +# Raviart-Thomas implementation +######################################################################## +function boundary_dofs(::Quad{RaviartThomas, 0}) + return reshape(collect(1:4), 1, 4) +end +interior_dofs(::Quad{RaviartThomas, 0}) = Int[] +num_cell_dofs(::Quad{RaviartThomas, 0}) = 4 +num_dofs_on_boundary(::Quad{RaviartThomas, 0}, ::Int) = 1 +num_interior_dofs(::Quad{RaviartThomas, 0}) = 0 + +function geometry_shape_function_value(::Quad{RaviartThomas, 0}, ξ) + return shape_function_value(Quad{Lagrange, 1}(), X, ξ) +end + +function geometry_shape_function_gradient(::Quad{RaviartThomas, 0}, ξ) + return shape_function_gradient(Quad{Lagrange, 1}(), ξ) +end + +# https://defelement.org/elements/examples/quadrilateral-raviart-thomas-lagrange-0.html +# but re-ordered for exodus +function shape_function_value(::Quad{RaviartThomas, 0}, ξ) + N = Matrix{Float64}(undef, 4, 2) + + # bottom + N[1, 1] = 0.0 + N[1, 2] = (1.0 - ξ[2]) / 2.0 + + # right + N[2, 1] = -(ξ[1] + 1.0) / 2.0 + N[2, 2] = 0.0 + + # top + N[3, 1] = 0.0 + N[3, 2] = (ξ[2] + 1.0) / 2.0 + + # left + N[4, 1] = (ξ[1] - 1.0) / 2.0 + N[4, 2] = 0.0 + + return N +end + +function shape_function_divergence(::Quad{RaviartThomas, 0}, ξ) + return [-0.5, -0.5, 0.5, 0.5] +end diff --git a/src/elements/Tet.jl b/src/elements/Tet.jl index 5fd8cc1..38edbd7 100644 --- a/src/elements/Tet.jl +++ b/src/elements/Tet.jl @@ -84,84 +84,6 @@ function boundary_dofs(e::Tet{Lagrange, PD}) where PD return faces end -function dof_coordinates(e::Tet{Lagrange, PD}) where PD - coords = vertex_coordinates(e) - - offset = 4 - - if PD < 2 - # return coords - # do nothing to coords - elseif PD == 2 - X = zeros(3, 10) - - # Vertices - X[:, 1] .= (0, 0, 0) # node 1 - X[:, 2] .= (1, 0, 0) # node 2 - X[:, 3] .= (0, 1, 0) # node 3 - X[:, 4] .= (0, 0, 1) # node 4 - - # Mid-edge nodes (ExodusII ordering) - X[:, 5] .= (1//2, 0, 0 ) # edge 1-2 - X[:, 6] .= (1//2, 1//2, 0 ) # edge 2-3 - X[:, 7] .= (0, 1//2, 0 ) # edge 3-1 - X[:, 8] .= (0, 0, 1//2 ) # edge 1-4 - X[:, 9] .= (1//2, 0, 1//2 ) # edge 2-4 - X[:, 10] .= (0, 1//2, 1//2 ) # edge 3-4 - - # return X - coords = X - else - - # edge DOFs - if PD ≥ 2 - for (v1, v2) in eachcol(edge_vertices(e)) - for i in 1:PD - 1 - t = i / PD - new_coord = (1.0 - t) * coords[:, v1] + t * coords[:, v2] - coords = hcat(coords, new_coord) - offset += 1 - end - end - end - - # face DOFs - if PD ≥ 3 - for (v1, v2, v3) in face_vertices(e) - for i in 1:PD - 1 - for j in 1:PD - 1 - i - t1 = i / PD - t2 = j / PD - t3 = 1.0 - t1 - t2 - new_coord = t1 * coords[:,v1] + t2 * coords[:,v2] + t3 * coords[:,v3] - coords = hcat(coords, new_coord) - offset += 1 - end - end - end - end - - # interior DOFs - if PD ≥ 4 - for i in 1:PD - 3 - for j in 1:PD - 2 - i - for k in 1:PD - 1 - i - j - t1 = i / PD - t2 = j / PD - t3 = k / PD - t4 = 1.0 - t1 - t2 - t3 - new_coord = t1 * coords[:,1] + t2 * coords[:,2] + t3 * coords[:,3] + t4 * coords[:,4] - coords = hcat(coords, new_coord) - offset += 1 - end - end - end - end - end - - return coords -end - function interior_dofs(::Tet{Lagrange, PD}) where PD if PD < 4 return Int[] @@ -178,120 +100,19 @@ end num_cell_dofs(::Tet{Lagrange, PD}) where PD = (PD + 1) * (PD + 2) * (PD + 3) ÷ 6 num_interior_dofs(::Tet{Lagrange, PD}) where PD = PD < 4 ? 0 : (PD - 1) * (PD - 2) * (PD - 3) ÷ 6 -function cell_quadrature_points_and_weights(::AbstractTet, q_rule::GaussLegendre) - deg = cell_quadrature_degree(q_rule) - if deg == 1 - # 1-point centroid rule (degree 1) - ξs = Matrix{Float64}(undef, 3, 1) - ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.] - ws = [1. / 6.] - elseif deg == 2 - # 4-point symmetric rule (degree 2) - s = sqrt(5.0) - a = (5. + 3. * s) / 20. - b = (5. - s) / 20. - ξs = Matrix{Float64}(undef, 3, 4) - ξs[:, 1] = [b, b, b] - ξs[:, 2] = [a, b, b] - ξs[:, 3] = [b, a, b] - ξs[:, 4] = [b, b, a] - ws = [1. / 24., 1. / 24., 1. / 24., 1. / 24.] - elseif deg == 3 - # 5-point rule (degree 3) - ξs = Matrix{Float64}(undef, 3, 5) - ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.] - ξs[:, 2] = [1. / 6., 1. / 6., 1. / 6.] - ξs[:, 3] = [1. / 6., 1. / 6., 1. / 2.] - ξs[:, 4] = [1. / 6., 1. / 2., 1. / 6.] - ξs[:, 5] = [1. / 2., 1. / 6., 1. / 6.] - ws = [-2. / 15., 3. / 40., 3. / 40., 3. / 40., 3. / 40.] - else - @assert false "GaussLegendre degree 1 through 3 supported for Tet." - end - return ξs, ws -end - -num_cell_quadrature_points(::AbstractTet, ::Type{GaussLegendre{1, SD}}) where SD = 1 -num_cell_quadrature_points(::AbstractTet, ::Type{GaussLegendre{2, SD}}) where SD = 4 -num_cell_quadrature_points(::AbstractTet, ::Type{GaussLegendre{3, SD}}) where SD = 5 - -function surface_quadrature_points_and_weights(e::AbstractTet, q_rule::GaussLegendre) - return surface_quadrature_points_and_weights(e, GaussLobattoLegendre(cell_quadrature_degree(q_rule), surface_quadrature_degree(q_rule))) -end - -function cell_quadrature_points_and_weights(e::AbstractTet, q_rule::GaussLobattoLegendre) - if cell_quadrature_degree(q_rule) == 1 - ξs = Matrix{Float64}(undef, 3, 1) - ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.] - ws = [1. / 6.] - elseif cell_quadrature_degree(q_rule) == 2 - ξs = Matrix{Float64}(undef, 3, 5) - ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.] - ξs[:, 2] = [1. / 6., 1. / 6., 1. / 6.] - ξs[:, 3] = [1. / 6., 1. / 6., 1. / 2.] - ξs[:, 4] = [1. / 6., 1. / 2., 1. / 6.] - ξs[:, 5] = [1. / 2., 1. / 6., 1. / 6.] - - # - ws = [ - -2. / 15. - 3. / 40. - 3. / 40. - 3. / 40. - 3. / 40. - ] - else - @assert false "Quadrature 1 through 2 currently supported." - end - return ξs, ws -end - -num_cell_quadrature_points(::AbstractTet, ::Type{GaussLobattoLegendre{1, SD}}) where SD = 1 -num_cell_quadrature_points(::AbstractTet, ::Type{GaussLobattoLegendre{2, SD}}) where SD = 5 - - -function surface_quadrature_points_and_weights(e::AbstractTet, q_rule::GaussLobattoLegendre) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) - - ξ_return = zeros(3, length(ws), 4) - w_return = zeros(length(ws), 4) - - ξ_return[1, :, 1] .= ξs[1, :] - ξ_return[2, :, 1] .= 0. - ξ_return[3, :, 1] .= ξs[2, :] - # - ξ_return[1, :, 2] .= ξs[1, :] - ξ_return[2, :, 2] .= ξs[2, :] - ξ_return[3, :, 2] .= 1. .- ξs[1, :] .- ξs[2, :] - # - ξ_return[1, :, 3] .= 0. - ξ_return[2, :, 3] .= ξs[1, :] - ξ_return[3, :, 3] .= ξs[2, :] - # - ξ_return[1, :, 4] .= ξs[1, :] - ξ_return[2, :, 4] .= ξs[2, :] - ξ_return[3, :, 4] .= 0. - - for n in 1:4 - w_return[:, n] .= ws - end - - return ξ_return, w_return -end - -function shape_function_value(::Tet{Lagrange, 0}, _, _) +function shape_function_value(::Tet{Lagrange, 0}, _) return ones(1) end -function shape_function_gradient(::Tet{Lagrange, 0}, _, _) +function shape_function_gradient(::Tet{Lagrange, 0}, _) return zeros(1, 3) end -function shape_function_hessian(::Tet{Lagrange, 0}, _, _) +function shape_function_hessian(::Tet{Lagrange, 0}, _) return zeros(1, 3, 3) end -function shape_function_value(::Tet{Lagrange, 1}, _, ξ) +function shape_function_value(::Tet{Lagrange, 1}, ξ) return [ 1. - ξ[1] - ξ[2] - ξ[3], ξ[1], @@ -300,7 +121,7 @@ function shape_function_value(::Tet{Lagrange, 1}, _, ξ) ] end -function shape_function_gradient(::Tet{Lagrange, 1}, _, ξ) +function shape_function_gradient(::Tet{Lagrange, 1}, ξ) grads = zeros(4, 3) grads[1, 1] = -1. @@ -315,11 +136,11 @@ function shape_function_gradient(::Tet{Lagrange, 1}, _, ξ) return grads end -function shape_function_hessian(::Tet{Lagrange, 1}, _, ξ) +function shape_function_hessian(::Tet{Lagrange, 1}, ξ) return zeros(4, 3, 3) end -function shape_function_value(::Tet{Lagrange, 2}, _, ξ) +function shape_function_value(::Tet{Lagrange, 2}, ξ) t0 = 1 - ξ[1] - ξ[2] - ξ[3] t1 = ξ[1] t2 = ξ[2] @@ -339,7 +160,7 @@ function shape_function_value(::Tet{Lagrange, 2}, _, ξ) return Ns end -function shape_function_gradient(::Tet{Lagrange, 2}, _, ξ) +function shape_function_gradient(::Tet{Lagrange, 2}, ξ) t0 = 1 - ξ[1] - ξ[2] - ξ[3] t1 = ξ[1] t2 = ξ[2] @@ -382,7 +203,7 @@ function shape_function_gradient(::Tet{Lagrange, 2}, _, ξ) return grads end -function shape_function_hessian(::Tet{Lagrange, 2}, _, _) +function shape_function_hessian(::Tet{Lagrange, 2}, _) hess = zeros(10, 3, 3) hess[:, 1, 1] .= [4., 4., 0., 0., -8., 0., 0., 0., 0., 0.] hess[:, 1, 2] .= [4., 0., 0., 0., -4., 4., -4., 0., 0., 0.] @@ -396,7 +217,7 @@ function shape_function_hessian(::Tet{Lagrange, 2}, _, _) return hess end -function shape_function_value(e::Tet{Lagrange, PD}, _, ξ) where PD +function shape_function_value(e::Tet{Lagrange, PD}, ξ) where PD # barycentric coordinates λ1 = 1 - ξ[1] - ξ[2] - ξ[3] λ2 = ξ[1] @@ -476,7 +297,7 @@ function shape_function_value(e::Tet{Lagrange, PD}, _, ξ) where PD return N end -function shape_function_gradient(e::Tet{Lagrange, PD}, _, ξ) where PD +function shape_function_gradient(e::Tet{Lagrange, PD}, ξ) where PD λ1 = 1 - ξ[1] - ξ[2] - ξ[3] λ2 = ξ[1] λ3 = ξ[2] diff --git a/src/elements/Tri.jl b/src/elements/Tri.jl index 9df7e84..142ddc4 100644 --- a/src/elements/Tri.jl +++ b/src/elements/Tri.jl @@ -1,148 +1,3 @@ -function cell_quadrature_points_and_weights(::AbstractTri, q_rule::GaussLegendre) - deg = cell_quadrature_degree(q_rule) - if deg == 1 - # 1-point centroid rule (degree 1) - ξs = Matrix{Float64}(undef, 2, 1) - ξs[:, 1] = [1. / 3., 1. / 3.] - ws = [0.5] - elseif deg == 2 - # 3-point rule (degree 2) - ξs = Matrix{Float64}(undef, 2, 3) - ξs[:, 1] = [1. / 6., 1. / 6.] - ξs[:, 2] = [4. / 6., 1. / 6.] - ξs[:, 3] = [1. / 6., 4. / 6.] - ws = [1. / 6., 1. / 6., 1. / 6.] - elseif deg == 3 - # 4-point rule (degree 3): centroid + 3 edge midpoints - ξs = Matrix{Float64}(undef, 2, 4) - ξs[:, 1] = [1. / 3., 1. / 3.] - ξs[:, 2] = [1. / 5., 3. / 5.] - ξs[:, 3] = [3. / 5., 1. / 5.] - ξs[:, 4] = [1. / 5., 1. / 5.] - ws = [-27. / 96., 25. / 96., 25. / 96., 25. / 96.] - else - @assert false "GaussLegendre degree 1 through 3 supported for Tri." - end - return ξs, ws -end - -num_cell_quadrature_points(::AbstractTri, ::Type{GaussLegendre{1, SD}}) where SD = 1 -num_cell_quadrature_points(::AbstractTri, ::Type{GaussLegendre{2, SD}}) where SD = 3 -num_cell_quadrature_points(::AbstractTri, ::Type{GaussLegendre{3, SD}}) where SD = 4 - -function surface_quadrature_points_and_weights(e::AbstractTri, q_rule::GaussLegendre) - return surface_quadrature_points_and_weights(e, GaussLobattoLegendre(cell_quadrature_degree(q_rule), surface_quadrature_degree(q_rule))) -end - -function cell_quadrature_points_and_weights(::AbstractTri, q_rule::GaussLobattoLegendre) - if cell_quadrature_degree(q_rule) == 1 - ξs = Matrix{Float64}(undef, 2, 1) - ξs[:, 1] = [1. / 3., 1. / 3.] - ws = [0.5] - elseif cell_quadrature_degree(q_rule) == 2 - ξs = Matrix{Float64}(undef, 2, 3) - ξs[:, 1] = [2. / 3., 1. / 6.] - ξs[:, 2] = [1. / 6., 2. / 3.] - ξs[:, 3] = [1. / 6., 1. / 6.] - ws = [1. / 6., 1. / 6., 1. / 6.] - elseif cell_quadrature_degree(q_rule) <= 4 - ξs = Matrix{Float64}(undef, 2, 6) - ξs[:, 1] = [1.081030181680700E-01, 4.459484909159650E-01] - ξs[:, 2] = [4.459484909159650E-01, 1.081030181680700E-01] - ξs[:, 3] = [4.459484909159650E-01, 4.459484909159650E-01] - ξs[:, 4] = [8.168475729804590E-01, 9.157621350977100E-02] - ξs[:, 5] = [9.157621350977100E-02, 8.168475729804590E-01] - ξs[:, 6] = [9.157621350977100E-02, 9.157621350977100E-02] - - ws = [ - 1.116907948390055E-01, - 1.116907948390055E-01, - 1.116907948390055E-01, - 5.497587182766100E-02, - 5.497587182766100E-02, - 5.497587182766100E-02 - ] - elseif cell_quadrature_degree(q_rule) <= 5 - ξs = Matrix{Float64}(undef, 2, 7) - ξs[:, 1] = [3.33333333333333E-01, 3.33333333333333E-01] - ξs[:, 2] = [5.97158717897700E-02, 4.70142064105115E-01] - ξs[:, 3] = [4.70142064105115E-01, 5.97158717897700E-02] - ξs[:, 4] = [4.70142064105115E-01, 4.70142064105115E-01] - ξs[:, 5] = [7.97426985353087E-01, 1.01286507323456E-01] - ξs[:, 6] = [1.01286507323456E-01, 7.97426985353087E-01] - ξs[:, 7] = [1.01286507323456E-01, 1.01286507323456E-01] - - ws = [ - 1.12500000000000E-01, - 6.61970763942530E-02, - 6.61970763942530E-02, - 6.61970763942530E-02, - 6.29695902724135E-02, - 6.29695902724135E-02, - 6.29695902724135E-02 - ] - elseif cell_quadrature_degree(q_rule) <= 6 - ξs = Matrix{Float64}(undef, 2, 12) - ξs[:, 1] = [5.01426509658179E-01, 2.49286745170910E-01] - ξs[:, 2] = [2.49286745170910E-01, 5.01426509658179E-01] - ξs[:, 3] = [2.49286745170910E-01, 2.49286745170910E-01] - ξs[:, 4] = [8.73821971016996E-01, 6.30890144915020E-02] - ξs[:, 5] = [6.30890144915020E-02, 8.73821971016996E-01] - ξs[:, 6] = [6.30890144915020E-02, 6.30890144915020E-02] - ξs[:, 7] = [5.31450498448170E-02, 3.10352451033784E-01] - ξs[:, 8] = [6.36502499121399E-01, 5.31450498448170E-02] - ξs[:, 9] = [3.10352451033784E-01, 6.36502499121399E-01] - ξs[:, 10] = [5.31450498448170E-02, 6.36502499121399E-01] - ξs[:, 11] = [6.36502499121399E-01, 3.10352451033784E-01] - ξs[:, 12] = [3.10352451033784E-01, 5.31450498448170E-02] - - ws = [ - 5.83931378631895E-02, - 5.83931378631895E-02, - 5.83931378631895E-02, - 2.54224531851035E-02, - 2.54224531851035E-02, - 2.54224531851035E-02, - 4.14255378091870E-02, - 4.14255378091870E-02, - 4.14255378091870E-02, - 4.14255378091870E-02, - 4.14255378091870E-02, - 4.14255378091870E-02 - ] - else - @assert false "Quadrature degree 1 through 6 currently supported." - end - - return ξs, ws -end - -num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{1, SD}}) where SD = 1 -num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{2, SD}}) where SD = 3 -num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{3, SD}}) where SD = 6 -num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{4, SD}}) where SD = 6 -num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{5, SD}}) where SD = 7 -num_cell_quadrature_points(::AbstractTri, ::Type{GaussLobattoLegendre{6, SD}}) where SD = 12 - -function surface_quadrature_points_and_weights(e::AbstractTri, q_rule::GaussLobattoLegendre) - ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule) - - ξ_return = zeros(2, length(ws), 3) - w_return = zeros(length(ws), 3) - - ξ_return[1, :, 1] .= ξs[1, :] - ξ_return[2, :, 1] .= -1. - ξ_return[1, :, 2] .= ξs[1, :] - ξ_return[2, :, 2] .= 1. .- ξs[1, :] - ξ_return[1, :, 3] .= -1. - ξ_return[2, :, 3] .= ξs[1, :] - - for n in 1:3 - w_return[:, n] .= ws - end - return ξ_return, w_return -end - """ $(TYPEDEF) """ @@ -167,44 +22,7 @@ function boundary_dofs(e::Tri{Lagrange, PD}) where PD return edges end end -function dof_coordinates(e::Tri{Lagrange, PD}) where PD - if PD == 0 - return zeros(2, 1) - end - - coords = [ - 0. 1. 0.; - 0. 0. 1. - ] - - if PD > 1 - # do edge midpoints - edge_coords = dof_coordinates(boundary_element(e, 0)) - - # face 1 - for n in 1:PD - 1 - coords = hcat(coords, [edge_coords[1, n + 2], -1.]) - end - # face 2 - for n in 1:PD - 1 - coords = hcat(coords, [edge_coords[1, n + 2], 1. - edge_coords[1, n + 2]]) - end - # face 3 - for n in 1:PD - 1 - coords = hcat(coords, [-1., edge_coords[1, n + 2]]) - end - # now for interiors - # TODO not correct yet - for n in 1:PD - 1 - for m in 1:PD - 1 - n - # @assert false "TODO" - coords = hcat(coords, [edge_coords[1, m + 2], edge_coords[1, n + 2]]) - end - end - end - return coords -end function interior_dofs(::Tri{Lagrange, PD}) where PD if PD < 3 return Int[] @@ -212,23 +30,24 @@ function interior_dofs(::Tri{Lagrange, PD}) where PD @assert false "TODO" end end + num_cell_dofs(::Tri{Lagrange, PD}) where PD = (PD + 1) * (PD + 2) ÷ 2 num_dofs_on_boundary(::Tri{Lagrange, PD}, ::Int) where PD = PD == 0 ? 2 : PD + 1 num_interior_dofs(::Tri{Lagrange, PD}) where PD = PD < 3 ? 0 : (PD - 1) * (PD - 2) ÷ 2 -function shape_function_value(::Tri{Lagrange, 0}, _, _) +function shape_function_value(::Tri{Lagrange, 0}, _) return ones(1) end -function shape_function_gradient(::Tri{Lagrange, 0}, _, _) +function shape_function_gradient(::Tri{Lagrange, 0}, _) return zeros(1, 2) end -function shape_function_hessian(::Tri{Lagrange, 0}, _, _) +function shape_function_hessian(::Tri{Lagrange, 0}, _) return zeros(1, 2, 2) end -function shape_function_value(::Tri{Lagrange, PD}, _, ξ) where PD +function shape_function_value(::Tri{Lagrange, PD}, ξ) where PD λ1 = 1 - ξ[1] - ξ[2] λ2 = ξ[1] λ3 = ξ[2] @@ -282,7 +101,7 @@ function shape_function_value(::Tri{Lagrange, PD}, _, ξ) where PD return N end -function shape_function_gradient(e::Tri{Lagrange, PD}, _, ξ) where PD +function shape_function_gradient(e::Tri{Lagrange, PD}, ξ) where PD λ1 = 1 - ξ[1] - ξ[2] λ2 = ξ[1] λ3 = ξ[2] @@ -363,7 +182,7 @@ function shape_function_gradient(e::Tri{Lagrange, PD}, _, ξ) where PD return dN end -function shape_function_hessian(e::Tri{Lagrange, PD}, _, ξ) where PD +function shape_function_hessian(e::Tri{Lagrange, PD}, ξ) where PD T = eltype(ξ) λ1 = one(T) - ξ[1] - ξ[2] @@ -445,29 +264,31 @@ interior_dofs(::Tri{RaviartThomas, 0}) = Int[] num_cell_dofs(::Tri{RaviartThomas, 0}) = 3 num_interior_dofs(::Tri{RaviartThomas, 0}) = 0 -function geometry_shape_function_value(::Tri{RaviartThomas, 0}, X, ξ) - return shape_function_value(Tri{Lagrange, 1}(), X, ξ) +function geometry_shape_function_value(::Tri{RaviartThomas, 0}, ξ) + return shape_function_value(Tri{Lagrange, 1}(), ξ) end -function geometry_shape_function_gradient(::Tri{RaviartThomas, 0}, X, ξ) - return shape_function_gradient(Tri{Lagrange, 1}(), X, ξ) +function geometry_shape_function_gradient(::Tri{RaviartThomas, 0}, ξ) + return shape_function_gradient(Tri{Lagrange, 1}(), ξ) end -function shape_function_value(::Tri{RaviartThomas, 0}, _, ξ) +# https://defelement.org/elements/examples/triangle-raviart-thomas-lagrange-0.html +# but re-ordered for exodus numbering +function shape_function_value(::Tri{RaviartThomas, 0}, ξ) N = Matrix{Float64}(undef, 3, 2) - - N[1, 1] = 1 - N[1, 2] = -ξ[2] + # + N[1, 1] = -ξ[1] + N[1, 2] = 1 - ξ[2] # N[2, 1] = -ξ[1] - N[2, 2] = 1 - ξ[2] + N[2, 2] = -ξ[2] # - N[3, 1] = ξ[1] + N[3, 1] = ξ[1] - 1 N[3, 2] = ξ[2] return N end -function shape_function_divergence(::Tri{RaviartThomas, 0}, _, ξ) +function shape_function_divergence(::Tri{RaviartThomas, 0}, ξ) return [-2., -2., 2.] end diff --git a/src/elements/Vertex.jl b/src/elements/Vertex.jl index 79b28e4..18abbe3 100644 --- a/src/elements/Vertex.jl +++ b/src/elements/Vertex.jl @@ -4,14 +4,14 @@ $(TYPEDEF) struct Vertex <: AbstractVertex end -function shape_function_value(::Vertex, _, _) +function shape_function_value(::Vertex, _) return [1.] end -function shape_function_gradient(::Vertex, _, _) +function shape_function_gradient(::Vertex, _) return Matrix{Float64}(undef, 0, 1) end -function shape_function_hessian(::Vertex, _, _) +function shape_function_hessian(::Vertex, _) return Array{Float64, 3}(undef, 0, 0, 1) end diff --git a/test/runtests.jl b/test/runtests.jl index 596dabb..9f82b55 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -53,7 +53,7 @@ function test_q_points_inside_element(re::ReferenceFE) # TODO re-enable # @test all(is_inside_element.((re.element,), surface_quadrature_points(re))) - # for f in 1:num_faces(re.element) + # for f in 1:num_faces_per_cell(re.element) # for n in 1:num_quadrature_points(surface_element(re.element)) # if typeof(re.element) <: ReferenceFiniteElements.AbstractTri # @test is_inside_element(ReferenceFiniteElements.surface_element(re.element), 2. * surface_quadrature_point(re, n, f) .- 1.) @@ -174,8 +174,8 @@ function test_topology_interface_vertex() @test edge_vertices(re) == Matrix{Int}(undef, 0, 0) @test face_vertices(re) == Matrix{Int}(undef, 0, 0) @test num_boundaries(re) == 0 - @test num_edges(re) == 0 - @test num_faces(re) == 0 + @test num_edges_per_cell(re) == 0 + @test num_faces_per_cell(re) == 0 @test num_vertices_per_cell(re) == 1 # @test num_vertices_per_edge(re) == 0 # @test num_vertices_per_face(re) == 0 @@ -192,8 +192,8 @@ function test_topology_interface_edge(interp_type, p, shifted) @test edge_vertices(re) == [1 2]' |> collect @test face_vertices(re) == Matrix{Int}(undef, 0, 0) @test num_boundaries(re) == 2 - @test num_edges(re) == 1 - @test num_faces(re) == 0 + @test num_edges_per_cell(re) == 1 + @test num_faces_per_cell(re) == 0 @test num_vertices_per_cell(re) == 2 # @test num_vertices_per_edge(re) == 2 # @test num_vertices_per_face(re) == 0 @@ -224,8 +224,8 @@ function test_topology_interface_quad(interp_type, p) ] @test face_vertices(re) == 1:4 |> collect @test num_boundaries(re) == 4 - @test num_edges(re) == 4 - @test num_faces(re) == 1 + @test num_edges_per_cell(re) == 4 + @test num_faces_per_cell(re) == 1 @test num_vertices_per_cell(re) == 4 # @test num_vertices_per_edge(re) == 2 # @test num_vertices_per_face(re) == 4 @@ -254,8 +254,8 @@ function test_topology_interface_tri(interp_type, p) ] @test face_vertices(re) == 1:3 |> collect @test num_boundaries(re) == 3 - @test num_edges(re) == 3 - @test num_faces(re) == 1 + @test num_edges_per_cell(re) == 3 + @test num_faces_per_cell(re) == 1 @test num_vertices_per_cell(re) == 3 # @test num_vertices_per_edge(re) == 2 # @test num_vertices_per_face(re) == 3 @@ -289,8 +289,8 @@ function test_topology_interface_hex(interp_type, p) 5 6 7 4 2 8 ] @test num_boundaries(re) == 6 - @test num_edges(re) == 12 - @test num_faces(re) == 6 + @test num_edges_per_cell(re) == 12 + @test num_faces_per_cell(re) == 6 @test num_vertices_per_cell(re) == 8 # @test num_vertices_per_edge(re) == 2 # @test num_vertices_per_face(re) == 4 @@ -323,8 +323,8 @@ function test_topology_interface_tet(interp_type, p) 4 4 3 2 ] @test num_boundaries(re) == 4 - @test num_edges(re) == 6 - @test num_faces(re) == 4 + @test num_edges_per_cell(re) == 6 + @test num_faces_per_cell(re) == 4 @test num_vertices_per_cell(re) == 4 # @test num_vertices_per_edge(re) == 2 # @test num_vertices_per_face(re) == 3 @@ -338,7 +338,7 @@ end function test_dof_interface_vertex() re = Vertex() @test boundary_dofs(re) == Matrix{Int}(undef, 0, 0) - @test dof_coordinates(re) ≈ [0. 0. 0.]' |> collect + # @test dof_coordinates(re) ≈ [0. 0. 0.]' |> collect @test num_cell_dofs(re) == 1 @test num_dofs_on_boundary(re, 0) == 0 for n in 1:num_boundaries(re) @@ -397,52 +397,52 @@ function test_dof_interface_quad(interp_type::Type{Lagrange}, p) offset += p - 1 end end - coords = dof_coordinates(re) - if p == 0 - @test coords ≈ [ - -1. 1. 1. -1. - -1. -1. 1. 1. - ] - elseif p == 1 - @test coords ≈ [ - -1. 1. 1. -1. - -1. -1. 1. 1. - ] - else - edge_coords = dof_coordinates(boundary_element(re, 0)) - # test faces - offset = 5 - # face 1 - for n in 1:p - 1 - @test coords[:, offset + n - 1] ≈ [edge_coords[1, n + 2], -1.] - end - offset += p - 1 - # face 2 - for n in 1:p - 1 - @test coords[:, offset + n - 1] ≈ [1., edge_coords[1, n + 2]] - end - offset += p - 1 - # face 3 - for n in 1:p - 1 - @test coords[:, offset + n - 1] ≈ [edge_coords[1, n + 2], 1.] - end - offset += p - 1 - # face 4 - for n in 1:p - 1 - @test coords[:, offset + n - 1] ≈ [-1., edge_coords[1, n + 2]] - end - offset += p - 1 - - # test interiors - offset = 4 + 4 * (p - 1) + 1 - carry = 1 - # for n in 1:p - 1 - # for m in 1:p - 1 - # @test coords[:, offset + carry - 1] ≈ [edge_coords[1, m + 2], edge_coords[1, n + 2]] - # carry += 1 - # end - # end - end + # coords = dof_coordinates(re) + # if p == 0 + # @test coords ≈ [ + # -1. 1. 1. -1. + # -1. -1. 1. 1. + # ] + # elseif p == 1 + # @test coords ≈ [ + # -1. 1. 1. -1. + # -1. -1. 1. 1. + # ] + # else + # edge_coords = dof_coordinates(boundary_element(re, 0)) + # # test faces + # offset = 5 + # # face 1 + # for n in 1:p - 1 + # @test coords[:, offset + n - 1] ≈ [edge_coords[1, n + 2], -1.] + # end + # offset += p - 1 + # # face 2 + # for n in 1:p - 1 + # @test coords[:, offset + n - 1] ≈ [1., edge_coords[1, n + 2]] + # end + # offset += p - 1 + # # face 3 + # for n in 1:p - 1 + # @test coords[:, offset + n - 1] ≈ [edge_coords[1, n + 2], 1.] + # end + # offset += p - 1 + # # face 4 + # for n in 1:p - 1 + # @test coords[:, offset + n - 1] ≈ [-1., edge_coords[1, n + 2]] + # end + # offset += p - 1 + + # # test interiors + # offset = 4 + 4 * (p - 1) + 1 + # carry = 1 + # # for n in 1:p - 1 + # # for m in 1:p - 1 + # # @test coords[:, offset + carry - 1] ≈ [edge_coords[1, m + 2], edge_coords[1, n + 2]] + # # carry += 1 + # # end + # # end + # end if p < 2 @test interior_dofs(re) == Int[] else @@ -455,45 +455,45 @@ end # TODO finish this up by testing boundary_dofs and interior_dofs function test_dof_interface_tri(interp_type::Type{Lagrange}, p) re = Tri{interp_type, p}() - coords = dof_coordinates(re) - if p == 0 - @test coords ≈ zeros(2, 1) - elseif p == 1 - @test coords ≈ [ - 0. 1. 0.; - 0. 0. 1. - ] - else - edge_coords = dof_coordinates(boundary_element(re, 0)) - # test faces - offset = 4 - # face 1 - for n in 1:p - 1 - @test coords[:, offset + n - 1] ≈ [edge_coords[1, n + 2], -1.] - end - offset += p - 1 - # face 2 - for n in 1:p - 1 - @test coords[:, offset + n - 1] ≈ [edge_coords[1, n + 2], 1. - edge_coords[1, n + 2]] - end - offset += p - 1 - # face 3 - for n in 1:p - 1 - @test coords[:, offset + n - 1] ≈ [-1., edge_coords[1, n + 2]] - end - offset += p - 1 - - # TODO fix this up - # test interiors - offset = 3 + 3 * (p - 1) + 1 - carry = 1 - for n in 1:p - 1 - for m in 1:p - 1 - n - @test coords[:, offset + carry - 1] ≈ [edge_coords[1, m + 2], edge_coords[1, n + 2]] - carry += 1 - end - end - end + # coords = dof_coordinates(re) + # if p == 0 + # @test coords ≈ zeros(2, 1) + # elseif p == 1 + # @test coords ≈ [ + # 0. 1. 0.; + # 0. 0. 1. + # ] + # else + # edge_coords = dof_coordinates(boundary_element(re, 0)) + # # test faces + # offset = 4 + # # face 1 + # for n in 1:p - 1 + # @test coords[:, offset + n - 1] ≈ [edge_coords[1, n + 2], -1.] + # end + # offset += p - 1 + # # face 2 + # for n in 1:p - 1 + # @test coords[:, offset + n - 1] ≈ [edge_coords[1, n + 2], 1. - edge_coords[1, n + 2]] + # end + # offset += p - 1 + # # face 3 + # for n in 1:p - 1 + # @test coords[:, offset + n - 1] ≈ [-1., edge_coords[1, n + 2]] + # end + # offset += p - 1 + + # # TODO fix this up + # # test interiors + # offset = 3 + 3 * (p - 1) + 1 + # carry = 1 + # for n in 1:p - 1 + # for m in 1:p - 1 - n + # @test coords[:, offset + carry - 1] ≈ [edge_coords[1, m + 2], edge_coords[1, n + 2]] + # carry += 1 + # end + # end + # end @test num_cell_dofs(re) == (p + 1) * (p + 2) ÷ 2 if p < 3 @test num_interior_dofs(re) == 0