Better operator evaluations + more operator tests

CHANGELOG.md

@@ -1,5 +1,10 @@
 # CHANGES
 
+## v1.7.0 May 26, 2026
+- improved operator evaluations, more operator tests
+- improved Examples, fixed Example210
+- scalarplot for FEVectorPlot now uses its name as default title (like vectorplot)
+
 ## v1.6.1 May 18, 2026
 - fixed bug in `PointEvaluator` initialization which broke automatic differentiation of the `evaluate!` function with respect to spatial coordinates
 

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
-version = "1.6.1"
+version = "1.7.0"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]

examples/Example200_LowLevelPoisson.jl

@@ -63,10 +63,10 @@ function main(;
         println(stdout, barplot(["Grid", "FaceNodes", "celldofs", "Assembly", "Solve"], [time_grid, time_facenodes, time_dofmap, time_assembly, time_solve], title = "Runtimes"))
 
         # plot
-        scalarplot!(plt[1, 1], solution[1])
-        scalarplot!(plt[1, 2], solution[1], Gradient; abs = true, clear = true)
-        vectorplot!(plt[1, 2], solution[1], Gradient; clear = false)
-        gridplot!(plt[1, 3], xgrid; markersize = 0)
+        scalarplot!(plt[1, 1], solution["u"])
+        scalarplot!(plt[1, 2], solution["u"], Gradient; abs = true, clear = true)
+        vectorplot!(plt[1, 2], solution["u"], Gradient; clear = false)
+        gridplot!(plt[1, 3], xgrid; markersize = 0, title = "grid")
         reveal(plt)
     end
 
@@ -75,7 +75,7 @@ end
 
 
 function solve_poisson_lowlevel(fe_space, mu, rhs)
-    solution = FEVector(fe_space)
+    solution = FEVector(fe_space; tags = ["u"])
     stiffness_matrix = FEMatrix(fe_space, fe_space)
     rhs_vector = FEVector(fe_space)
     println("Assembling...")
@@ -100,7 +100,7 @@ function solve_poisson_lowlevel(fe_space, mu, rhs)
     return solution, time_assembly, time_solve
 end
 
-function assemble!(A::ExtendableSparseMatrix, b::Vector, fe_space, rhs, mu = 1)
+function assemble!(A::ExtendableSparseMatrix, b::Vector, fe_space::FESpace{Tv, Ti}, rhs, mu = 1) where {Tv, Ti}
     xgrid = fe_space.xgrid
     EG = xgrid[UniqueCellGeometries][1]
     FEType = eltype(fe_space)
@@ -130,7 +130,7 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, fe_space, rhs, mu = 1)
         dof_j::Int, dof_k::Int = 0, 0
         x::Vector{Float64} = zeros(Float64, 2)
 
-        return loop_allocations += @allocated for cell in 1:ncells
+        return loop_allocations += @allocated for cell::Ti in 1:ncells
             ## update FE basis evaluators
             FEBasis_∇.citem[] = cell
             update_basis!(FEBasis_∇)

examples/Example205_LowLevelSpaceTimePoisson.jl

@@ -179,7 +179,7 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES_time, FES_space, f,
 
     ## ASSEMBLY LOOP
     loop_allocations = 0
-    function barrier(EG_time, EG_space, L2G_time::L2GTransformer, L2G_space::L2GTransformer)
+    function barrier(EG_time, EG_space, L2G_time::L2GTransformer{Tv, Ti}, L2G_space::L2GTransformer{Tv, Ti}) where {Tv, Ti}
         ## barrier function to avoid allocations by type dispatch
 
         ndofs4cell_time::Int = get_ndofs(ON_CELLS, FEType_time, EG_time)
@@ -192,7 +192,7 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES_time, FES_space, f,
         t::Vector{Float64} = zeros(Float64, 1)
 
         ## assemble Laplacian
-        loop_allocations += @allocated for cell in 1:ncells_space
+        loop_allocations += @allocated for cell::Ti in 1:ncells_space
             ## update FE basis evaluators for space
             FEBasis_space_∇.citem[] = cell
             update_basis!(FEBasis_space_∇)
@@ -208,7 +208,7 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES_time, FES_space, f,
             Aloc .*= cellvolumes_space[cell]
 
             ## add local matrix to global matrix
-            for time_cell in 1:ncells_time
+            for time_cell::Ti in 1:ncells_time
                 update_trafo!(L2G_time, time_cell)
                 for jT in 1:ndofs4cell_time, kT in 1:ndofs4cell_time
                     dofTj = celldofs_time[jT, time_cell]
@@ -237,7 +237,7 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES_time, FES_space, f,
             for qp in 1:nweights_space
                 ## evaluate coordinates of quadrature point in space
                 eval_trafo!(x, L2G_space, xref_space[qp])
