Skip to content

Commit

Permalink
use correct grid in interpolations
Browse files Browse the repository at this point in the history
  • Loading branch information
@chmerdon
chmerdon committed Feb 6, 2024
1 parent 366f9a4 commit e67a8d4
Show file tree
Hide file tree
Showing 22 changed files with 198 additions and 198 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "ExtendableFEMBase"
uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
authors = ["Christian Merdon <[email protected]>"]
version = "0.3.0"
version = "0.3.1"

[deps]
DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
Expand Down
6 changes: 3 additions & 3 deletions src/fedefs/h1_bubble.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,16 +35,16 @@ isdefined(FEType::Type{<:H1BUBBLE}, ::Type{<:Tetrahedron3D}) = true
interior_dofs_offset(::Type{ON_CELLS}, ::Type{H1BUBBLE{ncomponents}}, EG::Type{<:AbstractElementGeometry}) where {ncomponents} = 0

function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: H1BUBBLE, APT}
xCellVolumes = FE.xgrid[CellVolumes]
ncells = num_sources(FE.xgrid[CellNodes])
xCellVolumes = FE.dofgrid[CellVolumes]
ncells = num_sources(FE.dofgrid[CellNodes])
if items == []
items = 1:ncells
else
items = filter(!iszero, items)
end
ncomponents = get_ncomponents(FEType)
integrals4cell = zeros(T, ncomponents, ncells)
integrate!(integrals4cell, FE.xgrid, ON_CELLS, exact_function!; items = items, kwargs...)
integrate!(integrals4cell, FE.dofgrid, ON_CELLS, exact_function!; items = items, kwargs...)
for cell in items
if cell != 0
for c ∈ 1:ncomponents
Expand Down
14 changes: 7 additions & 7 deletions src/fedefs/h1_mini.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,35 +42,35 @@ get_ref_cellmoments(::Type{<:H1MINI}, ::Type{<:Quadrilateral2D}) = [1 // 4, 1 //
interior_dofs_offset(::Type{ON_CELLS}, ::Type{H1MINI{ncomponents, edim}}, EG::Type{<:AbstractElementGeometry}) where {ncomponents, edim} = num_nodes(EG)

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1MINI, APT}
nnodes = size(FE.xgrid[Coordinates], 2)
ncells = num_sources(FE.xgrid[CellNodes])
nnodes = size(FE.dofgrid[Coordinates], 2)
ncells = num_sources(FE.dofgrid[CellNodes])
point_evaluation!(Target, FE, AT_NODES, exact_function!; items = items, component_offset = nnodes + ncells, kwargs...)
end

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1MINI, APT}
# delegate edge nodes to node interpolation
subitems = slice(FE.xgrid[EdgeNodes], items)
subitems = slice(FE.dofgrid[EdgeNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
end

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1MINI, APT}
# delegate face nodes to node interpolation
subitems = slice(FE.xgrid[FaceNodes], items)
subitems = slice(FE.dofgrid[FaceNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
end

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1MINI, APT}
# delegate cell nodes to node interpolation
subitems = slice(FE.xgrid[CellNodes], items)
subitems = slice(FE.dofgrid[CellNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)

# fix cell bubble value by preserving integral mean
ensure_moments!(Target, FE, ON_CELLS, exact_function!; items = items, kwargs...)
end

function nodevalues!(Target::AbstractArray{<:Real, 2}, Source::AbstractArray{<:Real, 1}, FE::FESpace{<:H1MINI})
nnodes = num_sources(FE.xgrid[Coordinates])
ncells = num_sources(FE.xgrid[CellNodes])
nnodes = num_sources(FE.dofgrid[Coordinates])
ncells = num_sources(FE.dofgrid[CellNodes])
FEType = eltype(FE)
ncomponents = get_ncomponents(FEType)
offset4component = 0:(nnodes+ncells):ncomponents*(nnodes+ncells)
Expand Down
8 changes: 4 additions & 4 deletions src/fedefs/h1_p1.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,19 +34,19 @@ isdefined(FEType::Type{<:H1P1}, ::Type{<:Triangle2D}) = true
isdefined(FEType::Type{<:H1P1}, ::Type{<:Tetrahedron3D}) = true

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT}
nnodes = size(FE.xgrid[Coordinates], 2)
nnodes = size(FE.dofgrid[Coordinates], 2)
point_evaluation!(Target, FE, AT_NODES, exact_function; items = items, component_offset = nnodes, kwargs...)
end

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT}
# delegate edge nodes to node interpolation
subitems = slice(FE.xgrid[EdgeNodes], items)
subitems = slice(FE.dofgrid[EdgeNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function; items = subitems, kwargs...)
end

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT}
# delegate face nodes to node interpolation
subitems = slice(FE.xgrid[FaceNodes], items)
subitems = slice(FE.dofgrid[FaceNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function; items = subitems, kwargs...)
end

Expand All @@ -56,7 +56,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
point_evaluation_broken!(Target, FE, ON_CELLS, exact_function; items = items, time = time)
else
# delegate cell nodes to node interpolation
subitems = slice(FE.xgrid[CellNodes], items)
subitems = slice(FE.dofgrid[CellNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function; items = subitems, kwargs...)
end
end
Expand Down
22 changes: 11 additions & 11 deletions src/fedefs/h1_p2.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,13 +45,13 @@ interior_dofs_offset(::Type{<:AssemblyType}, ::Type{H1P2{ncomponents, edim}}, ::

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P2, APT}
edim = get_edim(FEType)
nnodes = size(FE.xgrid[Coordinates], 2)
nnodes = size(FE.dofgrid[Coordinates], 2)
if edim == 1
nedges = num_sources(FE.xgrid[CellNodes])
nedges = num_sources(FE.dofgrid[CellNodes])
elseif edim == 2
nedges = num_sources(FE.xgrid[FaceNodes])
nedges = num_sources(FE.dofgrid[FaceNodes])
elseif edim == 3
nedges = num_sources(FE.xgrid[EdgeNodes])
nedges = num_sources(FE.dofgrid[EdgeNodes])
end

point_evaluation!(Target, FE, AT_NODES, exact_function!; items = items, component_offset = nnodes + nedges, kwargs...)
Expand All @@ -61,7 +61,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
edim = get_edim(FEType)
if edim == 3
# delegate edge nodes to node interpolation
subitems = slice(FE.xgrid[EdgeNodes], items)
subitems = slice(FE.dofgrid[EdgeNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)

# perform edge mean interpolation
Expand All @@ -73,18 +73,18 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
edim = get_edim(FEType)
if edim == 2
# delegate face nodes to node interpolation
subitems = slice(FE.xgrid[FaceNodes], items)
subitems = slice(FE.dofgrid[FaceNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)

# perform face mean interpolation
ensure_moments!(Target, FE, ON_FACES, exact_function!; items = items, kwargs...)
elseif edim == 3
# delegate face edges to edge interpolation
subitems = slice(FE.xgrid[FaceEdges], items)
subitems = slice(FE.dofgrid[FaceEdges], items)
interpolate!(Target, FE, ON_EDGES, exact_function!; items = subitems, kwargs...)
elseif edim == 1
# delegate face nodes to node interpolation
subitems = slice(FE.xgrid[FaceNodes], items)
subitems = slice(FE.dofgrid[FaceNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
end
end
Expand All @@ -94,15 +94,15 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
edim = get_edim(FEType)
if edim == 2
# delegate cell faces to face interpolation
subitems = slice(FE.xgrid[CellFaces], items)
subitems = slice(FE.dofgrid[CellFaces], items)
interpolate!(Target, FE, ON_FACES, exact_function!; items = subitems, kwargs...)
elseif edim == 3
# delegate cell edges to edge interpolation
subitems = slice(FE.xgrid[CellEdges], items)
subitems = slice(FE.dofgrid[CellEdges], items)
interpolate!(Target, FE, ON_EDGES, exact_function!; items = subitems, kwargs...)
elseif edim == 1
# delegate cell nodes to node interpolation
subitems = slice(FE.xgrid[CellNodes], items)
subitems = slice(FE.dofgrid[CellNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)

# preserve cell integral
Expand Down
22 changes: 11 additions & 11 deletions src/fedefs/h1_p2b.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,12 +37,12 @@ get_ref_cellmoments(::Type{<:H1P2B}, ::Type{<:Triangle2D}) = [0 // 1, 0 // 1, 0

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P2B, APT}
edim = get_edim(FEType)
nnodes = size(FE.xgrid[Coordinates], 2)
offset = nnodes + num_sources(FE.xgrid[CellNodes])
nnodes = size(FE.dofgrid[Coordinates], 2)
offset = nnodes + num_sources(FE.dofgrid[CellNodes])
if edim == 2
offset += num_sources(FE.xgrid[FaceNodes])
offset += num_sources(FE.dofgrid[FaceNodes])
elseif edim == 3
offset += num_sources(FE.xgrid[EdgeNodes])
offset += num_sources(FE.dofgrid[EdgeNodes])
end

point_evaluation!(Target, FE, AT_NODES, exact_function!; items = items, component_offset = offset, kwargs...)
Expand All @@ -52,7 +52,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
edim = get_edim(FEType)
if edim == 3
# delegate edge nodes to node interpolation
subitems = slice(FE.xgrid[EdgeNodes], items)
subitems = slice(FE.dofgrid[EdgeNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)

# perform edge mean interpolation
Expand All @@ -64,21 +64,21 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
edim = get_edim(FEType)
if edim == 2
# delegate face nodes to node interpolation
subitems = slice(FE.xgrid[FaceNodes], items)
subitems = slice(FE.dofgrid[FaceNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)

# perform face mean interpolation
ensure_moments!(Target, FE, ON_FACES, exact_function!; items = items, kwargs...)
elseif edim == 3
# delegate face edges to edge interpolation
subitems = slice(FE.xgrid[FaceEdges], items)
subitems = slice(FE.dofgrid[FaceEdges], items)
interpolate!(Target, FE, ON_EDGES, exact_function!; items = subitems, kwargs...)

# perform face mean interpolation
# todo
elseif edim == 1
# delegate face nodes to node interpolation
subitems = slice(FE.xgrid[FaceNodes], items)
subitems = slice(FE.dofgrid[FaceNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
end
end
Expand All @@ -87,15 +87,15 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
edim = get_edim(FEType)
if edim == 2
# delegate cell faces to face interpolation
subitems = slice(FE.xgrid[CellFaces], items)
subitems = slice(FE.dofgrid[CellFaces], items)
interpolate!(Target, FE, ON_FACES, exact_function!; items = subitems, kwargs...)
elseif edim == 3
# delegate cell edges to edge interpolation
subitems = slice(FE.xgrid[CellEdges], items)
subitems = slice(FE.dofgrid[CellEdges], items)
interpolate!(Target, FE, ON_EDGES, exact_function!; items = subitems, kwargs...)
elseif edim == 1
# delegate cell nodes to node interpolation
subitems = slice(FE.xgrid[CellNodes], items)
subitems = slice(FE.dofgrid[CellNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
end

Expand Down
32 changes: 16 additions & 16 deletions src/fedefs/h1_p3.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,13 +49,13 @@ get_ref_cellmoments(::Type{<:H1P3}, ::Type{<:Triangle2D}) = [1 // 30, 1 // 30, 1

function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P3, APT}
edim = get_edim(FEType)
coffset = size(FE.xgrid[Coordinates], 2)
coffset = size(FE.dofgrid[Coordinates], 2)
if edim == 1
coffset += 2 * num_sources(FE.xgrid[CellNodes])
coffset += 2 * num_sources(FE.dofgrid[CellNodes])
elseif edim == 2
coffset += 2 * num_sources(FE.xgrid[FaceNodes]) + num_sources(FE.xgrid[CellNodes])
coffset += 2 * num_sources(FE.dofgrid[FaceNodes]) + num_sources(FE.dofgrid[CellNodes])
elseif edim == 3
coffset += 2 * num_sources(FE.xgrid[EdgeNodes]) + num_sources(FE.xgrid[FaceNodes])
coffset += 2 * num_sources(FE.dofgrid[EdgeNodes]) + num_sources(FE.dofgrid[FaceNodes])
end

point_evaluation!(Target, FE, AT_NODES, exact_function!; items = items, component_offset = coffset, kwargs...)
Expand All @@ -65,7 +65,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
edim = get_edim(FEType)
if edim == 3
# delegate edge nodes to node interpolation
subitems = slice(FE.xgrid[EdgeNodes], items)
subitems = slice(FE.dofgrid[EdgeNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)

# perform edge mean interpolation
Expand All @@ -77,21 +77,21 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
edim = get_edim(FEType)
if edim == 2
# delegate face nodes to node interpolation
subitems = slice(FE.xgrid[FaceNodes], items)
subitems = slice(FE.dofgrid[FaceNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)

# perform face mean interpolation
ensure_moments!(Target, FE, ON_FACES, exact_function!; items = items, order = 1, kwargs...)
elseif edim == 3
# delegate face edges to edge interpolation
subitems = slice(FE.xgrid[FaceEdges], items)
subitems = slice(FE.dofgrid[FaceEdges], items)
interpolate!(Target, FE, ON_EDGES, exact_function!; items = subitems, kwargs...)

# preserve face integral
ensure_moments!(Target, FE, ON_FACES, exact_function!; order = 0, items = items, kwargs...)
elseif edim == 1
# delegate face nodes to node interpolation
subitems = slice(FE.xgrid[FaceNodes], items)
subitems = slice(FE.dofgrid[FaceNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
end
end
Expand All @@ -101,18 +101,18 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
edim = get_edim(FEType)
if edim == 2
# delegate cell faces to face interpolation
subitems = slice(FE.xgrid[CellFaces], items)
subitems = slice(FE.dofgrid[CellFaces], items)
interpolate!(Target, FE, ON_FACES, exact_function!; items = subitems, kwargs...)

# fix cell bubble value by preserving integral mean
ensure_moments!(Target, FE, ON_CELLS, exact_function!; items = items, kwargs...)
elseif edim == 3
# delegate cell faces to face interpolation
subitems = slice(FE.xgrid[CellFaces], items)
subitems = slice(FE.dofgrid[CellFaces], items)
interpolate!(Target, FE, ON_FACES, exact_function!; items = subitems, kwargs...)
elseif edim == 1
# delegate cell nodes to node interpolation
subitems = slice(FE.xgrid[CellNodes], items)
subitems = slice(FE.dofgrid[CellNodes], items)
interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)

# preserve cell integral
Expand Down Expand Up @@ -198,7 +198,7 @@ end

# we need to change the ordering of the face dofs on faces that have a negative orientation sign
function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1P3{ncomponents, edim}, APT}, EG::Type{<:Triangle2D}) where {ncomponents, edim, Tv, Ti, APT}
xCellFaceSigns = FE.xgrid[CellFaceSigns]
xCellFaceSigns = FE.dofgrid[CellFaceSigns]
nfaces::Int = num_faces(EG)
function closure(subset_ids::Array{Int, 1}, cell)
for j ∈ 1:nfaces
Expand All @@ -220,7 +220,7 @@ end

# we need to change the ordering of the face dofs on faces that have a negative orientation sign
function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1P3{ncomponents, edim}, APT}, EG::Type{<:Tetrahedron3D}) where {ncomponents, edim, Tv, Ti, APT}
xCellEdgeSigns = FE.xgrid[CellEdgeSigns]
xCellEdgeSigns = FE.dofgrid[CellEdgeSigns]
nedges::Int = num_edges(EG)
function closure(subset_ids::Array{Int, 1}, cell)
for j ∈ 1:nedges
Expand All @@ -243,7 +243,7 @@ end

# we need to change the ordering of the face dofs on faces that have a negative orientation sign
function get_basissubset(::Type{ON_FACES}, FE::FESpace{Tv, Ti, H1P3{ncomponents, edim}, APT}, EG::Type{<:Triangle2D}) where {ncomponents, edim, Tv, Ti, APT}
xFaceEdgeSigns = FE.xgrid[FaceEdgeSigns]
xFaceEdgeSigns = FE.dofgrid[FaceEdgeSigns]
nedges::Int = num_edges(EG)
function closure(subset_ids::Array{Int, 1}, face)
for j ∈ 1:nedges
Expand All @@ -264,8 +264,8 @@ function get_basissubset(::Type{ON_FACES}, FE::FESpace{Tv, Ti, H1P3{ncomponents,
end

function get_basissubset(::Type{ON_BFACES}, FE::FESpace{Tv, Ti, H1P3{ncomponents, edim}, APT}, EG::Type{<:Triangle2D}) where {ncomponents, edim, Tv, Ti, APT}
xFaceEdgeSigns = FE.xgrid[FaceEdgeSigns]
xBFaceFaces = FE.xgrid[BFaceFaces]
xFaceEdgeSigns = FE.dofgrid[FaceEdgeSigns]
xBFaceFaces = FE.dofgrid[BFaceFaces]
nedges::Int = num_edges(EG)
function closure(subset_ids::Array{Int, 1}, bface)
for j ∈ 1:nedges
Expand Down
Loading

2 comments on commit e67a8d4

@chmerdon
Copy link
Member Author

Choose a reason for hiding this comment

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

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

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

Registration pull request created: JuliaRegistries/General/100349

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.3.1 -m "<description of version>" e67a8d4f3547a88125a95f8e9411656fdf32dbbc
git push origin v0.3.1

Please sign in to comment.