further extension of new interpolator structures (MomentInterpolator now also used for interior Hdiv dofs, new FluxEvaluator evaluates normal/tangential fluxes for Hdiv/Hcurl boundary dofs)

Author: chmerdon

src/fedefs/hcurl_n0.jl

@@ -38,41 +38,38 @@ isdefined(FEType::Type{<:HCURLN0}, ::Type{<:Triangle2D}) = true
 isdefined(FEType::Type{<:HCURLN0}, ::Type{<:Quadrilateral2D}) = true
 isdefined(FEType::Type{<:HCURLN0}, ::Type{<:Tetrahedron3D}) = true
 
-function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, data; items = [], kwargs...) where {T, Tv, Ti, FEType <: HCURLN0, APT}
+function N0_tangentflux_eval_2d!(result, f, qpinfo)
+    result[1] = -f[1] * qpinfo.normal[2] # rotated normal = tangent
+    return result[1] += f[2] * qpinfo.normal[1]
+end
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: HCURLN0{2}, APT} = FaceFluxEvaluator(N0_tangentflux_eval_2d!, FES, ON_FACES)
+
+function N0_tangentflux_eval_3d!(grid)
+    edgetangents = grid[EdgeTangents]
+    function closure(result, f, qpinfo)
+        result[1] = dot(f, view(edgetangents, :, qpinfo.item))
+    end
+    return closure
+end
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}) where {Tv, Ti, FEType <: HCURLN0{3}, APT} = FaceFluxEvaluator(N0_tangentflux_eval_3d!(FES.dofgrid), FES, ON_EDGES)
+
+function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HCURLN0, APT}
     edim = get_ncomponents(FEType)
     return if edim == 3
-        xEdgeTangents = FE.dofgrid[EdgeTangents]
-        if items == []
-            items = 1:size(xEdgeTangents, 2)
-        end
-        data_eval = zeros(T, 3)
-        function tangentflux_eval3d(result, qpinfo)
-            data(data_eval, qpinfo)
-            return result[1] = dot(data_eval, view(xEdgeTangents, :, qpinfo.item))
-        end
-        integrate!(Target, FE.dofgrid, ON_EDGES, tangentflux_eval3d; items = items, kwargs...)
+        get_interpolator(FE, ON_EDGES).evaluate!(Target, exact_function!, items; kwargs...)
     end
 end
 
-function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, data; items = [], kwargs...) where {T, Tv, Ti, FEType <: HCURLN0, APT}
+function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HCURLN0, APT}
     edim = get_ncomponents(FEType)
-    return if edim == 2
-        xFaceNormals = FE.dofgrid[FaceNormals]
-        if items == []
-            items = 1:size(xFaceNormals, 2)
-        end
-        data_eval = zeros(T, 2)
-        function tangentflux_eval2d(result, qpinfo)
-            data(data_eval, qpinfo)
-            result[1] = -data_eval[1] * xFaceNormals[2, qpinfo.item] # rotated normal = tangent
-            return result[1] += data_eval[2] * xFaceNormals[1, qpinfo.item]
-        end
-        integrate!(Target, FE.dofgrid, ON_FACES, tangentflux_eval2d; items = items, kwargs...)
+    if edim == 2
+        get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...)
     elseif edim == 3
         # delegate face edges to edge interpolation
         subitems = slice(FE.dofgrid[FaceEdges], items)
         interpolate!(Target, FE, ON_EDGES, data; items = subitems, kwargs...)
     end
+    return nothing
 end
 
 function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, data; items = [], kwargs...) where {Tv, Ti, FEType <: HCURLN0, APT}

src/fedefs/hcurl_n1.jl

@@ -27,30 +27,26 @@ get_dofmap_pattern(FEType::Type{<:HCURLN1{2}}, ::Union{Type{FaceDofs}, Type{BFac
 
 isdefined(FEType::Type{<:HCURLN1}, ::Type{<:Triangle2D}) = true
 
-function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, data; items = [], kwargs...) where {T, Tv, Ti, FEType <: HCURLN1, APT}
+function N1_tangentflux_eval_2d!(result, f, qpinfo)
+    result[1] = -f[1] * qpinfo.normal[2] # rotated normal = tangent
+    result[1] += f[2] * qpinfo.normal[1]
+    result[2] = result[1] * (qpinfo.xref[1] - 1 // 2)
+    return nothing
+end
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: HCURLN1{2}, APT} = FaceFluxEvaluator(N1_tangentflux_eval_2d!, FES, ON_FACES)
+
+
+function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HCURLN1, APT}
     edim = get_ncomponents(FEType)
     return if edim == 3
         # todo
     end
 end
 
-function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, data; items = [], kwargs...) where {T, Tv, Ti, FEType <: HCURLN1, APT}
+function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HCURLN1, APT}
     edim = get_ncomponents(FEType)
     return if edim == 2