-                for time_cell in 1:ncells_time
+                for time_cell::Ti in 1:ncells_time
                     update_trafo!(L2G_time, time_cell)
                     for qpT in 1:nweights_time
                         ## evaluate time coordinate

examples/Example210_LowLevelNavierStokes.jl

@@ -87,22 +87,20 @@ function main(; nref = 5, teval = 0, order = 2, Plotter = UnicodePlots)
     sol = solve_stokes_lowlevel(FES; teval = teval)
 
     ## move integral mean of pressure
-    pmean = sum(compute_error(sol[2], nothing, order, 1))
-    for j in 1:sol[2].FES.ndofs
-        sol[2][j] -= pmean
-    end
+    pmean = sum(compute_error(sol["p"], nothing, order, 1))
+    view(sol["p"]) .-= pmean
 
     ## calculate l2 error
-    error_u = sqrt(sum(compute_error(sol[1], u!, 2)))
-    error_p = sqrt(sum(compute_error(sol[2], p!, 2)))
-    println("\nl2 error velo = $(error_u)")
+    error_u = sqrt(sum(compute_error(sol["u"], u!, 2)))
+    error_p = sqrt(sum(compute_error(sol["p"], p!, 2)))
+    println("\nl2 error velocity = $(error_u)")
     println("l2 error pressure = $(error_p)")
 
     ## plot
     plt = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, resolution = (1200, 600))
-    scalarplot!(plt[1, 1], sol[1]; title = "|u| + quiver", abs = true)
-    vectorplot!(plt[1, 1], sol[1]; clear = false)
-    scalarplot!(plt[1, 2], sol[2]; title = "p")
+    scalarplot!(plt[1, 1], sol["u"]; abs = true)
+    vectorplot!(plt[1, 1], sol["u"]; clear = false)
+    scalarplot!(plt[1, 2], sol["p"])
     reveal(plt)
 
     return sol, plt
