diff --git a/CHANGELOG.md b/CHANGELOG.md
index 263c46d..1a72d5f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,7 +1,9 @@
# CHANGES
-## next version
+## v1.0.0
- doc improvements
+ - several bugfixes regarding finite elements defined on subgrids
+ - integrate now reacts on regions parameter and has proper docstrings
## v0.8.1 November 5, 2024
- fixed dimension in nodevals output (xdim was determined wrongly and caused huge arrays)
diff --git a/Project.toml b/Project.toml
index 7bdd0e8..55ecb83 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
name = "ExtendableFEMBase"
uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
authors = ["Christian Merdon <merdon@wias-berlin.de>"]
-version = "0.8.1"
+version = "1.0.0"
[deps]
DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
diff --git a/src/dofmaps.jl b/src/dofmaps.jl
index e14efcd..4704bd2 100644
--- a/src/dofmaps.jl
+++ b/src/dofmaps.jl
@@ -476,43 +476,70 @@ function init_dofmap_from_pattern!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Type{<
colstarts[end] = dofmap_totallength + 1
if FES.broken
- FES[DM] = SerialVariableTargetAdjacency(colstarts)
- return FES[DM]
- end
-
- colentries = zeros(Ti, dofmap_totallength)
- xItemDofs = VariableTargetAdjacency{Ti}(colentries, colstarts)
+ xItemDofs = SerialVariableTargetAdjacency(colstarts)
+ else
- cpattern::Array{DofMapPatternSegment, 1} = dofmap4EG[1].segments
- l::Int = 0
- offset::Int = 0
- pos::Int = 0
- q::Int = 0
- item_with_interiordofs::Int = 0
- for subitem in 1:nsubitems
- itemEG = xItemGeometries[subitem]
- item = sub2sup(subitem)
- if length(EG) > 1
- iEG = findfirst(isequal(itemEG), EG)
- end
- cpattern = dofmap4EG[iEG].segments
- l = length(cpattern)
- if has_interiordofs[iEG]
- item_with_interiordofs += 1
- end
- for c in 1:ncomponents
- offset = (c - 1) * offset4component
- for k in 1:l
- q = cpattern[k].ndofs
- if cpattern[k].type <: DofTypeNode && cpattern[k].each_component
- for n in 1:nnodes4EG[iEG]
+ colentries = zeros(Ti, dofmap_totallength)
+ xItemDofs = VariableTargetAdjacency{Ti}(colentries, colstarts)
+
+ cpattern::Array{DofMapPatternSegment, 1} = dofmap4EG[1].segments
+ l::Int = 0
+ offset::Int = 0
+ pos::Int = 0
+ q::Int = 0
+ item_with_interiordofs::Int = 0
+ for subitem in 1:nsubitems
+ itemEG = xItemGeometries[subitem]
+ item = sub2sup(subitem)
+ if length(EG) > 1
+ iEG = findfirst(isequal(itemEG), EG)
+ end
+ cpattern = dofmap4EG[iEG].segments
+ l = length(cpattern)
+ if has_interiordofs[iEG]
+ item_with_interiordofs += 1
+ end
+ for c in 1:ncomponents
+ offset = (c - 1) * offset4component
+ for k in 1:l
+ q = cpattern[k].ndofs
+ if cpattern[k].type <: DofTypeNode && cpattern[k].each_component
+ for n in 1:nnodes4EG[iEG]
+ for m in 1:q
+ pos += 1
+ xItemDofs.colentries[pos] = xItemNodes[n, item] + offset + (m - 1) * nnodes
+ end
+ end
+ offset += nnodes * q
+ elseif cpattern[k].type <: DofTypeFace && cpattern[k].each_component
+ for n in 1:nfaces4EG[iEG]
+ for m in 1:q
+ pos += 1
+ xItemDofs.colentries[pos] = xItemFaces[n, item] + offset + (m - 1) * nfaces
+ end
+ end
+ offset += nfaces * q
+ elseif cpattern[k].type <: DofTypeEdge && cpattern[k].each_component
+ for n in 1:nedges4EG[iEG]
+ for m in 1:q
+ pos += 1
+ xItemDofs.colentries[pos] = xItemEdges[n, item] + offset + (m - 1) * nedges
+ end
+ end
+ offset += nedges * q
+ elseif cpattern[k].type <: DofTypeInterior && cpattern[k].each_component
for m in 1:q
pos += 1
- xItemDofs.colentries[pos] = xItemNodes[n, item] + offset + (m - 1) * nnodes
+ xItemDofs.colentries[pos] = sub2sup(item_with_interiordofs) + offset
+ offset += nitems
end
end
- offset += nnodes * q
- elseif cpattern[k].type <: DofTypeFace && cpattern[k].each_component
+ end
+ end
+ offset = ncomponents * offset4component
+ for k in 1:l
+ q = cpattern[k].ndofs
+ if cpattern[k].type <: DofTypeFace && !cpattern[k].each_component
for n in 1:nfaces4EG[iEG]
for m in 1:q
pos += 1
@@ -520,7 +547,7 @@ function init_dofmap_from_pattern!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Type{<
end
end
offset += nfaces * q
- elseif cpattern[k].type <: DofTypeEdge && cpattern[k].each_component
+ elseif cpattern[k].type <: DofTypeEdge && !cpattern[k].each_component
for n in 1:nedges4EG[iEG]
for m in 1:q
pos += 1
@@ -528,7 +555,7 @@ function init_dofmap_from_pattern!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Type{<
end
end
offset += nedges * q
- elseif cpattern[k].type <: DofTypeInterior && cpattern[k].each_component
+ elseif cpattern[k].type <: DofTypeInterior && !cpattern[k].each_component
for m in 1:q
pos += 1
xItemDofs.colentries[pos] = sub2sup(item_with_interiordofs) + offset
@@ -537,37 +564,10 @@ function init_dofmap_from_pattern!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Type{<
end
end
end
- offset = ncomponents * offset4component
- for k in 1:l
- q = cpattern[k].ndofs
- if cpattern[k].type <: DofTypeFace && !cpattern[k].each_component
- for n in 1:nfaces4EG[iEG]
- for m in 1:q
- pos += 1
- xItemDofs.colentries[pos] = xItemFaces[n, item] + offset + (m - 1) * nfaces
- end
- end
- offset += nfaces * q
- elseif cpattern[k].type <: DofTypeEdge && !cpattern[k].each_component
- for n in 1:nedges4EG[iEG]
- for m in 1:q
- pos += 1
- xItemDofs.colentries[pos] = xItemEdges[n, item] + offset + (m - 1) * nedges
- end
- end
- offset += nedges * q
- elseif cpattern[k].type <: DofTypeInterior && !cpattern[k].each_component
- for m in 1:q
- pos += 1
- xItemDofs.colentries[pos] = sub2sup(item_with_interiordofs) + offset
- offset += nitems
- end
- end
- end
end
-
FES[DM] = xItemDofs
+
if FES.dofgrid !== FES.xgrid
## assume parent relation between xgrid and dofgrid
@assert FES.dofgrid[ParentGrid] == FES.xgrid "xgrid is not the parent grid of dofgrid !!!"
diff --git a/src/fedefs/h1_p3.jl b/src/fedefs/h1_p3.jl
index 4678f53..c60f854 100644
--- a/src/fedefs/h1_p3.jl
+++ b/src/fedefs/h1_p3.jl
@@ -201,8 +201,8 @@ function get_basis(::Type{<:AssemblyType}, ::Type{H1P3{ncomponents, edim}}, ::Ty
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.dofgrid[CellFaceSigns]
+function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1P3{ncomponents, edim}, APT}, EG::Type{<:Triangle2D}, xgrid) where {ncomponents, edim, Tv, Ti, APT}
+ xCellFaceSigns = xgrid[CellFaceSigns]
nfaces::Int = num_faces(EG)
return function closure(subset_ids::Array{Int, 1}, cell)
for j in 1:nfaces
@@ -223,8 +223,8 @@ function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1P3{ncomponents,
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.dofgrid[CellEdgeSigns]
+function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1P3{ncomponents, edim}, APT}, EG::Type{<:Tetrahedron3D}, xgrid) where {ncomponents, edim, Tv, Ti, APT}
+ xCellEdgeSigns = xgrid[CellEdgeSigns]
nedges::Int = num_edges(EG)
return function closure(subset_ids::Array{Int, 1}, cell)
for j in 1:nedges
diff --git a/src/fedefs/h1_pk.jl b/src/fedefs/h1_pk.jl
index 66fcd79..121695b 100644
--- a/src/fedefs/h1_pk.jl
+++ b/src/fedefs/h1_pk.jl
@@ -290,11 +290,11 @@ end
## if order > 2 on Triangl2D 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, H1Pk{ncomponents, edim, order}, APT}, EG::Type{<:Triangle2D}) where {ncomponents, edim, order, Tv, Ti, APT}
+function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1Pk{ncomponents, edim, order}, APT}, EG::Type{<:Triangle2D}, xgrid) where {ncomponents, edim, order, Tv, Ti, APT}
if order < 3
return NothingFunction # no reordering needed
end
- xCellFaceSigns = FE.dofgrid[CellFaceSigns]
+ xCellFaceSigns = xgrid[CellFaceSigns]
nfaces::Int = num_faces(EG)
ndofs_for_f::Int = order - 1
ndofs_for_c = get_ndofs(ON_CELLS, H1Pk{1, edim, order}, EG)
diff --git a/src/fedefs/h1v_br.jl b/src/fedefs/h1v_br.jl
index 9e8861c..9d6f39a 100644
--- a/src/fedefs/h1v_br.jl
+++ b/src/fedefs/h1v_br.jl
@@ -148,9 +148,9 @@ function get_basis(AT::Type{ON_CELLS}, FEType::Type{H1BR{2}}, EG::Type{<:Quadril
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{2}, APT}, ::Type{<:Triangle2D}) where {Tv, Ti, APT}
- xFaceNormals::Array{Tv, 2} = FE.dofgrid[FaceNormals]
- xCellFaces = FE.dofgrid[CellFaces]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{2}, APT}, ::Type{<:Triangle2D}, xgrid) where {Tv, Ti, APT}
+ xFaceNormals::Array{Tv, 2} = xgrid[FaceNormals]
+ xCellFaces = xgrid[CellFaces]
return function closure(coefficients::Array{<:Real, 2}, cell)
fill!(coefficients, 1.0)
coefficients[1, 7] = xFaceNormals[1, xCellFaces[1, cell]]
@@ -162,9 +162,9 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{2}, APT}, :
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{2}, APT}, ::Type{<:Quadrilateral2D}) where {Tv, Ti, APT}
- xFaceNormals::Array{Tv, 2} = FE.dofgrid[FaceNormals]
- xCellFaces = FE.dofgrid[CellFaces]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{2}, APT}, ::Type{<:Quadrilateral2D}, xgrid) where {Tv, Ti, APT}
+ xFaceNormals::Array{Tv, 2} = xgrid[FaceNormals]
+ xCellFaces = xgrid[CellFaces]
return function closure(coefficients::Array{<:Real, 2}, cell)
fill!(coefficients, 1.0)
coefficients[1, 9] = xFaceNormals[1, xCellFaces[1, cell]]
@@ -178,8 +178,8 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{2}, APT}, :
end
end
-function get_coefficients(::Type{<:ON_FACES}, FE::FESpace{Tv, Ti, H1BR{2}, APT}, ::Type{<:Edge1D}) where {Tv, Ti, APT}
- xFaceNormals::Array{Tv, 2} = FE.dofgrid[FaceNormals]
+function get_coefficients(::Type{<:ON_FACES}, FE::FESpace{Tv, Ti, H1BR{2}, APT}, ::Type{<:Edge1D}, xgrid) where {Tv, Ti, APT}
+ xFaceNormals::Array{Tv, 2} = xgrid[FaceNormals]
return function closure(coefficients::Array{<:Real, 2}, face)
# multiplication of face bubble with normal vector of face
fill!(coefficients, 1.0)
@@ -253,9 +253,9 @@ function get_basis(AT::Type{ON_CELLS}, ::Type{H1BR{3}}, EG::Type{<:Hexahedron3D}
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, ::Type{<:Tetrahedron3D}) where {Tv, Ti, APT}
- xFaceNormals::Array{Tv, 2} = FE.dofgrid[FaceNormals]
- xCellFaces::Adjacency{Ti} = FE.dofgrid[CellFaces]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, ::Type{<:Tetrahedron3D}, xgrid) where {Tv, Ti, APT}
+ xFaceNormals::Array{Tv, 2} = xgrid[FaceNormals]
+ xCellFaces::Adjacency{Ti} = xgrid[CellFaces]
return function closure(coefficients::Array{<:Real, 2}, cell)
# multiplication with normal vectors
fill!(coefficients, 1.0)
@@ -276,8 +276,8 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, :
end
-function get_coefficients(::Type{<:ON_FACES}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, ::Type{<:Triangle2D}) where {Tv, Ti, APT}
- xFaceNormals::Array{Tv, 2} = FE.dofgrid[FaceNormals]
+function get_coefficients(::Type{<:ON_FACES}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, ::Type{<:Triangle2D}, xgrid) where {Tv, Ti, APT}
+ xFaceNormals::Array{Tv, 2} = xgrid[FaceNormals]
return function closure(coefficients::Array{<:Real, 2}, face)
# multiplication of face bubble with normal vector of face
fill!(coefficients, 1.0)
@@ -287,9 +287,9 @@ function get_coefficients(::Type{<:ON_FACES}, FE::FESpace{Tv, Ti, H1BR{3}, APT},
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, ::Type{<:Hexahedron3D}) where {Tv, Ti, APT}
- xFaceNormals::Array{Tv, 2} = FE.dofgrid[FaceNormals]
- xCellFaces::Adjacency{Ti} = FE.dofgrid[CellFaces]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, ::Type{<:Hexahedron3D}, xgrid) where {Tv, Ti, APT}
+ xFaceNormals::Array{Tv, 2} = xgrid[FaceNormals]
+ xCellFaces::Adjacency{Ti} = xgrid[CellFaces]
return function closure(coefficients::Array{<:Real, 2}, cell)
# multiplication with normal vectors
fill!(coefficients, 1.0)
@@ -316,8 +316,8 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, :
end
-function get_coefficients(::Type{<:ON_FACES}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, ::Type{<:Quadrilateral2D}) where {Tv, Ti, APT}
- xFaceNormals::Array{Tv, 2} = FE.dofgrid[FaceNormals]
+function get_coefficients(::Type{<:ON_FACES}, FE::FESpace{Tv, Ti, H1BR{3}, APT}, ::Type{<:Quadrilateral2D}, xgrid) where {Tv, Ti, APT}
+ xFaceNormals::Array{Tv, 2} = xgrid[FaceNormals]
return function closure(coefficients::Array{<:Real, 2}, face)
# multiplication of face bubble with normal vector of face
fill!(coefficients, 1.0)
diff --git a/src/fedefs/h1v_p1teb.jl b/src/fedefs/h1v_p1teb.jl
index 9657628..504cf7f 100644
--- a/src/fedefs/h1v_p1teb.jl
+++ b/src/fedefs/h1v_p1teb.jl
@@ -253,9 +253,9 @@ function get_basis(AT::Type{ON_CELLS}, ::Type{H1P1TEB{3}}, EG::Type{<:Tetrahedro
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1P1TEB{3}, APT}, ::Type{<:Tetrahedron3D}) where {Tv, Ti, APT}
- xEdgeTangents::Array{Tv, 2} = FE.dofgrid[EdgeTangents]
- xCellEdges = FE.dofgrid[CellEdges]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1P1TEB{3}, APT}, ::Type{<:Tetrahedron3D}, xgrid) where {Tv, Ti, APT}
+ xEdgeTangents::Array{Tv, 2} = xgrid[EdgeTangents]
+ xCellEdges = xgrid[CellEdges]
return function closure(coefficients::Array{<:Real, 2}, cell)
fill!(coefficients, 1.0)
for e in 1:6, k in 1:2
@@ -265,9 +265,9 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, H1P1TEB{3}, APT}
end
end
-function get_coefficients(::Type{<:ON_FACES}, FE::FESpace{Tv, Ti, H1P1TEB{3}, APT}, ::Type{<:Triangle2D}) where {Tv, Ti, APT}
- xEdgeTangents::Array{Tv, 2} = FE.dofgrid[EdgeTangents]
- xFaceEdges = FE.dofgrid[FaceEdges]
+function get_coefficients(::Type{<:ON_FACES}, FE::FESpace{Tv, Ti, H1P1TEB{3}, APT}, ::Type{<:Triangle2D}, xgrid) where {Tv, Ti, APT}
+ xEdgeTangents::Array{Tv, 2} = xgrid[EdgeTangents]
+ xFaceEdges = xgrid[FaceEdges]
return function closure(coefficients::Array{<:Real, 2}, face)
fill!(coefficients, 1.0)
for e in 1:3, k in 1:2
diff --git a/src/fedefs/hcurl_n0.jl b/src/fedefs/hcurl_n0.jl
index 266bc66..6e5253e 100644
--- a/src/fedefs/hcurl_n0.jl
+++ b/src/fedefs/hcurl_n0.jl
@@ -144,8 +144,8 @@ function get_basis(::Type{ON_CELLS}, ::Type{HCURLN0{3}}, ::Type{<:Tetrahedron3D}
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HCURLN0, APT}, EG::Type{<:AbstractElementGeometry2D}) where {Tv, Ti, APT}
- xCellFaceSigns = FE.dofgrid[CellFaceSigns]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HCURLN0, APT}, EG::Type{<:AbstractElementGeometry2D}, xgrid) where {Tv, Ti, APT}
+ xCellFaceSigns = xgrid[CellFaceSigns]
nfaces = num_faces(EG)
return function closure(coefficients, cell)
# multiplication with normal vector signs
@@ -156,8 +156,8 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HCURLN0, APT},
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HCURLN0, APT}, EG::Type{<:AbstractElementGeometry3D}) where {Tv, Ti, APT}
- xCellEdgeSigns = FE.dofgrid[CellEdgeSigns]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HCURLN0, APT}, EG::Type{<:AbstractElementGeometry3D}, xgrid) where {Tv, Ti, APT}
+ xCellEdgeSigns = xgrid[CellEdgeSigns]
nedges = num_edges(EG)
return function closure(coefficients, cell)
# multiplication with normal vector signs
diff --git a/src/fedefs/hcurl_n1.jl b/src/fedefs/hcurl_n1.jl
index 739637e..39d7238 100644
--- a/src/fedefs/hcurl_n1.jl
+++ b/src/fedefs/hcurl_n1.jl
@@ -104,8 +104,8 @@ function get_basis(::Type{ON_CELLS}, ::Type{HCURLN1{2}}, ::Type{<:Triangle2D})
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HCURLN1, APT}, EG::Type{<:AbstractElementGeometry2D}) where {Tv, Ti, APT}
- xCellFaceSigns = FE.dofgrid[CellFaceSigns]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HCURLN1, APT}, EG::Type{<:AbstractElementGeometry2D}, xgrid) where {Tv, Ti, APT}
+ xCellFaceSigns = xgrid[CellFaceSigns]
nfaces = num_faces(EG)
return function closure(coefficients, cell)
# multiplication with normal vector signs (only RT0)
diff --git a/src/fedefs/hdiv_bdm1.jl b/src/fedefs/hdiv_bdm1.jl
index b653d67..f9eb12d 100644
--- a/src/fedefs/hdiv_bdm1.jl
+++ b/src/fedefs/hdiv_bdm1.jl
@@ -203,8 +203,8 @@ function get_basis(::Type{ON_CELLS}, ::Type{HDIVBDM1{3}}, ::Type{<:Tetrahedron3D
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM1, APT}, EG::Type{<:AbstractElementGeometry2D}) where {Tv, Ti, APT}
- xCellFaceSigns::Union{VariableTargetAdjacency{Int32}, Array{Int32, 2}} = FE.dofgrid[CellFaceSigns]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM1, APT}, EG::Type{<:AbstractElementGeometry2D}, xgrid) where {Tv, Ti, APT}
+ xCellFaceSigns::Union{VariableTargetAdjacency{Int32}, Array{Int32, 2}} = xgrid[CellFaceSigns]
nfaces::Int = num_faces(EG)
dim::Int = dim_element(EG)
return function closure(coefficients::Array{<:Real, 2}, cell::Int)
@@ -218,9 +218,8 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM1, APT}
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM1, APT}, EG::Type{<:AbstractElementGeometry3D}) where {Tv, Ti, APT}
- xCellFaceSigns::Union{VariableTargetAdjacency{Int32}, Array{Int32, 2}} = FE.dofgrid[CellFaceSigns]
- xCellFaceOrientations::Union{VariableTargetAdjacency{Int32}, Array{Int32, 2}} = FE.dofgrid[CellFaceOrientations]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM1, APT}, EG::Type{<:AbstractElementGeometry3D}, xgrid) where {Tv, Ti, APT}
+ xCellFaceSigns::Union{VariableTargetAdjacency{Int32}, Array{Int32, 2}} = xgrid[CellFaceSigns]
nfaces::Int = num_faces(EG)
dim::Int = dim_element(EG)
return function closure(coefficients::Array{<:Real, 2}, cell::Int)
@@ -238,8 +237,8 @@ end
# the RT0 and those two BDM1 face functions are chosen
# such that they reflect the two moments with respect to the second and third node
# of the global face enumeration
-function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM1, APT}, EG::Type{<:AbstractElementGeometry3D}) where {Tv, Ti, APT}
- xCellFaceOrientations = FE.dofgrid[CellFaceOrientations]
+function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM1, APT}, EG::Type{<:AbstractElementGeometry3D}, xgrid) where {Tv, Ti, APT}
+ xCellFaceOrientations = xgrid[CellFaceOrientations]
nfaces::Int = num_faces(EG)
orientation = xCellFaceOrientations[1, 1]
shift4orientation1::Array{Int, 1} = [1, 0, 1, 2]
diff --git a/src/fedefs/hdiv_bdm2.jl b/src/fedefs/hdiv_bdm2.jl
index 6746b2f..635b716 100644
--- a/src/fedefs/hdiv_bdm2.jl
+++ b/src/fedefs/hdiv_bdm2.jl
@@ -66,7 +66,7 @@ function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{T
xCellDofs::DofMapTypes{Ti} = FE[CellDofs]
qf = QuadratureRule{T, EG}(max(4, 2 + bonus_quadorder))
FEB = FEEvaluator(FE, Identity, qf; T = T)
- QP = QPInfos(FE.dofgrid)
+ QP = QPInfos(FE.dofgrid; kwargs...)
# evaluation of gradient of P1 functions
FE3 = H1P1{1}
@@ -206,8 +206,8 @@ function get_basis(::Type{ON_CELLS}, ::Type{HDIVBDM2{2}}, ::Type{<:Triangle2D})
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM2, APT}, EG::Type{<:AbstractElementGeometry2D}) where {Tv, Ti, APT}
- xCellFaceSigns::Union{VariableTargetAdjacency{Int32}, Array{Int32, 2}} = FE.dofgrid[CellFaceSigns]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM2, APT}, EG::Type{<:AbstractElementGeometry2D}, xgrid) where {Tv, Ti, APT}
+ xCellFaceSigns::Union{VariableTargetAdjacency{Int32}, Array{Int32, 2}} = xgrid[CellFaceSigns]
nfaces::Int = num_faces(EG)
dim::Int = dim_element(EG)
return function closure(coefficients::Array{<:Real, 2}, cell::Int)
diff --git a/src/fedefs/hdiv_rt0.jl b/src/fedefs/hdiv_rt0.jl
index 77dae04..073d2ab 100644
--- a/src/fedefs/hdiv_rt0.jl
+++ b/src/fedefs/hdiv_rt0.jl
@@ -117,8 +117,8 @@ function get_basis(::Type{ON_CELLS}, ::Type{HDIVRT0{3}}, ::Type{<:Hexahedron3D})
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRT0, APT}, EG::Type{<:AbstractElementGeometry}) where {Tv, Ti, APT}
- xCellFaceSigns = FE.dofgrid[CellFaceSigns]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRT0, APT}, EG::Type{<:AbstractElementGeometry}, xgrid) where {Tv, Ti, APT}
+ xCellFaceSigns = xgrid[CellFaceSigns]
nfaces = num_faces(EG)
return function closure(coefficients, cell)
# multiplication with normal vector signs
diff --git a/src/fedefs/hdiv_rt1.jl b/src/fedefs/hdiv_rt1.jl
index f9f2eba..c12b795 100644
--- a/src/fedefs/hdiv_rt1.jl
+++ b/src/fedefs/hdiv_rt1.jl
@@ -198,8 +198,8 @@ function get_basis(::Type{ON_CELLS}, ::Type{HDIVRT1{3}}, ::Type{<:Tetrahedron3D}
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRT1{2}, APT}, EG::Type{<:Triangle2D}) where {Tv, Ti, APT}
- xCellFaceSigns = FE.dofgrid[CellFaceSigns]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRT1{2}, APT}, EG::Type{<:Triangle2D}, xgrid) where {Tv, Ti, APT}
+ xCellFaceSigns = xgrid[CellFaceSigns]
nfaces = num_faces(EG)
return function closure(coefficients, cell)
fill!(coefficients, 1.0)
@@ -211,8 +211,8 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRT1{2}, AP
end
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRT1{3}, APT}, EG::Type{<:Tetrahedron3D}) where {Tv, Ti, APT}
- xCellFaceSigns = FE.dofgrid[CellFaceSigns]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRT1{3}, APT}, EG::Type{<:Tetrahedron3D}, xgrid) where {Tv, Ti, APT}
+ xCellFaceSigns = xgrid[CellFaceSigns]
nfaces = num_faces(EG)
return function closure(coefficients, cell)
fill!(coefficients, 1.0)
@@ -232,8 +232,8 @@ end
# the RT0 and those two BDM1 face functions are chosen
# such that they reflect the two moments with respect to the second and third node
# of the global face enumeration
-function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRT1{3}, APT}, EG::Type{<:Tetrahedron3D}) where {Tv, Ti, APT}
- xCellFaceOrientations = FE.dofgrid[CellFaceOrientations]
+function get_basissubset(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRT1{3}, APT}, EG::Type{<:Tetrahedron3D}, xgrid) where {Tv, Ti, APT}
+ xCellFaceOrientations = xgrid[CellFaceOrientations]
nfaces::Int = num_faces(EG)
orientation = xCellFaceOrientations[1, 1]
shift4orientation1::Array{Int, 1} = [1, 0, 1, 2]
diff --git a/src/fedefs/hdiv_rtk.jl b/src/fedefs/hdiv_rtk.jl
index 5f11a0b..be4ac7c 100644
--- a/src/fedefs/hdiv_rtk.jl
+++ b/src/fedefs/hdiv_rtk.jl
@@ -76,7 +76,7 @@ function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{T
xCellDofs::DofMapTypes{Ti} = FE[CellDofs]
qf = QuadratureRule{T, EG}(max(2 * order, order + 1 + bonus_quadorder))
FEB = FEEvaluator(FE, Identity, qf; T = T)
- QP = QPInfos(FE.dofgrid)
+ QP = QPInfos(FE.dofgrid; kwargs...)
# evaluation of gradient of P1 functions
FE3 = order == 1 ? L2P0{ncomponents} : H1Pk{ncomponents, 2, order - 1}
@@ -208,8 +208,8 @@ function get_basis(::Type{ON_CELLS}, ::Type{HDIVRTk{2, order}}, ::Type{<:Triangl
end
-function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRTk{2, order}, APT}, EG::Type{<:Triangle2D}) where {Tv, Ti, APT, order}
- xCellFaceSigns::Union{VariableTargetAdjacency{Int32}, Array{Int32, 2}} = FE.dofgrid[CellFaceSigns]
+function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRTk{2, order}, APT}, EG::Type{<:Triangle2D}, xgrid) where {Tv, Ti, APT, order}
+ xCellFaceSigns::Union{VariableTargetAdjacency{Int32}, Array{Int32, 2}} = xgrid[CellFaceSigns]
nfaces::Int = num_faces(EG)
dim::Int = dim_element(EG)
return function closure(coefficients::Array{<:Real, 2}, cell::Int)
diff --git a/src/feevaluator.jl b/src/feevaluator.jl
index 8eed8e6..a92a8fc 100644
--- a/src/feevaluator.jl
+++ b/src/feevaluator.jl
@@ -43,14 +43,14 @@ Note that matrix-valued operators evaluations, e.g. for Gradient, are given as a
function FEEvaluator(
FE::FESpace{TvG, TiG, FEType, FEAPT},
operator::Type{<:StandardFunctionOperator},
- qrule::QuadratureRule{TvR, EG};
+ qrule::QuadratureRule{TvR, EG},
+ xgrid = FE.dofgrid; # FE.xgrid also reasonable in certain situations
L2G = nothing,
T = Float64,
AT = ON_CELLS
) where {TvG, TiG, TvR, FEType <: AbstractFiniteElement, EG <: AbstractElementGeometry, FEAPT <: AssemblyType}
xref = qrule.xref
- xgrid = FE.xgrid
if L2G === nothing
L2G = L2GTransformer(EG, xgrid, AT)
end
@@ -76,7 +76,7 @@ function FEEvaluator(
end
end
edim = max(1, dim_element(EG))
- xdim = size(FE.xgrid[Coordinates], 1)
+ xdim = size(xgrid[Coordinates], 1)
resultdim = Int(Length4Operator(operator, edim, ncomponents))
# evaluate basis on reference domain
@@ -90,11 +90,11 @@ function FEEvaluator(
coeff_handler = NothingFunction
if FEType <: Union{AbstractH1FiniteElementWithCoefficients, AbstractHdivFiniteElement, AbstractHcurlFiniteElement}
coefficients = ones(T, ncomponents, ndofs4item)
- coeff_handler = get_coefficients(FEAT, FE, EG)
+ coeff_handler = get_coefficients(FEAT, FE, EG, xgrid)
end
# set subset handler (only relevant if ndofs4item_all > ndofs4item)
- subset_handler = get_basissubset(FEAT, FE, EG)
+ subset_handler = get_basissubset(FEAT, FE, EG, xgrid)
# prepare offsets and additional coefficients
offsets = 0:edim:((ncomponents - 1) * edim) # edim steps
@@ -104,8 +104,8 @@ function FEEvaluator(
@warn "ndofs = 0 for FEType = $FEType on EG = $EG"
offsets2 = []
end
- compressiontargets = _prepare_compressiontargets(operator, FE.xgrid, AT, edim)
- coefficients_op = _prepare_additional_coefficients(operator, FE.xgrid, AT, edim)
+ compressiontargets = _prepare_compressiontargets(operator, xgrid, AT, edim)
+ coefficients_op = _prepare_additional_coefficients(operator, xgrid, AT, edim)
current_eval = zeros(T, resultdim, ndofs4item, length(xref))
derivorder = NeededDerivative4Operator(operator)
@@ -178,7 +178,7 @@ _prepare_additional_coefficients(::Type{<:AbstractFunctionOperator}, xgrid, AT,
_prepare_additional_coefficients(::Type{<:NormalFlux}, xgrid, AT, edim) = AT == ON_BFACES ? view(xgrid[FaceNormals], :, xgrid[BFaceFaces]) : xgrid[FaceNormals]
_prepare_additional_coefficients(::Type{<:TangentialGradient}, xgrid, AT, edim) = xgrid[FaceNormals]
function _prepare_additional_coefficients(::Type{<:TangentFlux}, xgrid, AT, edim)
- if edim == 2
+ if edim == 1
return _prepare_additional_coefficients(NormalFlux, xgrid, AT, edim)
else
return AT == ON_BEDGES ? view(xgrid[EdgeTangents], :, xgrid[BEdgeEdges]) : xgrid[EdgeTangents]
diff --git a/src/feevaluator_h1.jl b/src/feevaluator_h1.jl
index f601720..ae6c7be 100644
--- a/src/feevaluator_h1.jl
+++ b/src/feevaluator_h1.jl
@@ -158,6 +158,21 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Norm
return nothing
end
+
+# TANGENTFLUX H1 (ON_FACES, ON_BFACES)
+function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:TangentFlux, <:AbstractH1FiniteElement})
+ # fetch normal of item
+ normal = view(FEBE.coefficients_op, :, FEBE.citem[])
+ subset = _update_subset!(FEBE)
+ cvals = FEBE.cvals
+ refbasisvals = FEBE.refbasisvals
+ fill!(cvals, 0)
+ for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), k in 1:length(normal)
+ cvals[1, dof_i, i] = refbasisvals[i][subset[dof_i], 1] * normal[2] - refbasisvals[i][subset[dof_i], 2] * normal[1]
+ end
+ return nothing
+end
+
# NORMALFLUX H1+COEFFICIENTS (ON_FACES, ON_BFACES)
function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:NormalFlux, <:AbstractH1FiniteElementWithCoefficients})
# fetch normal of item
diff --git a/src/finiteelements.jl b/src/finiteelements.jl
index 866c06d..99284e7 100644
--- a/src/finiteelements.jl
+++ b/src/finiteelements.jl
@@ -243,7 +243,7 @@ $(TYPEDSIGNATURES)
returns the coefficients for local evaluations of finite element functions
( see e.g. h1v_br.jl for a use-case)
"""
-function get_coefficients(::Type{<:AssemblyType}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{<:AbstractElementGeometry}) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
+function get_coefficients(::Type{<:AssemblyType}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{<:AbstractElementGeometry}, xgrid = FE.dofgrid) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
return NothingFunction
end
@@ -262,7 +262,7 @@ where different basis functions are chosen
depending on the face orientations (which in 3D is not just a sign)
"""
-function get_basissubset(::Type{<:AssemblyType}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{<:AbstractElementGeometry}) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
+function get_basissubset(::Type{<:AssemblyType}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{<:AbstractElementGeometry}, xgrid = FE.dofgrid) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
return NothingFunction
end
@@ -451,8 +451,8 @@ include("fedefs/hcurl_n0.jl");
include("fedefs/hcurl_n1.jl");
-function get_coefficients(::Type{ON_BFACES}, FE::FESpace{Tv, Ti, FEType, APT}, EG::Type{<:AbstractElementGeometry}) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
- get_coeffs_on_face = get_coefficients(ON_FACES, FE, EG)
+function get_coefficients(::Type{ON_BFACES}, FE::FESpace{Tv, Ti, FEType, APT}, EG::Type{<:AbstractElementGeometry}, xgrid) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
+ get_coeffs_on_face = get_coefficients(ON_FACES, FE, EG, xgrid)
xBFaceFaces = FE.xgrid[BFaceFaces]
return function closure(coefficients, bface)
get_coeffs_on_face(coefficients, xBFaceFaces[bface])
@@ -471,12 +471,12 @@ function get_reconstruction_matrix(T::Type{<:Real}, FE::FESpace, FER::FESpace)
ncells = num_sources(xgrid[CellNodes])
rhandlers = [get_reconstruction_coefficient(ON_CELLS, FE, FER, EG[1])]
- chandlers = [get_coefficients(ON_CELLS, FER, EG[1])]
+ chandlers = [get_coefficients(ON_CELLS, FER, EG[1], xgrid)]
shandlers = [get_basissubset(ON_CELLS, FER, EG[1])]
for j in 2:length(EG)
append!(rhandlers, [get_reconstruction_coefficients(ON_CELLS, FE, FER, EG[j])])
- append!(chandlers, [get_coefficients(ON_CELLS, FER, EG[j])])
- append!(shandlers, [get_basissubset(ON_CELLS, FER, EG[j])])
+ append!(chandlers, [get_coefficients(ON_CELLS, FER, EG[j], xgrid)])
+ append!(shandlers, [get_basissubset(ON_CELLS, FER, EG[j],)])
end
ndofs_FE = zeros(Int, length(EG))
diff --git a/src/quadrature.jl b/src/quadrature.jl
index 4e55504..22e771e 100644
--- a/src/quadrature.jl
+++ b/src/quadrature.jl
@@ -645,7 +645,22 @@ end
"""
$(TYPEDSIGNATURES)
-Integration that writes result on every item into integral4items.
+Integration of an integrand of the signature
+
+ integrand!(result, qpinfo)
+
+over the entities AT of the grid. The result on every item is written into integral4items
+(at least of size nitems x length of result).
+As usual, qpinfo allows to access information
+at the current quadrature point.
+
+Keyword arguments:
+- quadorder = specifies the quadrature order (default = 0)
+- regions = restricts integration to these regions (default = [] meaning all regions)
+- items = restricts integration to these item numbers (default = [] meaning all items)
+- time = sets the time to be passed to qpinfo (default = 0)
+- params = sets the params array to be passed to qpinfo (default = [])
+- offset = offset for writing into result array integral4items (default = [0]),
"""
function integrate!(
integral4items::AbstractArray{T},
@@ -655,6 +670,7 @@ function integrate!(
offset = [0],
bonus_quadorder = 0,
quadorder = 0,
+ regions = [],
time = 0,
items = [],
force_quadrature_rule = nothing,
@@ -702,6 +718,14 @@ function integrate!(
local2global[j] = L2GTransformer(EG[j], grid, AT)
end
+ ## prepare regions
+ visit_region = zeros(Bool, maximum(xItemRegions))
+ if length(regions) > 0
+ visit_region[regions] .= true
+ else
+ visit_region .= true
+ end
+
# loop over items
if items == []
items = 1:nitems
@@ -734,6 +758,16 @@ function integrate!(
fill!(integral4items, 0)
for item::Int in items
+ if xItemRegions[item] > 0
+ if !(visit_region[xItemRegions[item]]) || AT == ON_IFACES
+ continue
+ end
+ else
+ if length(regions) > 0
+ continue
+ end
+ end
+
# find index for CellType
if ngeoms > 1
itemET = xItemGeometries[item]
@@ -764,6 +798,16 @@ function integrate!(
fill!(integral4items, 0)
for item::Int in items
+ if xItemRegions[item] > 0
+ if !(visit_region[xItemRegions[item]]) || AT == ON_IFACES
+ continue
+ end
+ else
+ if length(regions) > 0
+ continue
+ end
+ end
+
# find index for CellType
if ngeoms > 1
itemET = xItemGeometries[item]
@@ -779,7 +823,20 @@ end
"""
$(TYPEDSIGNATURES)
-Integration that returns total integral.
+Integration that returns total integral of an integrand of the signature
+
+ integrand!(result, qpinfo)
+
+over the entities AT of the grid. Here, qpinfo allows to access information
+at the current quadrature point.
+The length of the result needs to be specified via resultdim.
+
+Keyword arguments:
+- quadorder = specifies the quadrature order (default = 0)
+- regions = restricts integration to these regions (default = [] meaning all regions)
+- items = restricts integration to these item numbers (default = [] meaning all items)
+- time = sets the time to be passed to qpinfo (default = 0)
+- params = sets the params array to be passed to qpinfo (default = [])
"""
function integrate(
grid::ExtendableGrid,
diff --git a/src/reconstructionoperators.jl b/src/reconstructionoperators.jl
index 4c5d07a..4270253 100644
--- a/src/reconstructionoperators.jl
+++ b/src/reconstructionoperators.jl
@@ -31,7 +31,8 @@ end
function FEEvaluator(
FE::FESpace{TvG, TiG, FEType, FEAPT},
operator::Type{<:Reconstruct{FETypeReconst, stdop}},
- qrule::QuadratureRule{TvR, EG};
+ qrule::QuadratureRule{TvR, EG},
+ xgrid = FE.xgrid;
L2G = nothing,
T = Float64,
AT = ON_CELLS
@@ -40,7 +41,6 @@ function FEEvaluator(
@debug "Creating FEBasisEvaluator for reconstruction of $stdop operator of $FEType into $FETypeReconst on $EG"
## generate FESpace for reconstruction
- xgrid = FE.xgrid
FE2 = FESpace{FETypeReconst}(xgrid)
## collect dimension information
@@ -62,7 +62,7 @@ function FEEvaluator(
if L2G === nothing
L2G = L2GTransformer(EG, xgrid, AT)
end
- FEB = FEEvaluator(FE2, stdop, qrule; L2G = L2G, T = T, AT = AT)
+ FEB = FEEvaluator(FE2, stdop, qrule, xgrid; L2G = L2G, T = T, AT = AT)
## reconstruction coefficient handler
reconst_handler = ReconstructionHandler(FE, FE2, AT, EG)