-        xFaceNormals = FE.dofgrid[FaceNormals]
-        nfaces = num_sources(xFaceNormals)
-        if items == []
-            items = 1:size(xFaceNormals, 2)
-        end
-        # integrate normal flux of exact_function over edges
-        data_eval = zeros(T, 2)
-        function tangentflux_eval2d(result, qpinfo)
-            data(data_eval, qpinfo)
-            result[1] = -data_eval[1] * xFaceNormals[2, qpinfo.item] # rotated normal = tangent
-            result[1] += data_eval[2] * xFaceNormals[1, qpinfo.item]
-            return result[2] = result[1] * (qpinfo.xref[1] - 1 // 2)
-        end
-        integrate!(Target, FE.dofgrid, ON_FACES, tangentflux_eval2d; quadorder = 2, items = items, offset = [0, nfaces], kwargs...)
+        get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...)
     elseif edim == 3
         # delegate face edges to edge interpolation
         subitems = slice(FE.dofgrid[FaceEdges], items)

src/fedefs/hdiv_bdm1.jl

@@ -40,26 +40,20 @@ isdefined(FEType::Type{<:HDIVBDM1}, ::Type{<:Triangle2D}) = true
 isdefined(FEType::Type{<:HDIVBDM1}, ::Type{<:Quadrilateral2D}) = true
 isdefined(FEType::Type{<:HDIVBDM1}, ::Type{<:Tetrahedron3D}) = true
 
-
-function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, data; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVBDM1, APT}
-    ncomponents = get_ncomponents(FEType)
-    xFaceNormals = FE.dofgrid[FaceNormals]
-    nfaces = num_sources(xFaceNormals)
-    if items == []
-        items = 1:nfaces
-    end
-
-    # integrate normal flux of exact_function over edges
-    data_eval = zeros(T, ncomponents)
-    function normalflux_eval(result, qpinfo)
-        data(data_eval, qpinfo)
-        result[1] = dot(data_eval, view(xFaceNormals, :, qpinfo.item))
-        result[2] = result[1] * (qpinfo.xref[1] - 1 // ncomponents)
-        return if ncomponents == 3
-            result[3] = result[1] * (qpinfo.xref[2] - 1 // ncomponents)
+function BDM1_normalflux_eval!(dim)
+    function closure(result, f, qpinfo)
+        result[1] = dot(f, qpinfo.normal)
+        result[2] = result[1] * (qpinfo.xref[1] - 1 // dim)
+        if dim == 3
+            result[3] = result[1] * (qpinfo.xref[2] - 1 // dim)
         end
     end
-    return integrate!(Target, FE.dofgrid, ON_FACES, normalflux_eval; quadorder = 2, items = items, offset = 0:nfaces:((ncomponents - 1) * nfaces), kwargs...)
+end
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: HDIVBDM1, APT} = FaceFluxEvaluator(BDM1_normalflux_eval!(FEType.parameters[1]), FES, ON_FACES)
+
+
+function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVBDM1, APT}
+    get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
 function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, data; items = [], kwargs...) where {Tv, Ti, FEType <: HDIVBDM1, APT}

src/fedefs/hdiv_bdm2.jl

@@ -30,138 +30,30 @@ isdefined(FEType::Type{<:HDIVBDM2}, ::Type{<:Triangle2D}) = true
 interior_dofs_offset(::Type{<:ON_CELLS}, ::Type{<:HDIVBDM2{2}}, ::Type{<:Triangle2D}) = 9
 
 
-function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, data; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVBDM2, APT}
-    ncomponents = get_ncomponents(FEType)
-    xFaceNormals = FE.dofgrid[FaceNormals]
-    nfaces = num_sources(xFaceNormals)
-    if items == []
-        items = 1:nfaces
+function BDM2_normalflux_eval!(dim)
+    @assert dim == 2 "BDM3 for dim=3 not available yet"
+    function closure(result, f, qpinfo)
+        result[1] = dot(f, qpinfo.normal)
+        result[2] = result[1] * (qpinfo.xref[1] - 1 // dim)
+        result[3] = result[1] * (qpinfo.xref[1]^2 - qpinfo.xref[1] + 1 // 6)
     end
+end
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: HDIVBDM2, APT} = FaceFluxEvaluator(BDM2_normalflux_eval!(FEType.parameters[1]), FES, ON_FACES)
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: HDIVBDM2{2}, APT} = MomentInterpolator(FES, ON_CELLS; FEType_moments = H1MINI{1,2}, moments_dofs = [1,2,4], moments_operator = CurlScalar)
 