@@ -120,12 +118,12 @@ function compute_error(uh::FEVectorBlock, u, order = get_polynomialorder(get_FET
     uhval = zeros(Float64, ncomponents)
     uval = zeros(Float64, ncomponents)
     L2G = L2GTransformer(EG, xgrid, ON_CELLS)
-    QP = QPInfos(xgrid)
+    QP = QPInfos(xgrid; time = 0)
     qf = VertexRule(EG, order)
     FEB = FEEvaluator(FES, Identity, qf)
 
-    function barrier(L2G::L2GTransformer)
-        for cell in 1:num_cells(xgrid)
+    function barrier(L2G::L2GTransformer{Tv, Ti}) where {Tv, Ti}
+        for cell::Ti in 1:num_cells(xgrid)
             update_trafo!(L2G, cell)
             update_basis!(FEB, cell)
             for (qp, weight) in enumerate(qf.w)
@@ -156,7 +154,7 @@ end
 function solve_stokes_lowlevel(FES; teval = 0)
 
     println("Initializing system...")
-    sol = FEVector(FES)
+    sol = FEVector(FES; tags = ["u", "p"])
     A = FEMatrix(FES)
     b = FEVector(FES)
     @time update_system! = prepare_assembly!(A, b, FES[1], FES[2], sol)
@@ -172,6 +170,13 @@ function solve_stokes_lowlevel(FES; teval = 0)
         push!(fixed_dofs, FES[1].ndofs + 1) ## fix one pressure dof
     end
 
+    ## fix boundary dofs
+    for dof in fixed_dofs
+        A.entries[dof, dof] = 1.0e60
+        b.entries[dof] = 1.0e60 * u_init.entries[dof]
+    end
+    ExtendableSparse.flush!(A.entries)
+
     for it in 1:20
         ## solve
         println("\nITERATION $it\n=============")
@@ -262,7 +267,7 @@ function prepare_assembly!(A, b, FESu, FESp, sol; teval = 0)
     value = DiffResults.value(Dresult)
 
     ## ASSEMBLY LOOP
-    function barrier(EG, L2G::L2GTransformer, linear::Bool, nonlinear::Bool)
+    function barrier(EG, L2G::L2GTransformer{Tv, Ti}, linear::Bool, nonlinear::Bool) where {Tv, Ti}
         ## barrier function to avoid allocations caused by L2G
 
         ndofs4cell_u::Int = get_ndofs(ON_CELLS, FEType_u, EG)
@@ -273,7 +278,7 @@ function prepare_assembly!(A, b, FESu, FESp, sol; teval = 0)
         fval::Vector{Float64} = zeros(Float64, 2)
         x::Vector{Float64} = zeros(Float64, 2)
 
-        for cell in 1:ncells
+        for cell::Ti in 1:ncells
             ## update FE basis evaluators
             update_basis!(FEBasis_∇u, cell)
             update_basis!(FEBasis_idu, cell)

src/fedefs/h1v_br.jl

@@ -312,7 +312,7 @@ end
 function get_reconstruction_coefficients!(xgrid::ExtendableGrid{Tv, Ti}, ::Union{Type{<:ON_FACES}, Type{<:ON_BFACES}}, FE::Type{<:H1BR{2}}, FER::Type{<:HDIVRT0{2}}, ::Type{<:Edge1D}) where {Tv, Ti}
     xFaceVolumes::Array{<:Real, 1} = xgrid[FaceVolumes]
     xFaceNormals::Array{<:Real, 2} = xgrid[FaceNormals]
-    return function closure(coefficients::Array{<:Real, 2}, face::Int)
+    return function closure(coefficients::Array{<:Real, 2}, face)
         coefficients[1, 1] = 1 // 2 * xFaceVolumes[face] * xFaceNormals[1, face]
         coefficients[2, 1] = 1 // 2 * xFaceVolumes[face] * xFaceNormals[1, face]
         coefficients[3, 1] = 1 // 2 * xFaceVolumes[face] * xFaceNormals[2, face]
@@ -350,7 +350,7 @@ function get_reconstruction_coefficients!(xgrid::ExtendableGrid{Tv, Ti}, ::Union
     xFaceVolumes::Array{Tv, 1} = xgrid[FaceVolumes]
     xFaceNormals::Array{Tv, 2} = xgrid[FaceNormals]
     nfacenodes::Int = num_nodes(EG)
-    return function closure(coefficients::Array{<:Real, 2}, face::Int)
+    return function closure(coefficients::Array{<:Real, 2}, face)
         for j in 1:nfacenodes, k in 1:3
             coefficients[(j - 1) * 3 + j, 1] = 1 // nfacenodes * xFaceVolumes[face] * xFaceNormals[k, face]
             coefficients[(j - 1) * 3 + j, 2] = -1 // 36 * xFaceVolumes[face] * xFaceNormals[k, face]

src/fedefs/hdiv_bdm1.jl

@@ -201,7 +201,7 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM1, 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)
+    return function closure(coefficients::Array{<:Real, 2}, cell)
         fill!(coefficients, 1.0)
         # multiplication with normal vector signs (only RT0)
         for j in 1:nfaces, k in 1:dim
@@ -216,7 +216,7 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM1, 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)
+    return function closure(coefficients::Array{<:Real, 2}, cell)
         fill!(coefficients, 1.0)
         for j in 1:nfaces, k in 1:dim
             coefficients[k, 3 * j - 2] = xCellFaceSigns[j, cell] # RT0

src/fedefs/hdiv_bdm2.jl

@@ -102,7 +102,7 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVBDM2, 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)
+    return function closure(coefficients::Array{<:Real, 2}, cell)
         fill!(coefficients, 1.0)
         # multiplication with normal vector signs (only RT0)
         for j in 1:nfaces, k in 1:dim

src/fedefs/hdiv_rtk.jl

@@ -121,7 +121,7 @@ function get_coefficients(::Type{ON_CELLS}, FE::FESpace{Tv, Ti, <:HDIVRTk{2, ord
     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)
+    return function closure(coefficients::Array{<:Real, 2}, cell)
         fill!(coefficients, 1.0)
         # multiplication with normal vector signs (only RT0, RT2, RT4 etc.)
         for j in 1:nfaces, k in 1:dim

src/feevaluator.jl

@@ -1,7 +1,7 @@
 abstract type FEEvaluator{T <: Real, TvG <: Real, TiG <: Integer} end
 
 struct SingleFEEvaluator{T <: Real, TvG <: Real, TiG <: Integer, operator, FEType, EG, FType_basis <: Function, FType_coeffs <: Function, FType_subset <: Function, FType_jac <: Function} <: FEEvaluator{T, TvG, TiG}
-    citem::Base.RefValue{Int}                   # current item
+    citem::Base.RefValue{TiG}                   # current item
     FE::FESpace{TvG, TiG, FEType}                 # link to full FE (e.g. for coefficients)
     L2G::L2GTransformer{TvG, TiG, EG}           # local2global mapper
     L2GAinv::Array{TvG, 2}                       # 2nd heap for transformation matrix (e.g. Piola + mapderiv)
@@ -189,7 +189,7 @@ function FEEvaluator(
     end
 
     return SingleFEEvaluator{T, TvG, TiG, operator, FEType, EG, typeof(refbasis), typeof(coeff_handler), typeof(subset_handler), typeof(jacobian_wrap)}(
-        Ref(0),
+        Ref(TiG(0)),
         FE,
         L2G,
         L2GAinv,