Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implementation of Nedelec interpolation on tetrahedral and hexahedral elements #1162

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Implementation of Nedelec interpolation on tetrahedral and hexahedral elements #1162

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
321ff6d
Add 1st order Nedelec interpolation on RefTetrahedron
gijswl Feb 24, 2025
67483df
Update Nedelec{RefTetrahedron,1} interpolation based on continuity un…
gijswl Feb 25, 2025
a20d735
Add 1st order Nedelec interpolation on RefHexahedron
gijswl Feb 25, 2025
d6648bf
Update reference_shape_value based on code review feedback
gijswl Feb 26, 2025
8f47022
1st order Raviart-Thomas interpolations on tetrahedrons & hexahedrons
gijswl Mar 1, 2025
286c6ec
Fix mistake in test_interpolations.jl
gijswl Mar 1, 2025
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
114 changes: 114 additions & 0 deletions src/interpolations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1839,6 +1839,64 @@ function get_direction(::RaviartThomas{RefTriangle, 2}, j, cell)
return ifelse(edge[2] > edge[1], 1, -1)
end

# RefTetrahedron, 1st order Lagrange
# https://defelement.com/elements/examples/tetrahedron-raviart-thomas-lagrange-1.html
function reference_shape_value(ip::RaviartThomas{RefTetrahedron, 1}, ξ::Vec{3}, i::Int)
x, y, z = ξ
i == 1 && return Vec(2 * x, 2 * y, 2 * (z - 1)) # Changed sign, follow Ferrite's sign convention
i == 2 && return Vec(2 * x, 2 * (y - 1), 2 * z)
i == 3 && return Vec(2 * x, 2 * y, 2 * z)
i == 4 && return Vec(2 * (x - 1), 2 * y, 2 * z) # Changed sign, follow Ferrite's sign convention
throw(ArgumentError("no shape function $i for interpolation $ip"))
end

getnbasefunctions(::RaviartThomas{RefTetrahedron, 1}) = 4
edgedof_indices(ip::RaviartThomas{RefTetrahedron, 1}) = edgedof_interior_indices(ip)
facedof_indices(ip::RaviartThomas{RefTetrahedron, 1}) = facedof_interior_indices(ip)
edgedof_interior_indices(::RaviartThomas{RefTetrahedron, 1}) = ((), (), (), (), (), ())
facedof_interior_indices(::RaviartThomas{RefTetrahedron, 1}) = ((1), (2), (3), (4))
adjust_dofs_during_distribution(::RaviartThomas{RefTetrahedron, 1}) = false

function get_direction(::RaviartThomas{RefTetrahedron, 1}, j, cell)
face = faces(cell)[j]
_, idx_min = findmin(face)

idx_next = mod(idx_min, length(face)) + 1
idx_prev = mod(idx_min - 2, length(face)) + 1
return face[idx_next] < face[idx_prev] ? 1 : -1
end

# RefHexahedron, 1st order Lagrange
# https://defelement.com/elements/examples/hexahedron-raviart-thomas-lagrange-1.html
function reference_shape_value(ip::RaviartThomas{RefHexahedron, 1}, ξ::Vec{3,T}, i::Int) where {T}
x, y, z = ξ
nil = zero(T)

i == 1 && return Vec(nil, nil, (z - 1) / 2) # Changed sign, follow Ferrite's sign convention
i == 2 && return Vec(nil, (y - 1) / 2, nil)
i == 3 && return Vec((x + 1) / 2, nil, nil)
i == 4 && return Vec(nil, (y + 1) / 2, nil) # Changed sign, follow Ferrite's sign convention
i == 5 && return Vec((x - 1) / 2, nil, nil) # Changed sign, follow Ferrite's sign convention
i == 6 && return Vec(nil, nil, (z + 1) / 2)
throw(ArgumentError("no shape function $i for interpolation $ip"))
end

getnbasefunctions(::RaviartThomas{RefHexahedron, 1}) = 6
edgedof_indices(ip::RaviartThomas{RefHexahedron, 1}) = edgedof_interior_indices(ip)
facedof_indices(ip::RaviartThomas{RefHexahedron, 1}) = facedof_interior_indices(ip)
edgedof_interior_indices(::RaviartThomas{RefHexahedron, 1}) = ((), (), (), (), (), (), (), (), (), (), (), ())
facedof_interior_indices(::RaviartThomas{RefHexahedron, 1}) = ((1), (2), (3), (4), (5), (6))
adjust_dofs_during_distribution(::RaviartThomas{RefHexahedron, 1}) = false

function get_direction(::RaviartThomas{RefHexahedron, 1}, j, cell)
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not very elegant

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that is quite reasonable, but we could possibly abstract the logic to be more generally applicable.

Just realized that we have

Ferrite.jl/src/Grid/grid.jl