-    # integrate normal flux of exact_function over edges
-    data_eval = zeros(T, ncomponents)
-    function normalflux_eval(result, qpinfo)
-        data(data_eval, qpinfo)
-        result[1] = dot(data_eval, view(xFaceNormals, :, qpinfo.item))
-        result[2] = result[1] * (qpinfo.xref[1] - 1 // ncomponents)
-        return result[3] = result[1] * (qpinfo.xref[1]^2 - qpinfo.xref[1] + 1 // 6)
-    end
-    return integrate!(Target, FE.dofgrid, ON_FACES, normalflux_eval; quadorder = 4, items = items, offset = [0, nfaces, 2 * nfaces], kwargs...)
+
+function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVBDM2, APT}
+    get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
-function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, data; items = [], time = 0, bonus_quadorder = 0, kwargs...) where {T, Tv, Ti, FEType <: HDIVBDM2, APT}
+function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, exact_function!; items = [], time = 0, bonus_quadorder = 0, kwargs...) where {T, Tv, Ti, FEType <: HDIVBDM2, APT}
     # delegate cell faces to face interpolation
     subitems = slice(FE.dofgrid[CellFaces], items)
-    interpolate!(Target, FE, ON_FACES, data; items = subitems, kwargs...)
-
-    # set values of interior BDM2 functions as piecewise best-approximation
-    ncomponents = get_ncomponents(FEType)
-    EG = (ncomponents == 2) ? Triangle2D : Tetrahedron3D
-    ndofs = get_ndofs(ON_CELLS, FEType, EG)
-    interior_offset::Int = 9
-    nidofs::Int = ndofs - interior_offset
-    ncells = num_sources(FE.dofgrid[CellNodes])
-    xCellVolumes::Array{Tv, 1} = FE.dofgrid[CellVolumes]
-    xCellRegions = FE.dofgrid[CellRegions]
-    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; kwargs...)
-
-    # evaluation of gradient of P1 functions
-    FE3 = H1P1{1}
-    FES3 = FESpace{FE3, ON_CELLS}(FE.dofgrid)
-    FEBP1 = FEEvaluator(FES3, Gradient, qf; T = T)
-    # evaluation of curl of bubble functions
-    FE4 = H1BUBBLE{1}
-    FES4 = FESpace{FE4, ON_CELLS}(FE.dofgrid)
-    FEBB = FEEvaluator(FES4, CurlScalar, qf; T = T)
-
-    if items == []
-        items = 1:ncells
-    end
-
-    interiordofs = zeros(Int, nidofs)
-    basisvals::Array{T, 3} = FEB.cvals
-    basisvalsP1::Array{T, 3} = FEBP1.cvals
-    basisvalsB::Array{T, 3} = FEBB.cvals
-    IMM_face = zeros(T, nidofs, interior_offset)
-    IMM = zeros(T, nidofs, nidofs)
-    lb = zeros(T, nidofs)
-    temp::T = 0
-    data_eval = zeros(T, ncomponents)
-    for cell in items
-        # update basis
-        update_basis!(FEB, cell)
-        update_basis!(FEBP1, cell)
-        update_basis!(FEBB, cell)
-        fill!(IMM, 0)
-        fill!(IMM_face, 0)
-        fill!(lb, 0)
-
-        QP.item = cell
-        QP.cell = cell
-        QP.region = xCellRegions[cell]
-
-        # quadrature loop
-        for i in 1:length(qf.w)
-            # right-hand side : f times grad(P1),curl(bubble)
-            eval_trafo!(QP.x, FEB.L2G, FEB.xref[i])
-            QP.xref = FEB.xref[i]
-            data(data_eval, QP)
-            data_eval .*= xCellVolumes[cell] * qf.w[i]
-            for dof in 1:nidofs
-                for k in 1:ncomponents
-                    if dof < 3
-                        lb[dof] += data_eval[k] * basisvalsP1[k, dof, i]
-                    elseif dof == 3
-                        lb[dof] += data_eval[k] * basisvalsB[k, 1, i]
-                    end
-                end
-
-                # mass matrix of interior basis functions
-                for dof2 in 1:(nidofs - 1)
-                    temp = 0
-                    for k in 1:ncomponents
-                        temp += basisvals[k, interior_offset + dof, i] * basisvalsP1[k, dof2, i]
-                    end
-                    IMM[dof2, dof] += temp * xCellVolumes[cell] * qf.w[i]
-                end
-                temp = 0
-                for k in 1:ncomponents
-                    temp += basisvals[k, interior_offset + dof, i] * basisvalsB[k, 1, i]
-                end
-                IMM[3, dof] += temp * xCellVolumes[cell] * qf.w[i]
-
-                # mass matrix of face basis functions
-                for dof2 in 1:interior_offset
-                    temp = 0
-                    if dof < 3
-                        for k in 1:ncomponents
-                            temp += basisvalsP1[k, dof, i] * basisvals[k, dof2, i]
-                        end
-                    elseif dof == 3
-                        for k in 1:ncomponents
-                            temp += basisvalsB[k, 1, i] * basisvals[k, dof2, i]
-                        end
-                    end
-                    IMM_face[dof, dof2] += temp * xCellVolumes[cell] * qf.w[i]
-                end
-            end
-        end
-
-        # subtract face interpolation from right-hand side
-        for dof in 1:nidofs, dof2 in 1:interior_offset
-            lb[dof] -= Target[xCellDofs[dof2, cell]] * IMM_face[dof, dof2]
-        end
+    interpolate!(Target, FE, ON_FACES, exact_function!; items = subitems, kwargs...)
 