Lines 763 to 779 in bd65884

function OrientationInfo(path::NTuple{2, Int})
flipped = first(path) < last(path)
return OrientationInfo(flipped, 0)
end
function OrientationInfo(surface::NTuple{N, Int}) where {N}
min_idx = argmin(surface)
shift_index = min_idx - 1
if min_idx == 1
flipped = surface[2] < surface[end]
elseif min_idx == length(surface)
flipped = surface[1] < surface[end - 1]
else
flipped = surface[min_idx + 1] < surface[min_idx - 1]
end
return OrientationInfo(flipped, shift_index)
end

But that gives that "flipped" is positive according to the definition as I gave above - @termi-official, @AbdAlazezAhmed, @fredrikekre - is this intentional or did this use a different definition (essentially reversed)?

IMO we shouldn't use this function directly, but perhaps it would make more sense to do something like,

get_face_direction(cell, facenr) = get_face_direction(faces(cell)[facenr])

get_face_direction(facenodes::Tuple)
    min_idx = argmin(facenodes)
    if min_idx == 1
        positive = surface[2] < surface[end]
    elseif min_idx == length(surface)
        positive = surface[1] < surface[end - 1]
    else
        positive = surface[min_idx + 1] < surface[min_idx - 1]
    end
    return positive ? 1 : -1
end

or your code, which is equally good (except that I think argmin is cleaner (and potentially faster?) than findmin).

Then we should call get_face_direction in OrientationInfo.

We could do the same get_edge_direction abstraction as well, even if the logic is easier, maybe it makes the code easier to read to have this logic in the same place.

face = faces(cell)[j]
_, idx_min = findmin(face)

idx_next = mod(idx_min, length(face)) + 1
idx_prev = mod(idx_min - 2, length(face)) + 1
return face[idx_next] < face[idx_prev] ? 1 : -1
end

#####################################
# Brezzi-Douglas–Marini, H(div) #
#####################################
Expand Down Expand Up @@ -1942,3 +2000,59 @@ function get_direction(::Nedelec{RefTriangle, 2}, j, cell)
edge = edges(cell)[(j + 1) ÷ 2]
return ifelse(edge[2] > edge[1], 1, -1)
end

# RefTetrahedron, 1st order Lagrange
# https://defelement.org/elements/examples/tetrahedron-nedelec1-lagrange-1.html
function reference_shape_value(ip::Nedelec{RefTetrahedron, 1}, ξ::Vec{3,T}, i::Int) where {T}
x, y, z = ξ
nil = zero(T)

i == 1 && return Vec(1 - y - z, x, x)
i == 2 && return Vec(-y, x, nil)
i == 3 && return Vec(-y, x + z - 1, -y) # Changed sign, follow Ferrite's sign convention
i == 4 && return Vec(z, z, 1 - x - y)
i == 5 && return Vec(-z, nil, x)
i == 6 && return Vec(nil, -z, y)
throw(ArgumentError("no shape function $i for interpolation $ip"))
end

getnbasefunctions(::Nedelec{RefTetrahedron, 1}) = 6
edgedof_interior_indices(::Nedelec{RefTetrahedron, 1}) = ((1,), (2,), (3,), (4,), (5,), (6,))
facedof_indices(::Nedelec{RefTetrahedron, 1}) = ((1, 2, 3), (1, 4, 5), (2, 5, 6), (3, 4, 6))
adjust_dofs_during_distribution(::Nedelec{RefTetrahedron, 1}) = false

function get_direction(::Nedelec{RefTetrahedron, 1}, j, cell)
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is this expected to return?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It should return +1 if the edge associated with shape function j is positive, and -1 if negative. A positive edge is defined by the global node number of the first vertex being smaller than that of the second.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about shape functions associated with the cell faces (such as would be necessary for a Raviart-Thomas interpolation)?

Copy link
Member

@KnutAM KnutAM Feb 26, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we have defined what a positive face is in Ferrite, but maybe the definition that @fredrikekre suggested is already in place in some parts:
A face is positive if, starting from the lowest node number, the next vertex' node number is smaller than the previous vertex number.

Here, the vertices are defined by Ferrite.faces(cell), e.g. for cell = Tetrahedron((3,5,4,7)), we have
Ferrite.faces(Tetrahedron((3,5,4,2)) = ((3, 4, 5), (3, 5, 2), (5, 4, 2), (3, 2, 4))

Edit (was tired):

  • face 1: (3, 4, 5) - positive: 5-3-4 (5 > 4)
  • face 2: (3, 5, 2) - positive: 5-2-3 (5 > 3)
  • face 3: (5, 4, 2) - negative: 4-2-5 (4 < 5)
  • face 4: (3, 2, 4) - negative: 3-2-4 (3 < 4)

(Please correct me if I'm wrong with the convention @fredrikekre ! (Edit: Tag correct Fredrik :) )

So far, that was not required to use since we didn't have faces as facets for H(div) interpolations.

Copy link
Author

@gijswl gijswl Mar 1, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That appears to be right. Implementing this reasonsing (for RaviartThomas{RefTetrahedron,1} and RaviartThomas{RefHexahedron,1}), I get passing continuity tests.
I initially thought we should be able to default all the faces to positive, since the description of Ferrite.faces says the tuples it returns are oriented (normals pointing outward). That also appears to be the approach taken for the 2D elements.

edge = edges(cell)[j]
return ifelse(edge[2] > edge[1], 1, -1)
end

# RefHexahedron, 1st order Lagrange
# https://defelement.org/elements/examples/hexahedron-nedelec1-lagrange-1.html
function reference_shape_value(ip::Nedelec{RefHexahedron, 1}, ξ::Vec{3,T}, i::Int) where {T}
x, y, z = ξ
nil = zero(T)

i == 1 && return Vec((y * z - y - z + 1) / 8, nil, nil)
i == 2 && return Vec(nil, -(x * z - x + z - 1) / 8, nil)
i == 3 && return Vec((y * z - y + z - 1) / 8, nil, nil) # Changed sign, follow Ferrite's sign convention
i == 4 && return Vec(nil, -(x * z - x - z + 1) / 8, nil) # Changed sign, follow Ferrite's sign convention
i == 5 && return Vec(-(z * y - z + y - 1) / 8, nil, nil)
i == 6 && return Vec(nil, (x * z + x + z + 1) / 8, nil)
i == 7 && return Vec(-(y * z + y + z + 1) / 8, nil, nil) # Changed sign, follow Ferrite's sign convention
i == 8 && return Vec(nil, (z * x - z + x - 1) / 8, nil) # Changed sign, follow Ferrite's sign convention
i == 9 && return Vec(nil, nil, (x * y - x - y + 1) / 8)
i == 10 && return Vec(nil, nil, -(x * y - x + y - 1) / 8)
i == 11 && return Vec(nil, nil, (x * y + x + y + 1) / 8)
i == 12 && return Vec(nil, nil, -(y * x - y + x - 1) / 8)
throw(ArgumentError("no shape function $i for interpolation $ip"))
end

getnbasefunctions(::Nedelec{RefHexahedron, 1}) = 12
edgedof_interior_indices(::Nedelec{RefHexahedron, 1}) = ((1,), (2,), (3,), (4,), (5,), (6,), (7,), (8,), (9,), (10,), (11,), (12,))
facedof_indices(::Nedelec{RefHexahedron, 1}) = ((1, 2, 3, 4), (1, 5, 9, 10), (2, 6, 10, 11), (3, 7, 11, 12), (4, 8, 9, 12), (5, 6, 7, 8))
adjust_dofs_during_distribution(::Nedelec{RefHexahedron, 1}) = false

function get_direction(::Nedelec{RefHexahedron, 1}, j, cell)
edge = edges(cell)[j]
return ifelse(edge[2] > edge[1], 1, -1)
end
5 changes: 3 additions & 2 deletions test/test_continuity.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,8 +130,9 @@

test_ips = [
Lagrange{RefTriangle, 2}(), Lagrange{RefQuadrilateral, 2}(), Lagrange{RefHexahedron, 2}()^3, # Test should also work for identity mapping
Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(),
RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), BrezziDouglasMarini{RefTriangle, 1}(),
Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(), Nedelec{RefTetrahedron, 1}(), Nedelec{RefHexahedron, 1}(),
RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), RaviartThomas{RefTetrahedron, 1}(), RaviartThomas{RefHexahedron, 1}(),
BrezziDouglasMarini{RefTriangle, 1}(),
]

for ip in test_ips
Expand Down
46 changes: 42 additions & 4 deletions test/test_interpolations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -351,17 +351,28 @@ end
end

@testset "H(curl) and H(div)" begin
Hcurl_interpolations = [Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}()] # Nedelec{3, RefTetrahedron, 1}(), Nedelec{3, RefHexahedron, 1}()]
Hdiv_interpolations = [RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), BrezziDouglasMarini{RefTriangle, 1}()]
Hcurl_interpolations = [
Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(),
Nedelec{RefTetrahedron, 1}(), Nedelec{RefHexahedron, 1}()
]
Hdiv_interpolations = [
RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(),
RaviartThomas{RefTetrahedron, 1}(), RaviartThomas{RefHexahedron, 1}(),
BrezziDouglasMarini{RefTriangle, 1}()
]
test_interpolation_properties.(Hcurl_interpolations) # Requires PR1136
test_interpolation_properties.(Hdiv_interpolations) # Requires PR1136