-        # solve local system
-        for dof in 1:nidofs
-            interiordofs[dof] = xCellDofs[interior_offset + dof, cell]
-        end
-        Target[interiordofs] = IMM \ lb
-    end
-    return
+    # set interior dofs such that moments with ∇P1 and curl(b) are preserved
+    # (b is the cell bubble)
+    get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
 ## only normalfluxes on faces

src/fedefs/hdiv_rt0.jl

@@ -35,25 +35,19 @@ isdefined(FEType::Type{<:HDIVRT0}, ::Type{<:Quadrilateral2D}) = true
 isdefined(FEType::Type{<:HDIVRT0}, ::Type{<:Tetrahedron3D}) = true
 isdefined(FEType::Type{<:HDIVRT0}, ::Type{<:Hexahedron3D}) = true
 
-function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, data; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVRT0, APT}
-    xFaceNormals = FE.dofgrid[FaceNormals]
-    if items == []
-        items = 1:size(xFaceNormals, 2)
-    end
+function RT0_normalflux_eval!(result, f, qpinfo)
+    return result[1] = dot(f, qpinfo.normal)
+end
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: HDIVRT0, APT} = FaceFluxEvaluator(RT0_normalflux_eval!, FES, ON_FACES)
 
-    # compute exact face means
-    data_eval = zeros(T, get_ncomponents(FEType))
-    function normalflux_eval(result, qpinfo)
-        data(data_eval, qpinfo)
-        return result[1] = dot(data_eval, view(xFaceNormals, :, qpinfo.item))
-    end
-    return integrate!(Target, FE.dofgrid, ON_FACES, normalflux_eval; items = items, kwargs...)
+function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVRT0, APT}
+    return get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
-function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, data; items = [], kwargs...) where {Tv, Ti, FEType <: HDIVRT0, APT}
+function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: HDIVRT0, APT}
     # delegate cell faces to face interpolation
     subitems = slice(FE.dofgrid[CellFaces], items)
-    return interpolate!(Target, FE, ON_FACES, data; items = subitems, kwargs...)
+    return interpolate!(Target, FE, ON_FACES, exact_function!; items = subitems, kwargs...)
 end
 
 # only normalfluxes on faces