# These reference moments define the functionals that an interpolation should fulfill
reference_moment(::RaviartThomas{RefTriangle, 1}, s, facet_shape_nr) = 1
reference_moment(::RaviartThomas{RefTriangle, 2}, s, facet_shape_nr) = facet_shape_nr == 1 ? (1 - s) : s
reference_moment(::RaviartThomas{RefTetrahedron, 1}, s0, s1, facet_shape_nr) = 1
reference_moment(::RaviartThomas{RefHexahedron, 1}, s0, s1, facet_shape_nr) = 1
reference_moment(::BrezziDouglasMarini{RefTriangle, 1}, s, facet_shape_nr) = facet_shape_nr == 1 ? (1 - s) : s
reference_moment(::Nedelec{RefTriangle, 1}, s, edge_shape_nr) = 1
reference_moment(::Nedelec{RefTriangle, 2}, s, edge_shape_nr) = edge_shape_nr == 1 ? (1 - s) : s
reference_moment(::Nedelec{RefTetrahedron, 1}, s, edge_shape_nr) = 1
reference_moment(::Nedelec{RefHexahedron, 1}, s, edge_shape_nr) = 1

function_space(::RaviartThomas) = Val(:Hdiv)
function_space(::BrezziDouglasMarini) = Val(:Hdiv)
Expand All @@ -373,7 +384,7 @@ end
return test_interpolation_functionals(function_space(ip), Val(Ferrite.getrefdim(ip)), ip)
end

# 2D, H(div) -> facet
# 2D, H(div)
function test_interpolation_functionals(::Val{:Hdiv}, ::Val{2}, ip::Interpolation)
RefShape = getrefshape(ip)
ipg = Lagrange{RefShape, 1}()
Expand All @@ -398,7 +409,34 @@ end
end
end

function test_interpolation_functionals(::Val{:Hcurl}, ::Val{2}, ip::Interpolation)
# 3D, H(div)
function test_interpolation_functionals(::Val{:Hdiv}, ::Val{3}, ip::Interpolation)
RefShape = getrefshape(ip)
ipg = Lagrange{RefShape, 1}()
for facetnr in 1:nfacets(RefShape)
facet_coords = getindex.((Ferrite.reference_coordinates(ipg),), Ferrite.reference_facets(RefShape)[facetnr])
dof_inds = Ferrite.facetdof_interior_indices(ip)[facetnr]
ξ(s0, s1) = facet_coords[1] + (facet_coords[2] - facet_coords[1]) * s0 + (facet_coords[3] - facet_coords[1]) * s1
weighted_normal = reference_normals(RefShape)[facetnr] * reference_face_area(ip, facetnr)
for (facet_shape_nr, shape_nr) in pairs(dof_inds)
moment_fun(s0, s1) = reference_moment(ip, s0, s1, facet_shape_nr)
f(s0, s1) = moment_fun(s0, s1) * (reference_shape_value(ip, ξ(s0, s1), shape_nr) ⋅ weighted_normal)

if(length(facet_coords) == 3)
val, _ = quadgk(s0 -> quadgk(s1 -> f(s0, s1), 0, 1 - s0; atol = 1.0e-8)[1], 0, 1; atol = 1.0e-8) # TODO Replace quadgk by more suitable 2D cubature
@test val ≈ 0.5
elseif(length(facet_coords) == 4)
val, _ = quadgk(s0 -> quadgk(s1 -> f(s0, s1), 0, 1; atol = 1.0e-8)[1], 0, 1; atol = 1.0e-8) # TODO Replace quadgk by more suitable 2D cubature
@test val ≈ 4
else
error("Cubature not defined for facets that are not triangles or quadrilaterals")
end
end
end
end

# 2D and 3D, H(curl)
function test_interpolation_functionals(::Val{:Hcurl}, ::Union{Val{2}, Val{3}}, ip::Interpolation)
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this a valid generalization, or do we need to implement additional tests for interpolation functionals with rdim > 2

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For $H(\mathrm{curl})$ (edge elements) I think this should be sufficient, but haven't validated - need to check what the actual functionals defined on DefElement says if the same can be used - but I don't see why not.

For $H(\text{div})$, we need two parameters in general and have to integrate a face, so havent' implemented how that should be done.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've made an attempt at validating Raviart-Thomas interpolation on tetrahedral and hexahedral elements, by computing the functional through integration over the faces. quadgk is not very well suited for that, so improvement is probably possible, but I couldn't find another method which allows the integration boundaries to depend on an integration variable (for the triangular faces).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could also use FacetValues I guess, so far I decided to avoid since to make the test purer and avoid using these before they have been tested.

Will look at the code later this week and give feedback!

Really nice work so far!

RefShape = getrefshape(ip)
ipg = Lagrange{RefShape, 1}()
for edgenr in 1:Ferrite.nedges(RefShape)
Expand Down